Fractals using just modulo operation Let us calculate the remainder after division of $27$ by $10$.
$27 \equiv 7 \pmod{10}$
We have $7$. So let's calculate the remainder after divison of $27$ by $7$.
$ 27 \equiv 6 \pmod{7}$
Ok, so let us continue with $6$ as the divisor...
$ 27 \equiv 3 \pmod{6}$
Getting closer...
$ 27 \equiv 0 \pmod{3}$
Good! We have finally reached $0$ which means we are unable to continue the procedure. 
Let's make a function that counts the modulo operations we need to perform until we finally arrive at $0$.
So we find some remainder $r_{1}$ after division of some $a$ by some $b$, then we find some remainder $r_{2}$ after division of $a$ by $r_{1}$ and we repeat the procedure until we find such index $i$ that $r_{i} = 0$.
Therefore, let $$ M(a, b) = i-1$$
for $a, b \in \mathbb{N}, b \neq 0 $
(I like to call it "modulity of a by b", thence M)
For our example: $M(27, 10) = 3$.
Notice that $M(a, b) = 0 \Leftrightarrow  b|a $ (this is why $i-1$ feels nicer to me than just $i$)
Recall what happens if we put a white pixel at such $(x, y)$ that $y|x$:

This is also the plot of $M(x, y) = 0$.
(the image is reflected over x and y axes for aesthetic reasons. $(0, 0)$ is exactly in the center)
What we see here is the common divisor plot that's already been studied extensively by prime number researchers.
Now here's where things start getting interesting:
What if we put a pixel at such $(x, y)$ that $M(x, y) = 1$?

Looks almost like the divisor plot... but get a closer look at the rays. It's like copies of the divisor plot are growing on each of the original line!
How about $M(x, y) = 2$?

Copies are growing on the copies!
Note that I do not overlay any of the images, I just follow this single equation.
Now here is my favorite.
Let us determine luminosity ($0 - 255$) of a pixel at $(x, y)$ by the following equation:
$$255 \over{ M(x,y) + 1 }$$
(it is therefore full white whenever $y$ divides $x$, half-white if $M(x, y) = 1$ and so on)

The full resolution version is around ~35 mb so I couldn't upload it here  (I totally recommend seeing this in 1:1):
https://drive.google.com/file/d/0B_gBQSJQBKcjakVSZG1KUVVoTmM/view?usp=sharing
What strikes me the most is that some black stripes appear in the gray area and they most often represent prime number locations.
Trivia


*

*The above plot with and without prime numbers marked with red stripes:
http://i.imgur.com/E9YIIbd.png
http://i.imgur.com/vDgkT8j.png

*The above plot considering only prime $x$:

Formula: $255 \over{ M(p_{x},y) }$ (note I do not add $1$ to the denominator because it would be full white only at $y$ equal $1$ or the prime. Therefore, the pixel is fully white when $p_{x}$ mod $y = 1$ )
Full 1:1 resolution: https://drive.google.com/file/d/0B_gBQSJQBKcjTWMzc3ZHWmxERjA/view?usp=sharing
Interestingly, these modulities form a divisor plot of their own. 

*Notice that for $ M(a, b) = i-1, r_{i-1}$ results in either $1$ or a divisor of $a$ (which is neither $1$ nor $a$).
I put a white pixel at such $(x, y)$ that for $M(x, y) = i - 1$, it is true that $r_{i-1}\neq 1 \wedge r_{i-1} | x$ (the one before last iteration results in a remainder that divides $x$ and is not $1$ (the uninteresting case))
http://i.imgur.com/I85rlH5.png
It is worth our notice that growth of $M(a, b)$ is rather slow and so if we could discover a rule by which to describe a suitable $b$ that most often leads to encountering a proper factor of $a$, we would discover a primality test that works really fast (it'd be $O(M(a, b))$ because we'd just need to calculate this $r_{i-1}$).

*Think of $M'(a, b)$ as a function that does not calculate $a$ mod $b$ but instead does $M(a, b)$ until a zero is found.
These two are plots of $M'''(x, y)$ with and without primes marked:
http://i.imgur.com/gE0Bvwg.png
http://i.imgur.com/vb5YxVP.png

*Plot of $M(x, 11)$, enlarged 5 times vertically:
http://i.imgur.com/K2ghJqe.png
Can't notice any periodicity in the first 1920 values even though it's just 11.
For comparison, plot of $x$ mod $11$ (1:1 scale):
http://i.imgur.com/KM6SCF3.png

*As it's been pointed out in the comments, subsequent iterations of $M(a, b)$ look very much like Euclidean algorithm for finding the greatest common divisors using repeated modulo. A strikingly similar result can be obtained if for $(x, y)$ we plot the number of steps of $gcd(x, y)$: 

I've also found similar picture on wikipedia:

This is basically the plot of algorithmic efficiency of $gcd$.
Somebody even drew a density plot here on stackexchange.
The primes, however, are not so clearly visible in GCD plots. Overall, they seem more orderly and stripes do not align vertically like they do when we use $M(a, b)$ instead.
Here's a convenient comparative animation between GCD timer (complexity plot) and my Modulity function ($M(x, y)$). Best viewed in 1:1 zoom. $M(x, y)$ appears to be different in nature from Euclid's GCD algorithm.

Questions


*

*Where is $M(a, b)$ used in mathematics? 

*Is it already named somehow?

*How could one estimate growth of $M(a, b)$ with relation to both $a$ and $b$, or with just $a$ increasing?

*What interesting properties could $M(a, b)$ possibly have and could it be of any significance to number theory?

 A: Some additional notes, which I cannot add to my previous answer, because apparently I am close to a 30K character limit and MathJax complains.

Addendum
The fundamental pattern which emerges in $\phi(n)$ then, is that of the Farey series dissection of the continuum. This pattern is naturally related to Euclid's Orchard.
Euclid's Orchard is basically a listing of the Farey sequence of (all irreducible) fractions $p_k/q_k$ for the unit interval, with height equal to $1/q_k$, at level $k$:

Euclid's Orchard on [0,1].
In turn, Euclid's Orchard is related to Thomae's Function and to Dirichlet's Function:

Thomae's Function on [0,1].
The emergence of this pattern can be seen easier in a combined plot, that of the GCD timer and Euclid's Orchard on the unit interval:

Farey series dissection of the continuum of [0,1].
Euclid's Orchard is a fractal. It is the most "elementary" fractal in a sense, since it characterizes the structure of the unit interval, which is essential in understanding the continuum of the real line.
Follow some animated gifs which show zoom-in at specific numbers:

Zoom-in at $\sqrt{2}-1$.

Zoom-in at $\phi-1$.
The point of convergence of the zoom is marked by a red indicator.
White vertical lines which show up during the zoom-in, are (a finite number of open covers of) irrationals. They are the irrational limits of the convergents of the corresponding continued fractions which are formed by considering any particular tower-top path that descends down to the irrational limit.
In other words, a zoom-in such as those shown, displays some specific continued fraction decompositions for the (irrational on these two) target of the zoom.
The corresponding continued fraction decomposition (and its associated convergents) are given by the tops of the highest towers, which descend to the limits.

Addendum #2 (for your last comment to my previous answer)
For the difference between the two kinds of graphs you are getting - because I am fairly certain you are still confused - what you are doing produces two different kinds of graphs. If you set M(x,y) to their natural value, you are forcing a smooth graph like the the GCD timer. If you start modifying M(x,y) or set it to other values (($M(x,y)=k$ or if you calculate as $M(x,p^k)$), you will begin reproducing vertical artifacts which are characteristic of $\phi$. And that, because as you correctly observe, doing so, you start dissecting the horizontal continuum as well (in the above case according to $p^k$, etc). In this case, the appropriate figure which reveals the vertical cuts, would be like the following:


Appendix:
Some Maple procedures for numerical verification of some of the theorems and for the generation of some of the figures.
generate fig.1:

with(numtheory): with(plots): N:=10000;
  liste:=[seq([n,phi(n)],n=1..N)]: with(plots):#Generate fig.1
  p:=plot(liste,style=point, symbol=POINT,color=BLACK): display(p);

Generate fig.2:

q:=plot({x-1,x/2,x/3,2*x/3,2*x/5, 4*x/5,4*x/15,8*x/15,2*x/7,3*x/7,
  4*x/7,6*x/7,8*x/35,12*x/35,16*x/35,24*x/35},x=0..N,color=grey):
  display({p,q});#p as in example 1.

Generate fig.3:

F:=proc(n) #Farey series local a,a1,b,b1,c,c1,d,d1,k,L;
  a:=0;b:=1;c:=1;d:=n;L:=[0]; while (c < n) do  k:=floor((n+b)/d); 
  a1:=c;b1:=d;c1:=kc-a;d1:=kd-b;a:=a1;b:=b1;c:=c1;d:=d1; 
  L:=[op(L),a/b];  od:  L; end:
  n:=10;
  for m from 1 to nops(F(n)) do  f:=(m,x)->F(n)[m]*x; od:
  q:={}; with(plots): for m from 1 to nops(F(n)) do 
  qn:=plot(f(m,x),x=0..10000,color=grey);  q:=q union {qn}; od: 
  display(p,q);

Implements Theorem 4.1:

S:=proc(L,N) local LS,k,ub; LS:=nops(L);#find how many arguments if
  LS=1 then  floor(logL[LS]); else  ub:=floor(logL[LS]); 
  add(S(L[1..LS-1],floor(N/L[LS]^k)),k=1..ub); fi; end:

Brute force approach for Theorem 4.1:

search3:=proc(N,a1,a2,a3,s) local cp,k1,k2,k3; cp:=0; for k1 from 1 to
  s do  for k2 from 1 to s do   for k3 from 1 to s do    if
  a1^k1*a2^k2*a3^k3 <= N then
  cp:=cp+1;fi;od;  od; od; cp; end:

Verify Theorem 4.1:

L:=[5,6,10];N:=1738412;S(L,N); 37 s:=50 #maximum exponent for brute
  force search search3(N,5,6,10,s); 37 #identical

Times GCD algorithm:

reduce:=proc(m,n) local T,M,N,c; M:=m/gcd(m,n);#GCD(km,kn)=GCD(m,n)
  N:=n/gcd(m,n); c:=0; while M>1 do  T:=M;  M:=N;  N:=T;#flip  M:=M mod
  N;#reduce  c:=c+1; od; c; end:

Generate fig.6:

max:=200; nmax:=200; rt:=array(1..mmax,1..nmax); for m from 1 to mmax
  do  for n from 1 to nmax do   rt[m,n]:=reduce(n,m); # assign GCD steps
  to array  od; od;
  n:='n';m:='m'; rz:=(m,n)->rt[m,n]; # convert GCD steps to function
  p:=plot3d(rz(m,n),
      m=1..mmax,n=1..nmax,
      grid=[mmax,nmax],
      axes=NORMAL,
      orientation=[-90,0],
      shading= ZGREYSCALE,
      style=PATCHNOGRID,
      scaling=CONSTRAINED): display(p);


