I am trying to figure out what the following function $f:\Bbb{R}^3-\{\mathbf{0}\}\to\Bbb{R}$ defined below (in pseudocode) does:
function $f(\mathbf{v}\in \Bbb{R}^3-\{\mathbf{0}\})$
{
$a=b=0;$ $(a,b\in \Bbb{R})$
for ${(i=0; i<I; i=i+1)}$
{
$\mathbf{v}=\frac{|\mathbf{v}|}{||\mathbf{v}||^2}-\mathbf{A};$
$a=a+|||\mathbf{v}||-b|;$
$b=||\mathbf{v}||;$
}
return $ a;$
}
Where $I$ is some positive integer that determines the number of iterations the loop will run, and $\mathbf{A}$ is some vector in the set $\{(t,t,t)\in\Bbb{R}^3|t\in[-1,1]\}$.
(The notation $|\mathbf{v}|$ denotes component wise absolute value of the vector elements)
I do not know how this function looks like and what it does and I am only interested in vectors $\mathbf{v}\in[-1,1]^3-\{\mathbf{0}\}$.
Also what is $lim_{\mathbf{v}\to \mathbf{0}}f(v)$?
Here are some picks that were generated by my shader (maybe they can help):
When $I = 1$ and $\mathbf{A}$ is any vector in the set $\{(t,t,t)\in\Bbb{R}^3|t\in[-1,1]\}$ the images will be something like this

As we increase the value of $I$ things get interesting and we get imaged like those:
When $I = 10$ and $\mathbf{A}=\mathbf{(0.5, 0.5, 0.5)}$

When $I = 10$ and $\mathbf{A}=\mathbf{(0.63, 0.63, 0.63)}$

When $I = 10$ and $\mathbf{A}=\mathbf{(0.74, 0.74, 0.74)}$

When $I = 10$ and $\mathbf{A}=\mathbf{(1, 1, 1)}$

Since the center of every image is solid white I can conjecture that $lim_{\mathbf{v}\to\mathbf{0}}f(\mathbf{v})=\infty$
the main part of my GLSL fragment shader that produces the images and implements the function $f$ defind above is
vec3 color = vec3(0.);
vec3 v = vec3(2 * (TexCoord.xy - .5), 0.); // TexCoord.xy in [0,1]^2
float a = 0.;
float b = 0.;
// The loop that I do not understand
for (uint i = 0U; i < I; i++)
{
v = abs(v) / dot(v,v) - A;
a += abs(length(v) - b);
b = length(v);
}
a = a * a * a; // Adding some contrast
float brightness = 0.015;
color = vec3(.1, .01, .0001) * a * brightness;
FragColor = vec4(color, 1.);
Any help will be appreciated.
(I've encountered this function while trying to achieve some interesting fractal like visual effects while doing graphics programming)
EDIT
In all pictures, the parameter $Z$ denotes the $Z$ coordinate of the vector $\mathbf{v}$
$I=90$, $A=0.5$ and $Z=0$
$I=90$, $A=1$ and $Z=0$
$I=100$, $A=0$ and $Z=0$
$I=100$, $A=0.5$ and $Z=0$
$I=100$, $A=1$ and $Z=0$

$I=100$, $A= 0.042373 $ and $Z=0$
$I=100$, $A= 0.906780 $ and $Z=0$
$I=100$, $A= 0.932203 $ and $Z=0$

$I=100$, $A= 0.046610 $ and $Z= 0.894068 $ (cool dotted effect)

$I=16$, $A= 0.194915 $ and $Z=0$
$I=29$, $A= 0.194915 $ and $Z=0$

$I=500$, $A=0$ and $Z=0$
$I=500$, $A=0.5$ and $Z=0$ (nearly all white)
$I=500$, $A=1$ and $Z=0$

$I=1000$, $A=0.5$ and $Z=0$ (nearly all white)

Here are some images that view the first iteration of the loop on the xy plane with $A=(0.5,0.5,0)$:
(I will try to generate images for greater iterations)








I feel the need to take this to a more mathematical realm. So here's a commentary to long for comments.
Set-Up
We have a system of recurrence relations defined by,
$$(1) \quad \vec v_{n+1}={{|\vec v_n|} \over {||\vec v_n||^2}}-\vec A$$ $$(2) \quad a_{n+1}=a_n+| \ ||\vec v_n||-b_n|$$ $$(3) \quad b_{n+1}=||\vec v_n||$$
First let's reduce this to a system of two relations. $(3)$ is redundant, we can simply use $||v_{n+1}||$. We now have,
$$(1) \quad \vec v_{n+1}={{|\vec v_n|} \over {||\vec v_n||^2}}-\vec A$$ $$(4) \quad a_{n+1}=a_n+| \ ||\vec v_n||-||\vec v_{n-1}|| \ |$$
$(1)$ happens to be independent of $(4)$ so I'd suggest that we understand $(1)$ before trying to figure out $(4)$, which is just the picture generator.
Interpretations
It may be hard to interpret what $(1)$ and $(4)$ are actually doing. Here's a briefer. $(4)$ is easy. It measures how volatile the path taken from $v_0$ to $v_n$ was. A higher value means the path was longer and more volatile. A lower value means it either didn't have a value explosion, or that it converged to some value. The interpretation of $(1)$ is more subtle. I still don't know exactly what it does, but it seems to normalize the vector to a unit vector, then it does that again, then it subtracts a constant and we start over. Because the vector is normalized, most of the interesting behavior happens in the unit circle around the origin.
This system is Chaotic
We can start by constructing a toy model. We'll do this by collapsing $(1)$ to only a recurrence in $1$-d.
$$(5) \quad v_{n+1}={1 \over {|v_n|}}-c$$
Which should retain some of the same characteristics. You seem interested in the case where $v \rightarrow 0$. Using the toy model $(5)$ we see that,
$$v_0=0$$ $$v_1=\infty$$ $$v_2=-c$$ $$v_3={1 \over c}-c$$ $$v_4={1 \over {|{1 \over c}-c|}}-c$$ And so on...
The general behavior of the iterations is very similar to continued fractions. The only difference is the absolute value sign. This changes the behavior dramatically. For instance here's the behavior of $v_0=0.7$ and $c=0.5$ under iteration.
The behavior is chaotic. So basically, this toy model shows you why your getting incredibly complex behavior in the $3$-d model. It's the inclusion of the absolute value that's doing it! Compare this to the version without it here.
Note that it has an explicit formula for the evolution.
Singularity Analysis
Let's move to the $2$-d case, equivalent to setting all of the $z$ components to $0$. I'll skip the algebra, but here's the evolution of the vector system.
$$(6) \quad x_{n+1}={{|x_n|} \over {x_n^2+y_n^2}}-v_x$$ $$(7) \quad y_{n+1}={{|y_n|} \over {x_n^2+y_n^2}}-v_y$$
Once again, we'll ignore the detailed analysis of the path indicator function. (I hope to show why that is later). By inspection we can see that the function will become infinite whenever a vector passes through infinity. By inspection we can see that this happens if the initial components satisfy,
$$x_0^2+y_0^2=0$$
So it really wouldn't be difficult to prove that for the $3$-d case with $z_0$, the same equation must be satisfied. So this answers your question. $\lim_{\vec v \to 0} a_n=\infty$. So the path indicator function becomes infinite. That's why the origin is white in your pictures. The path indicator function also becomes infinite at iteration $n=2$ if,
$$v_x^2 \cdot (x_0^2+y_0^2)+v_y^2 \cdot (x_0^2+y_0^2)+1=2 \cdot v_x \cdot |x_0|+2 \cdot v_y \cdot |y_0|$$
If one could set up a suitable renormalization scheme, it might be possible to find all the zeros within some suitable error tolerance.
Show and Tell
Here's my attempt at providing pictures. If you wondering about the quality, it was made on a DS. We have a logarithmic path indicator function. It's identical to yours, except I take a log after we're done calculating it. The range of the graphic is the square from $-1 \le x \le 1$ and $-1 \le y \le 1$. Use the origin to base off what the colors mean, as far the path indicator function is concerned.
Below is what results when the constant vector is $v_x=v_y=0.1$.
Using $v_x=v_y=0.5$, we get qualitatively different behavior.
Let's add some asymmetry with $v_x=0.25, \ v_y=0$. The line at $y=0$ corresponds to the $1$-D case we discussed above.
Fixed and Singularity Points
So the pictures I posted have very intricate and complex patterns. However, there are two pieces of these pictures that can be readily explained. These are the fixed and singularity points. These values are pretty much what they sound like. Fixed points are points where the path indicator function eventually becomes fixed in value. Singularity points are points where the path indicator function becomes infinite after some iteration. Finding these points becomes very difficult so I'll be only giving an intro to this topic. In the one dimensional case, the recurrence $(5)$ has fixed points at,
$$(8) \pm \cfrac{-c+\sqrt{c^2+4}}{2}$$
If your being keen, you'll notice that the negative solution is false. However, if your very keen you'll notice that the negative solution will evolve to the real fixed point in one iteration. Basically a fixed point of $-2$ is false, but after one iteration, it will become $2$ which is a fixed point.
To find the singularity points, we have to find values such that the iteration becomes infinite. In our case, we can do this by finding values such that the iteration is $0$. If we find these values, we know that the very next value will be infinite. Looking at $(5)$ again, we find that for an initial point $v_0$ to evolve to an infinite point after $1$ iteration, $v_0$ must equal $0$. For, the case of two iterations, we solve,
$$0={1 \over {|v_0|}}-c$$
Which implies that $v_0=\pm {1 \over c}$. We can then move on to finding the values of $v_0$ that give infinite values after three iterations and so forth.
Let's use the third picture for an example. We have $v_x=0.25, \ v_y=0$ so we can our above analysis quite nicely. Using $(8)$ we find the fixed points to be
$$\pm {{\sqrt{65}-1} \over 8}=\pm 0.8828...$$
That corresponds to the dark areas on the $y=0$ line in the third picture. Similarly we can find the first three singularity points to be, $0$, $\pm 4$, and $\pm 4/17=0.235...$. By the way, formulas exist to find the next $6$ singularity points.
There are similar, albeit more complicated, formulas for the two and three dimensional cases. For instance, I provide a formula that can be utilized for finding singularity points in another section.
Potpourri
The final interesting thing to note is that the iteration of $v$, in the two dimensional case, can be turned into a complex iteration equivalent to,
$$z_{n+1}=\cfrac{\operatorname{Re}(z_n)+\operatorname{Im}(z_n)}{|z_n|^2}$$
Where we have the real and imaginary part operators involved.
I may or may not add something else...