http://www.bugman123.com/Fractals/Fractals.html
(* runtime: 0.06 second *)
ta = 1.958591030 -
0.011278560I; tb = 2; tab = (ta tb + Sqrt[ta^2tb^2 - 4(ta^2 +
tb^2)])/2; z0 = (tab - 2)tb/(tb tab - 2ta + 2I tab);
b = {{tb - 2I, tb}, {tb, tb + 2I}}/2; B = Inverse[b]; a = {{tab, (tab -
2)/z0}, {(tab + 2)z0, tab}}.B; A = Inverse[a];
Fix[{{a_, b_}, {c_, d_}}] := (a - d - Sqrt[4 b c + (a - d)^2])/(2 c);
ToMatrix[{z_, r_}] := (I/r){{z, r^2 - z Conjugate[z]}, {1,
-Conjugate[z]}};
MotherCircle[M1_, M2_, M3_] := ToMatrix[{x0 +I y0, r}] /.
Solve[Map[(Re[#] - x0)^2 + (Im[#] - y0)^2 == r^2 &, Fix /@ {M1, M2,
M3}], {x0, y0, r}][[2]];
C1 = MotherCircle[b, a.b.A, a.b.A.B]; C2 =
MotherCircle[b.a.a.a.a.a.a.a.a.a.a.a.a.a.a.a,
a.b.a.a.a.a.a.a.a.a.a.a.a.a.a.a, a.b.A.B];
Reflect[C_, M_] := M.C.Inverse[Conjugate[M]];
orbits = Join[Reverse[NestList[Reflect[#, a] &, C1, 63]],
Drop[NestList[Reflect[#, A] &, C1, 63], 1],
Reverse[NestList[Reflect[#, a] &, C2, 71]],
Drop[NestList[Reflect[#, A] &, C2, 56], 1]];
Show[Graphics[MapIndexed[({{a, b}, {c, d}} = #1; {Hue[#2[[1]]/15],
Disk[{Re[a/c], Im[a/c]}, Re[I/c]]}) &, orbits]], PlotRange ->
35{{-1, 1}, {-1, 1}}, AspectRatio -> Automatic];