Asmus Schmidt's Complex Continued Fraction Algorithm

Implementing the Schmidt algorithm for approximation of complex numbers by elements of Q(i)

by Doug Hensley

Explanation

This implementation is restricted to the case of Gaussian rational numbers w inside the square with corners at 0,1,1+i, and i. It returns  a list containing the approximation target w, the first Schmidt matrix, global and zoomed-in graphics structures, and  the sequence of T's needed to build the successive M's.

Asmus Schmidt's Complex Continued Fraction algorithm takes as input a complex number which we shall call w, and returns a sequence of matrices ((Subscript[T, n])) that encode a sequence of essentially best-possible Gaussian rational approximations ((Subscript[r, n])) to w, each of the form (a+i b)/(c+i d). 

In this implementation, we restrict attention to w drawn from a basic square, with opposite corners at 0 and 1+i. The T's are drawn from a list of 7 choices detailed below.  Each of these 7 has Gaussian integer entries and determinant 1,-1,i, or -i, and thus has an inverse with Gaussian integer entries. With these T's, one constructs a further sequence ((Subscript[M, n])), where (Subscript[M, n]) = (Subscript[T, 1]) (Subscript[T, 2]) ...  (Subscript[T, n]). From each (Subscript[M, n])=(({{(p^((1))), (p^((2)))}, {(q^((1))), (q^((2)))}})), one constructs three Gaussian rationals (r^((1)))=(p^((1)))/(q^((1))), (r^((2))) = (p^((2)))/(q^((2))), and (r^((3))) = ((p^((1))) + 
                (p^((2))))/((q^((1))) + (q^((2)))) = 
          (p^((3)))/(q^((3))), say.  Equivalently, {{p1,p2,p3},{q1,q2,q3}}=(Subscript[M, n]).{{1,0,1},{0,1,1}}. Whichever of these three candidate Gaussian rational approximations is nearest w is our (Subscript[r, n]). 

The T's are chosen so that certain regions (Subscript[F, n]) shall contain w. If det((Subscript[M, n]))=1 or -1, the region is the {Schmidt disk} of (Subscript[M, n]), the solid disk bounded by the circle through (r^((1))), (r^((2))), and (r^((3))). If instead, det((Subscript[M, n]))=i or -i, the region is a {Schmidt triangle} of (Subscript[M, n]) , a triangle-like region bounded by arcs, and with vertices (r^((1))), (r^((2))), and (r^((3))). 

The Schmidt triangle of (Subscript[M, n]) is the closed, bounded region contained within three arcs, from  (r^((1))) to  (r^((2))), (r^((2)))to (r^((3))), and (r^((3)))to  (r^((1))). These arcs can be either arcs of a circle, or, exceptionally, segments of a straight line. Each arc must be perpendicular at both ends to the circle through  (r^((1))), (r^((2))), and (r^((3))), and each lies within the disk rather than through its complement. 

Each Schmidt disk has six distinguished points: the first three are  (r^((1))), (r^((2))), and (r^((3))), and the next three are termed (r^((1)))', (r^((2)))', and (r^((3)))'. These last three are given by {{p1',p2',p3'},{q1',q2',q3'}}=(Subscript[M, n]).c.{{1,0,1},{0,1,1}}, (c={{1,-1+i},{1-i,1}}), and are inside the Schmidt disk and on the edges of the Schmidt triangle, depending on whether M has real or imaginary determinant. 

Each Schmidt disk D is partitioned into 7 regions, consisting of 3 smaller Schmidt disks D1, D2, D3 tangent to D at r1, r2, r3 respectively, and tangent to each other as well, and 4 Schmidt triangles, one occupying the central region bounded by arcs of D1, D2, and D3.

D1 and D2 are tangent at r3', D2 and D3 at r1', and D3 and D1 tangent at r2'.  The four Schmidt triangles G1,G2,G3, and GC, say, have vertices r3r1'r2, r1r2'r3, r2r3'r1, and r1'r2'r3' respectively. 


The figure immediately below shows the situation after the third step in the execution of the Schmidt algorithm when w=E^i =0.540302+0.841471 i. We have 
{T1,T2,T3}={c,c,v3}, so that \(TraditionalForm\`M\_3 = \(\({{3 + 
                3\ i, \(-4\) - i}, {4 - i, \(-3\) + 3\ i}}\)\(.\)\(\ \)\)\)Thus (Subscript[F, 3])is the disk shown. Now w sits inside the central triangle of the picture, so that T4=c, and (Subscript[M, 4])={{-2+6 i,-5-4 i},{4+5 i,-6+2 i}}. The figures next below show the situation at the sixth generation. Now, the T-sequence is {c,c,v3,c,c,v3}. The left picture shows the overall layout of the nested zones bracketing w, the middle  picture zooms in to show the top picture and further iterations, and the right picture shows a tighter zoom on the sixth generation situation. The figures at bottom show the situation as of the 20th generation; the window sizes vary from 10^(-4)*2^(-4), to 10^(-4)*2^(-9) by powers of 2. 

It should be noted that E^i, our approximation target in this series of examples, has a highly atypical Schmidt expansion. Schmidt was able to show that the T-sequence for this number is \(TraditionalForm\`\(C\^2\) \(VC\^2\) \(V\^3\) \(C\^2\) \(V\^5\) \(C\^2\) V\^7 ... \)

The algorithm makes somewhat slower progress when it encounters a string of (Subscript[E, j])'s, with fixed j. This can happen if w finds itself nestled nearly at the vertex of a Schmidt triangle at iteration n; several iterations of a particular (Subscript[E, j]) will be needed. Happily, the area of such regions is small so the average performance of the algorithm is not much affected.

Notation tracks closely that used by Schmidt where possible.The matrices v1,v2,v3,e1,e2,e3,and c are the same as in[Diophantine Approximation of Complex Numbers].Other notation tries to be helpful.frax[a,b] gives ∞ if b=0,else a/b.ctrnradsqrd gives a list,{center,square of the radius} for the center and radius-squared of the circle through three complex numbers.Degenerate cases return {∞,∞}.The function "inside" takes as input a complex number w and a circle defined by three complex numbers on its perimeter,and returns either 1,0,or-1 depending on whether w lies inside,on the rim,or outside the circle.It should not be given a degenerate circle.

Off[General::spell]

(Many of the names assigned in this notebook are variations on a theme. )

Initialization

v1 = {{1, I}, {0, 1}}; v2 = {{1, 0}, {-I, 1}}; 
    v3 = {{1 - I, I}, {-I, 1 + I}}

e1 = {{1, 0}, {1 - I, I}}; e2 = {{1, -1 + I}, {0, I}}; 
    e3 = {{I, 0}, {0, 1}}; c = {{1, -1 + I}, {1 - I, I}}

schmidtmatrix[1] = v1; schmidtmatrix[2] = v2; schmidtmatrix[3] = v3; 
    schmidtmatrix[4] = e1; schmidtmatrix[5] = e2; schmidtmatrix[6] = e3; 
    schmidtmatrix[7] = c;
schmidtmatrix[8] = v1;
schmidtmatrix[9] = v2;
schmidtmatrix[10] = v3;
schmidtmatrix[11] = c;

The T's of the Schmidt algorithm are drawn from this list of choices, v1-v3,e1-e3, and c.

frax[a_, b_] := If[b == 0, ComplexInfinity, a/b]

Now we give the code which got the pictures. The function "rvals" converts a matrix M into the list {{r1,r2,r3},{r1',r2',r3'}}.

rvals[m_] := 
      Module[{p1, q1, p2, q2, p3, q3, p1prime, q1prime, p2prime, q2prime, 
          p3prime, q3prime, r1, r2, r3, r1prime, r2prime, 
          r3prime}, 
{{p1, p2, p3}, {q1, q2, q3}} = 
          m . {{1, 0, 1}, {0, 1, 1}}; 
{{p1prime, p2prime, 
              p3prime}, {q1prime, q2prime, q3prime}} = 
          m . c . {{1, 0, 1}, {0, 1, 1}}; 
r1 = 
          frax[p1, q1]; r2 = frax[p2, q2]; 
        r3 = frax[p3, q3]; 
r1prime = 
          frax[p1prime, q1prime]; r2prime = frax[p2prime, q2prime]; 
        r3prime = 
          frax[p3prime, q3prime]; 
{{r1, r2, r3}, {r1prime, 
            r2prime, r3prime}}]

The functions cleanup and cleanupb remove text and point elements from graphics, for cleaner display.

cleanup[pix_] := 
      Select[pix, 
        Head[#] == Circle || Head[#] == Line || 
            Head[#] == Point || 
            Head[#] == PointSize &]

cleanupb[pix_] := 
      Select[pix, 
        Head[#1] == Circle || Head[#1] == Line &]

nearest[z_] := Floor[z + 1/2]

The function "linedup" determines whether or not three points are collinear.

linedup[zs_] := 
      Module[{z1, z2, z3, s}, 
{z1, z2, z3} = 
          zs; 
s = (z3 - z1)/(z2 - 
                z1); 
(Im[s] == 0 && 
            Re[s] ≥ 1)]

reim[z_] := {Re[z], Im[z]};

midpt[u_, v_] := (u + v)/2;
diff[v1_, v2_] := v2 - v1

\(centernradsqrd[{z1_, z2_, z3_}] := 
      Module[{x1, y1, x2, y2, x3, y3, s, t, d, u, w, v1, v2, ctr, rr, 
          ck}, {{x1, y1}, {x2, y2}, {x3, y3}} = 
          Map[reim, {z1, z2, 
              z3}]; 
 (*\(Print["centernradsqrd called on inputs", {z1, z2, z3}];\)*) 
v1 = {x2, y2} - {x1, y1}; 
        v2 = {x3, y3} - {x2, y2}; 
d = 
          Det[{v1, v2}]; 
If[d == 0, 
          Return[{ComplexInfinity, 
              ComplexInfinity (*\(,\)\(True\)*) }]]; 
u = {x3, y3} - {x1, y1}; 
s = v2 . u/2; s = s/d; t = v1 . u/2; 
        t = t/d; 
ctr = \((z1 + z2)\)/2 + 
            s\ I \((z2 - z1)\); 
rr = 
          FullSimplify[
            Abs[z1 - 
                  ctr]^2]; 
 (*\(ck = \((FullSimplify[
                      rr - Abs[z2 - ctr]^2] == 0 && 
                  FullSimplify[rr - Abs[z3 - ctr]^2] == 
                    0)\);\)*) 
{ctr, 
          rr (*\(,\)\(ck\)*) }]\)

The function "inside" returns 1 if the point is in the interior of the circle through the given points, 0 if it lies on the boundary, and -1 if it lies outside.

cutup[m_] := 
      Module[{t}, 
        t = Det[m]; 
If[
          Re[t] ≠ 0, {m . c, m . v1, m . v2, m . v3, m . e1, 
            m . e2, m . e3}, {m . c, m . v1, m . v2, m . v3}]]

showcircle[m_] := 
      Module[{c0, r1, r2, r3}, 
        If[Im[Det[m]] ≠ 0, 
          Return[Point[{0, 0}]]]; {r1, r2, r3} = rvals[m][[1]]; 
        If[r1 == ComplexInfinity || r2 == ComplexInfinity || 
            r3 == ComplexInfinity, 
          Print["complex infinity, lions and tigers and bears oh my"]; 
          Return[Point[{0, 0}]]]; 
If[
          collinear[{r1, r2, r3}], Return[Line[{r1, r2, r3}]]]; 
        c0 = centernradsqrd[{r1, r2, r3}]; 
c0 = {reim[
              c0[[1]]], Sqrt[c0[[2]]]}; 
Circle[
          c0[[1]], c0[[2]]]]

collinear[{r1_, r2_, r3_}] := 
      Module[{u, v}, 
        If[r1 == ComplexInfinity || r2 == ComplexInfinity || 
            r3 == ComplexInfinity, Print["bears"]; Return[True]]; 
        u = reim[r2 - r1]; v = reim[r3 - r1]; 
        Det[{u, v}] == 0]

bigenough[m_, eps_] := 
      Module[{c0, r1, r2, r3, mc, u, v, w}, 
        mc = If[Im[Det[m]] ≠ 0, m . c, m]; {r1, r2, 
            r3} = rvals[mc][[1]]; 
If[
          collinear[{r1, r2, r3}], 
          Print["m=", m, "has collinear r values but circle type"]; 
          Return[False]]; 
c0 = 
          centernradsqrd[{r1, r2, 
              r3}]; 
c0[[2]] ≥ 
          eps^2]

nextgeneration[mlist_, eps_] := 
      Module[{newms, goodms}, newms = Flatten[Map[cutup, mlist], 1]; 
        goodms = Select[newms, bigenough[#, eps] &]; goodms]

bigpicture[g0_, eps_] := 
      Module[{gg, hh}, hh = g0; 
        gg = g0; 
While[gg ≠ {}, 
          gg = nextgeneration[gg, eps]; hh = Union[gg, hh]]; hh]

Visualization - Schmidt Circle

g0 = {v3}

g1 = nextgeneration[g0, 0]

g2 = nextgeneration[g1, 0];
g3 = nextgeneration[g2, 0];
g4 = nextgeneration[g3, 0];

Length /@  {g0,  g1,  g2,  g3,  g4}

To keep successive generations from exploding, we make the smallest circle radius 1/60

cutoff  =  1/300;
g2 = nextgeneration[g1, cutoff];
g3 = nextgeneration[g2, cutoff];
g4 = nextgeneration[g3, cutoff];
g5 = nextgeneration[g4, cutoff];
g6 = nextgeneration[g5, cutoff];
g7 = nextgeneration[g6, cutoff];
g8 = nextgeneration[g7, cutoff];
g9 = nextgeneration[g8, cutoff];
g10 = nextgeneration[g9, cutoff];
g11 = nextgeneration[g10, cutoff];
g12 = nextgeneration[g11, cutoff];
g13 = nextgeneration[g12, cutoff];
g14 = nextgeneration[g13, cutoff];
g15 = nextgeneration[g14, cutoff];
g16 = nextgeneration[g15, cutoff];
g17 = nextgeneration[g16, cutoff];
g18 = nextgeneration[g17, cutoff];
g19 = nextgeneration[g18, cutoff];
g20 = nextgeneration[g19, cutoff];
g21 = nextgeneration[g20, cutoff];
g22 = nextgeneration[g21, cutoff];
g23 = nextgeneration[g22, cutoff];
g24 = nextgeneration[g23, cutoff];
g25 = nextgeneration[g24, cutoff];
g26 = nextgeneration[g25, cutoff];
g27 = nextgeneration[g26, cutoff];

allgenerations  =  {g0,  g1,  g2,  g3,  g4,  g5,  g6,  g7,  
          g8,  g9,  g10,  g11,  g12,  g13,  g14,  g15,  g16,  g17,  
          g18,  g19,  g20, g21,  g22,  g23,  g24,  g25,  g26,  
          g27};

Length /@  allgenerations

pixygen = 
        Drop[Union[
            Map[showcircle, 
              Union[Flatten[allgenerations, 1]]]], -1];

Show[Graphics[pixygen], AspectRatio -> Automatic]

Code - deeper triangles

inside[w_, zs_] := 
      Module[{ctr, rr, s}, 
        If[linedup[zs], 
s = 
            Im[(w - zs[[1]])/(zs[[2]] - 
                    zs[[1]])]; 
Return[
            Sign[s]]]; 
{ctr, rr  } = 
          centernradsqrd[zs]; 
Sign[
          FullSimplify[rr - Abs[w - ctr]^2]]]

schmidttriangle[w_, m_] := 
      Module[{r, rp, ck0, ck1, ck2, ck3}, {r, rp} = 
          rvals[m]; 
ck0 = 
          inside[w, r]; 
ck1 = 
          inside[w, {r[[1]], rp[[2]], 
              r[[3]]}]; 
ck2 = 
          inside[w, {r[[2]], rp[[3]], 
              r[[1]]}]; 
ck3 = 
          inside[w, {r[[3]], rp[[2]] , 
              r[[1]]}]; 
(ck0 ≥ 0 && 
            ck1 ≤ 0 && ck2 ≤ 0 && 
            ck3 ≤ 0)]

This isn't used. But it's the core of the next bit of code. The second version of schmidttriangle says whether w lies inside a specific triangle zone associated with the matrix m.

\(schmidttriangle[w_, m_, k_] := 
      Module[{r, rp, ck0, ck1, ck2, ck3, c1, c2, c3, sk1, sk2, sk3, z1, z2, 
          z3, z4, z5, z6, sk4, sk5, sk6, ck4, ck5, ck6, 
          type}, 
type = \((Re[Det[m]] ≠ 
              0)\); 
If[type && k > 7, 
          Print["circle type m, use k from 1 to 7"]; 
          Return[False]]; 
If[
          Not[type] && k ≤ 7, 
          Print["triangle type m, use k from 8 to 11"]; 
Return[False]]; 
{r, rp} = 
          rvals[m]; 
ck0 = inside[w, r]; 
        If[ck0 < 0, Return[False]]; 
ck1 = 
          inside[w, {r[\([1]\)], rp[\([3]\)], 
              rp[\([2]\)]}]; 
If[k == 1, 
          Return[\((ck1 ≥ 0)\)]]; 
ck2 = 
          inside[w, {r[\([2]\)], rp[\([1]\)], 
              rp[\([3]\)]}]; 
If[k == 2, 
          Return[\((ck2 ≥ 0)\)]]; 
ck3 = 
          inside[w, {r[\([3]\)], rp[\([2]\)], 
              rp[\([1]\)]}]; 
If[k == 3, 
          Return[\((ck3 ≥ 
                0)\)]]; 
z1 = \(centernradsqrd[{r[\([1]\)], 
                rp[\([3]\)], 
                rp[\([2]\)]}]\)[\([1]\)]; 
z2 = \(centernradsqrd[{r[\([2]\)], rp[\([1]\)], 
                rp[\([3]\)]}]\)[\([1]\)]; 
z3 = \(centernradsqrd[{r[\([3]\)], rp[\([2]\)], 
                rp[\([1]\)]}]\)[\([1]\)]; 
sk1 = 
          Sign[Im[\((w - z3)\)/\((z2 - z3)\)]]; 
sk2 = 
          Sign[Im[\((w - z1)\)/\((z3 - z1)\)]]; 
sk3 = 
          Sign[Im[\((w - z2)\)/\((z1 - z2)\)]]; 
If[
          k == 4, 
          Return[\((ck2 ≤ 0 && ck3 ≤ 0 && 
                sk1 ≥ 0)\)]]; 
If[
          k == 5, 
          Return[\((ck1 ≤ 0 && ck3 ≤ 0 && 
                sk2 ≥ 0)\)]]; 
If[
          k == 6, 
          Return[\((ck1 ≤ 0 && ck2 ≤ 0 && 
                sk3 ≥ 0)\)]]; 
If[
          k == 7, Return[True]]; 
ck4 = 
          inside[w, {r[\([3]\)], rp[\([1]\)], 
              r[\([2]\)]}]; 
If[ck4 > 0, 
          Return[False]]; 
ck5 = 
          inside[w, {r[\([1]\)], rp[\([2]\)], 
              r[\([3]\)]}]; 
If[ck5 > 0, 
          Return[False]]; 
ck6 = 
          inside[w, {r[\([2]\)], rp[\([3]\)], 
              r[\([1]\)]}]; 
If[ck6 > 0, 
          Return[False]]; 
If[k == 11, 
          Return[inside[w, rp] ≥ 0]]; 
        If[k < 7 && inside[w, rp] > 0, Return[False]]; 
If[
          k == 8, 
          Return[Im[\((w - rp[\([3]\)])\)/\((rp[\([2]\)] - 
                      rp[\([3]\)])\)] ≥ 
              0]]; 
If[k == 9, 
          Return[Im[\((w - rp[\([1]\)])\)/\((rp[\([3]\)] - 
                      rp[\([1]\)])\)] ≥ 
              0]]; 
If[k == 10, 
          Return[Im[\((w - rp[\([2]\)])\)/\((rp[\([1]\)] - 
                      rp[\([2]\)])\)] ≥ 
              0]]; 
Print["search for zone failed, algorithm defective"]; Return[False]]\)

circlepix[w_, m_] := 
      Module[{rs, u, v, r1, r2, r3, r1p, r2p, r3p, c0, c1, c2, c3, d1, d2, 
          d3, det}, 
det = Det[m]; 
        If[Re[det] == 
            0,   
            Return[Null]]; 
rs = 
          rvals[m]; 
  
{{r1, r2, r3}, {r1p, r2p, r3p}} = 
          rs; 
If[linedup[{r1, r2, r3}], 
          c0 = Line[{reim[r1], reim[r2], reim[r3]}]]; 
c0 = 
          centernradsqrd[{r1, r2, r3}]; 
c0 = {reim[
              c0[[1]]], Sqrt[c0[[2]]]}; 
c0 = 
          Circle[c0[[1]], c0[[2]]]; 
If[
          linedup[{r1, r3p, r2p}], 
          c1 = Line[{reim[r1], reim[r3p], 
                reim[r2p]}]]; 
c1 = 
          centernradsqrd[{r1, r3p, r2p}]; 
c1 = {reim[
              c1[[1]]], Sqrt[c1[[2]]]}; 
        d1 = c1; 
c1 = 
          Circle[c1[[1]], c1[[2]]]; 
If[
          linedup[{r2, r1p, r3p}], 
          c2 = Line[{reim[r2], reim[r1p], 
                reim[r3p]}]]; 
c2 = 
          centernradsqrd[{r2, r1p, r3p}]; 
c2 = {reim[
              c2[[1]]], Sqrt[c2[[2]]]}; 
        d2 = c2; 
c2 = 
          Circle[c2[[1]], c2[[2]]]; 
If[
          linedup[{r3, r2p, r1p}], 
          c3 = Line[{reim[r3], reim[r2p], 
                reim[r1p]}]]; 
c3 = 
          centernradsqrd[{r3, r2p, r1p}]; 
c3 = {reim[
              c3[[1]]], Sqrt[c3[[2]]]}; 
        d3 = c3; 
c3 = 
          Circle[c3[[1]], c3[[2]]]; 
  ]

circlepix2[w_, m_] := 
      Module[{rs, u, v, r1, r2, r3, r1p, r2p, r3p, c0, c1, c2, c3, d1, d2, 
          d3, det}, 
det = Det[m]; 
        If[Re[det] == 
            0,   
            Return[Null]]; 
rs = 
          rvals[m]; 
  
{{r1, r2, r3}, {r1p, r2p, r3p}} = 
          rs; 
If[linedup[{r1, r2, r3}], 
          c0 = Line[{reim[r1], reim[r2], reim[r3]}]]; 
c0 = 
          centernradsqrd[{r1, r2, r3}]; 
c0 = {reim[
              c0[[1]]], Sqrt[c0[[2]]]}; 
c0 = 
          Circle[c0[[1]], c0[[2]]]; 
If[
          linedup[{r1, r3p, r2p}], 
          c1 = Line[{reim[r1], reim[r3p], 
                reim[r2p]}]]; 
c1 = 
          centernradsqrd[{r1, r3p, r2p}]; 
c1 = {reim[
              c1[[1]]], Sqrt[c1[[2]]]}; 
        d1 = c1; 
c1 = 
          Circle[c1[[1]], c1[[2]]]; 
If[
          linedup[{r2, r1p, r3p}], 
          c2 = Line[{reim[r2], reim[r1p], 
                reim[r3p]}]]; 
c2 = 
          centernradsqrd[{r2, r1p, r3p}]; 
c2 = {reim[
              c2[[1]]], Sqrt[c2[[2]]]}; 
        d2 = c2; 
c2 = 
          Circle[c2[[1]], c2[[2]]]; 
If[
          linedup[{r3, r2p, r1p}], 
          c3 = Line[{reim[r3], reim[r2p], 
                reim[r1p]}]]; 
c3 = 
          centernradsqrd[{r3, r2p, r1p}]; 
c3 = {reim[
              c3[[1]]], Sqrt[c3[[2]]]}; 
        d3 = c3; 
c3 = 
          Circle[c3[[1]], 
            c3[[2]]]; 
{PointSize[ .05 Sqrt[
                c0[[2]]]], Point[reim[r1]], Point[reim[r2]], 
          Point[reim[r3]], 
          Text["w", reim[w]], 
Text["1", 
            d1[[1]]], 
          Text["2", d2[[1]]], 
Text["3", 
            d3[[1]]], c0, c1, c2, c3}]

\(\(\(ark[zs_] := 
      Module[{z1, z2, z3, theta1, theta2, ctr, rad, 
          cr}, 
If[
          linedup[zs],  (*\(Print["ark fixing to return a line thru{", 
                reim[zs[\([1]\)]], ",", 
                reim[zs[\([2]\)]], "}"];\)*) 
Return[
Line[
{reim[zs[\([1]\)]], 
                reim[zs[\([2]\)]], 
                reim[zs[\([3]\)]]}
	]
]]; 
cr = 
          centernradsqrd[zs]; 
ctr = 
          reim[cr[\([1]\)]]; 
rad = 
          Sqrt[cr[\([2]\)]]; 
theta1 = 
          Arg[zs[\([1]\)] - cr[\([1]\)]]; 
theta2 = 
          Arg[zs[\([3]\)] - cr[\([1]\)]]; 
If[
          theta2 < theta1, 
          theta2 += 
            2\ Pi]; 
 (*Print["{theta1,theta2}=", {theta1, theta2}]; 
Print["ark fixing to return a circle with parameters {ctr,rad,{theta1,theta2}}=", 
{ctr, 
              rad, {theta1, theta2}}];*) 
Circle[ctr, 
          rad, {theta1, theta2}]]\)\(
\)
    \)\)

arc[zs_] := 
      Module[{u, v, x, y, z}, 
If[
          linedup[zs],   
Return[
            Line[Map[reim, zs]]]]; 
{x, y, z} = 
          zs; 
u = y - x; 
        v = z - x; 
v = v/u; 
If[
          Im[v] > 0,   Return[
            ark[zs]]]; ark[Reverse[zs]]]

insideschmidttriangle[w_, {z1_, z2_, z3_}] := 
      Module[{u, s, t}, 
If[
          w == z1 || w == z2 || w == z3, 
          Return[True]]; 
s = (z2 - z1)/(w - 
                z1); 
u = (z2 - z1)/(z3 - 
                z1); 
t = (s - 1)/(u - 
                1); 
t = 
          reim[t]; 
(t[[1]] - 1/2)^2 + 
              t[[2]]^2 ≥ 1/4 && 
          t[[2]] ≥ 0 && t[[1]] ≥ 0 && 
          t[[1]] ≤ 1]

trianglepix3[w_, zs_] := 
      Module[{z1, z2, z3, u, v, a1, a2, 
          a3}, 
{z1, z2, z3} = zs; 
If[
          linedup[zs], Print["no triangle, the inputs are collinear"]; 
          Return[False]]; 
u = (z3 - z1)/(z2 - 
                z1); 
v = 
          1/u; 
a1 = {1, 1/(1 - I (1 - v)), 
            0}; 
a2 = {u, 1/(v - I (1 - v)), 
            0}; 
a3 = {u, 
            1/((1 + v)/2 - I (1 - v)/2), 
            1}; 
{a1, a2, 
            a3} = (z2 - z1) {a1, a2, a3}; 
{a1, a2, 
            a3} = z1 + {a1, a2, a3}; 
{Text["w", 
            reim[w]], arc[a1], arc[a2], arc[a3]}]

trianglepix4[w_, m_] := 
      Module[{rs, circ, triangle}, 
        rs = rvals[m]; 
triangle = 
          trianglepix3[w, rs[[1]]]; 
circ = 
          centernradsqrd[
            rs[[2]]]; 
circ = {reim[circ[[1]]], 
            Sqrt[circ[[2]]]}; 
circ = 
          Circle[circ[[1]], 
            circ[[2]]]; 
  
  {Text["w", reim[w]], triangle, 
          circ}  ]

trianglepix5[w_, m_] := 
      Module[{rs, circ, triangle}, 
        rs = rvals[m]; 
triangle = 
          trianglepix3[w, rs[[1]]]; 
circ = 
          centernradsqrd[
            rs[[2]]]; 
circ = {reim[circ[[1]]], 
            Sqrt[circ[[2]]]}; 
circ = 
          Circle[circ[[1]], 
            circ[[2]]]; 
  
{Text["w", reim[w]], 
          triangle, circ}]

zset[m_, k_] := 
      Module[{rs, zs, r1, r2, r3, r1p, r2p, r3p, 
          type}, 
type = (Re[Det[m]] ≠ 0); 
        rs = rvals[m]; 
{{r1, r2, r3}, {r1p, r2p, r3p}} = 
          rs; 
If[
          type, 
If[k > 7, 
            Print["circle type matrix, k should be between 1 and 7"]; 
            Return[False]]; 
If[k == 1, 
            Return[{r1, r3p, r2p}]]; 
If[k == 2, 
            Return[{r2, r1p, r3p}]]; 
If[k == 3, 
            Return[{r3, r2p, r1p}]]; 
If[k == 4, 
            Return[{r3, r1p, r2}]]; 
If[k == 5, 
            Return[{r1, r2p, r3}]]; 
If[k == 6, 
            Return[{r2, r3p, r1}]]; 
If[k == 7, 
            Return[{r1p, r2p, 
                r3p}]]
]; 
If[
          Not[type], 
If[k ≤ 7 || k > 11, 
            Print["triangle type matrix, k should be between 8 and 11"]; 
            Return[False]]; 
If[k == 8, 
            Return[{r1, r3p, r2p}]]; 
If[k == 9, 
            Return[{r2, r1p, r3p}]]; 
If[k == 10, 
            Return[{r3, r2p, r1p}]]; 
If[k == 11, 
            Return[{r1p, r2p, r3p}]]]]

ei = N[E^I, 100]

ei = 10^(-100) nearest[10^100  ei]

N[ei]

c

m = c; w = ei;

inside[2 + I/2, {1, 1 + I, (1 + I)/2}]

rangy[ww_, eps_] := 
      Module[{x, y}, {x, y} = 
          reim[ww]; {{x - eps, x + eps}, {y - eps, y + eps}}]

The function schmidtstart should be given a Gaussian rational within the starting box. It then returns a list with 5 elements, suitable for continuation of the algorithm by iterating the next function, schmidtstep.

schmidtstart[z_] := 
        Module[{w, x, y, m, slist}, 
If[
            Re[z] > 1 || Re[z] < 0 || Im[z] > 1 || 
              Im[z] < 0, 
Print["z lies outside the \
starting square"]; Return[False]]; 
If[
            inside[z, {0, (1 + I)/2, I}] ≥ 0, 
            Return[{z, v2, {}, {}, {{"v2", v2}}}]]; 
If[
            inside[z, {1, 1 + I, (1 + I)/2}] ≥ 0, 
            Return[{z, v3, {}, {}, {{"v3", v3}}}]]; 
If[
            Im[w] ≤ 1/2, 
            Return[{z, 
                e1, {}, {}, {{"e1", e1}}}]]; 
{z, 
            c, {}, {}, {{"c", c}}}];

schmidtstart[w2]

\(schmidtstep[{w_, m_, pictures_, localpix_, schmidtlist_}] := 
      Module[{type, k, zs, rs, pix, zoom, tmatrix, tlist}, zoom = {}; 
        tlist = schmidtlist; 
pix = pictures; 
        type = \((Re[Det[m]] ≠ 0)\); 
        If[type, zoom = circlepix2[w, m]; 
          pix = Join[pix, 
              circlepix2[w, m]], 
 (*\(zoom = 
                trianglepix3[
                  w, \(rvals[m]\)[\([1]\)]];\)*) 
zoom = 
            trianglepix5[w, m]; 
pix = 
            Join[pix, 
              trianglepix3[
                w, \(rvals[
                    m]\)[\([1]\)]]]]; 
 (*\(Print["pix length is ", Length[pix]];\)*) 
Show[
          Graphics[cleanupb[pix]], AspectRatio -> Automatic, 
          PlotRange -> All]; 
circlepix[w, m]; 
        rs = rvals[
            m]; 
If[\(rs[\([1]\)]\)[\([1]\)] == w, 
          Print["w is equal to the first rval of m", {w, \(rs[\([1]\)]\)[\([1]\)]}]; 
          Return[rs]]; 
If[\(rs[\([1]\)]\)[\([2]\)] == w, 
          Print["w is the second rval of m", {w, \(rs[\([1]\)]\)[\([2]\)]}]; 
          Return[rs]]; 
If[\(rs[\([1]\)]\)[\([3]\)] == w, 
          Print["w is the third rval of m", {w, \(rs[\([1]\)]\)[\([3]\)]}]; 
Return[
            rs]]; 
 (*\(Print["type=", 
              type];\)*) 
If[
          type, 
 (*\(Print["Circle case, Det[m] is real"];\)*) 
k = 1; 
tf = 
            schmidttriangle[w, m, 1]; 
If[
            tf,  (*\(Print["k=1,new matrix=", m . v1];\)*) AppendTo[
              tlist, {"v1", 
                v1}]; 
 (*\(circlepix[w, m . v1];\)*) 
              Return[{w, m . v1, pix, zoom, tlist}]]; 
k = 
            2; tf = schmidttriangle[w, m, 2]; 
If[
            tf,  (*\(Print["k=2,new matrix=", m . v2];\)*) AppendTo[
              tlist, {"v2", 
                v2}]; 
 (*\(circlepix[w, 
                  m . v2];\)*) 
Return[{w, m . v2, pix, 
                zoom, tlist}]]; 
k = 3; 
          tf = schmidttriangle[w, m, 3]; 
If[
            tf,  (*\(Print["k=3,new matrix=", m . v3];\)*) AppendTo[
              tlist, {"v3", 
                v3}]; 
 (*\(circlepix[w, 
                  m . v3];\)*) 
Return[{w, m . v3, pix, 
                zoom, tlist}]]; 
k = 4; zs = zset[m, k]; 
          tf = insideschmidttriangle[w, zs]; 
If[
            tf,  (*\(Print["k=4, new matrix=", 
                  m . schmidtmatrix[4]];\)*) 
AppendTo[
              tlist, {"e1", 
                schmidtmatrix[4]}]; 
 (*\(trianglepix3[w, 
                  zs];\)*) 
Return[{w, 
                m . schmidtmatrix[4], pix, zoom, 
                tlist}]]; 
k = 5; 
          zs = zset[m, k]; 
tf = 
            insideschmidttriangle[w, zs]; 
If[
            tf,  (*\(Print["k=5, new matrix=", 
                  m . schmidtmatrix[5]];\)*) 
AppendTo[
              tlist, {"e2", 
                schmidtmatrix[5]}]; 
 (*\(trianglepix3[w, 
                  zs];\)*) 
Return[{w, 
                m . schmidtmatrix[5], pix, zoom, 
                tlist}]]; 
k = 6; 
          zs = zset[m, k]; 
tf = 
            insideschmidttriangle[w, zs]; 
If[
            tf,  (*\(Print["k=6, new matrix=", 
                  m . schmidtmatrix[6]];\)*) 
AppendTo[
              tlist, {"e3", 
                schmidtmatrix[6]}]; 
 (*\(trianglepix3[w, 
                  zs];\)*) 
Return[{w, 
                m . schmidtmatrix[6], pix, zoom, 
                tlist}]]; 
k = 7; 
          zs = zset[m, k]; 
tf = 
            insideschmidttriangle[w, zs]; 
If[
            tf,  (*\(Print["k=7, new matrix=", 
                  m . schmidtmatrix[7]];\)*) 
AppendTo[
              tlist, {"c", 
                schmidtmatrix[7]}]; 
 (*\(trianglepix3[w, 
                  zs];\)*) 
Return[{w, 
                m . schmidtmatrix[7], pix, zoom, 
                tlist}]]; 
Print["w didn't fit anywhere", Return[False]]]; 
If[
          Not[
            type], 
 (*\(Print["triangle type matrix", 
                m];\)*) 
rs = 
            rvals[m]; 
 (*\(trianglepix4[w, 
                m];\)*) 
k = 8; 
          zs = zset[m, 
              k];  (*\(Print["zs=", zs];\)*) 
tf = 
            insideschmidttriangle[w, 
              zs];  (*\(Print["tf=", tf];\)*) 
If[
            tf, 
 (*\(Print["k=8, new matrix=", 
                  m . schmidtmatrix[8]];\)*) 
AppendTo[
              tlist, {"v1", 
                v1}]; 
 (*\(trianglepix3[w, 
                  zs];\)*) 
Return[{w, 
                m . schmidtmatrix[8], pix, zoom, 
                tlist}]]; 
k = 9; 
          zs = zset[m, k]; 
tf = 
            insideschmidttriangle[w, zs]; 
If[
            tf, 
 (*\(Print["k=9, new matrix=", 
                  m . schmidtmatrix[9]];\)*) 
AppendTo[
              tlist, {"v2", 
                v2}]; 
 (*\(trianglepix3[w, 
                  zs];\)*) 
Return[{w, 
                m . schmidtmatrix[9], pix, zoom, 
                tlist}]]; 
k = 10; 
          zs = zset[m, k]; 
tf = 
            insideschmidttriangle[w, zs]; 
If[
            tf, 
 (*\(Print["k=10, new matrix=", 
                  m . schmidtmatrix[10]];\)*) 
AppendTo[
              tlist, {"v3", 
                v3}]; 
 (*\(trianglepix3[w, 
                  zs];\)*) 
Return[{w, 
                m . schmidtmatrix[10], pix, zoom, 
                tlist}]]; 
k = 
            11; 
zs = \(tf = \((inside[w, 
                    rs[\([2]\)]] ≥ 
                  0)\)\); 
If[
            tf, 
 (*\(Print["k=1, new matrix=", 
                  m . c];\)*) 
AppendTo[
              tlist, {"c", 
                c}]; 
 (*\(circlepix[w, 
                  m . c];\)*) 
Return[{w, m . c, pix, zoom, 
                tlist}]]; 
Print["no zone worked, algorithm failed"]; Return[False]]]\)

And now, we execute the algorithm on some examples and see what we get. Outputs have been suppressed, as they tend to be graphics-intensive and thus fatten the file into several extra megs of space.

w = 10^(-100) Floor[10^(100) N[E^I, 105]]

m0 = schmidtstart[w]

Visualization - Deeper Triangles

m1 = schmidtstep[m0];

m10 = Nest[schmidtstep, m0, 10];

Show[Graphics[cleanupb[m10[[3]]]], AspectRatio -> Automatic, 
      PlotRange -> rangy[w, 1/4]]

w2 = 7/11 + 14  I/31

mm0 = schmidtstart[w2]

mm1 = schmidtstep[mm0];

mm2 = schmidtstep[mm1];

mm3 = schmidtstep[mm2];

mm4 = schmidtstep[mm3];

mm5 = schmidtstep[mm4];

mm6 = schmidtstep[mm5];

Show[
      Graphics[cleanupb[mm6[[3]]], AspectRatio -> Automatic, 
        PlotRange -> rangy[w2, 1/3]]]

mm12 = Nest[schmidtstep, mm0, 12];

Show[
      Graphics[cleanupb[mm12[[3]]], AspectRatio -> Automatic, 
        PlotRange -> rangy[w2, 1/200]]]

mm13 = schmidtstep[mm12];

mm14 = schmidtstep[mm13];