<<  c:/drorbn/projects/KAtlas/KnotTheory.m

EasyKnots = {Knot[4, 1],  Knot[5, 2],  Knot[6, 1],  Knot[7, 2],  
        Knot[8, 1],  Knot[9, 2],  Knot[9, 42]}

Lob[t_] := -NIntegrate[Log[Abs[2  Sin[u]]], {u, 0, t}];  
    Vol[z_] := Plus  @@  (Lob[Arg[#]] &  /@  {z, 1 - 1/z,  
              1/(1 - z)})

2  Vol[omega  =  1/2 + Sqrt[3]/2*I]

PachnerTest[z_, w_]  :=  {Vol[z],  Vol[w],  
        Vol[z] + Vol[w] + Vol[1/(z*w)] - 
          Vol[z (1 - w)/(1 - z*w)] - 
          Vol[w (1 - z)/(1 - z*w)]}

PachnerTest[Random[Real, {-1, 1}] + I*Random[Real, {-1, 1}],  
      Random[Real, {-1, 1}] + I*Random[Real, {-1, 1}]]

SetAttributes[p,  Orderless];  
    SetAttributes[{EQs,  EQ},  {Orderless,  Flat,  
        OneIdentity}]

ClearAll[Rels]

Rels[K_]  :=  Rels[K]  =  
          Rels[PD[K]];
\(Rels[pd_PD]\  := \ 
      Module[
{n, \ dt, z, e, u, v, j, im4, dm4, w, zz, xx, 
          yy, ab, merge, rels, x, y, a, b, c, d, d1, d2, d3, 
          eq}, 
im4[a_]\  := \ a + 1\  /. \ 5 -> 1; \ 
        dm4[a_]\  := \ a - 1\  /. \ 0\  -> 4; 
n = 
          Length[pd]; 
dt\  = \ 
          EQs\  @@ \ 
              Flatten[
{
Table[
t[z[v, j], 3, {1, 2, 4}]~e~
                        t[z[v, im4@j], 4, {2, 1, 3}], 
{v, 
                        n}, {j, 
                        4}
], 
\(Table[
                          e @@ Position[pd, \ i], \ {i, 2  n}]\  /. \ 
                        e[x_, y_]\  -> \ {e[x, y], 
                            e[y, x]}\)\  /. \ {
e[{v_, 
                              a_}, {u_, b_}]\  :> \ 
                          t[z[v, dm4@a], 2, {1, 3, 4}]~e~
                              t[z[u, b], 1, {2, 3, 4}] /; \ 
                            OddQ[a + b], 
e[{v_, a_}, {u_, 
                              b_}]\  :> \ 
                          t[z[v, dm4@a], 2, {1, 3, 4}]~e~
                              t[z[u, b], 1, {2, 4, 3}] /; \ 
                            EvenQ[a + 
                                b]
}
}\  /. \ e[t[z_, a_, {c1_, c2_, c3_}], \ 
                      t[w_, b_, {d1_, d2_, 
                          d3_}]]\  :> \ {
t[z, 
                          a, c1]~p~t[w, b, d1], \ 
                      t[z, a, c2]~p~t[w, b, d2], \ 
                      t[z, a, c3]~p~
                        t[w, b, 
                          d3]
}
]\  /. \ 
            t[z[v_, j_], \ a_, \ b_]\  :> \ 
              t[z[v, j]]~
                  Join~\((t[a, b]\  /. \ {3 -> 4, \ 4 -> 3})\) /; \ 
                EvenQ[j]; 
rels\  = \ 
          List\  @@ \ Union[
\(\(\(dt\  /. \ 
                        p[t[z_, a_, b_], \ t[w_, c_, d_]]\  :> \ 
                          EQs[t[z, a, b]~to~t[w, d, c], \ 
                            t[w, c, d]~to~
                              t[z, b, a]]\)\  //. \ \((
EQs[
                            xx_. *\((t[x_, \ ab__]~to~t[y_, cd__])\), \ 
                            zz_. *\((t[y_, cd__]~to~
                                  t[z_, 
                                    ef__])\)]
 :> \ EQs[xx*y[cd]*
                              zz*\((t[x, ab]~to~
                                  t[z, ef])\)]
)\)\)\  /. \ \((t[z_, ab__]~to~t[z_, ab__])\)\  :> \ z[ab]\)\  /. \ \((w_z)\)[
                      ab__]\  :> \ 
                    Switch[Sort[{ab}], 
{1, 2} | {3, 4}, \ 
                      w[0], 
{1, 3} | {2, 4}, \ \ w[
                        1], 
{1, 4} | {2, 3}, \ 
                      w[2]
]
]\  /. \ 
            eq_Times\  :> \ \((EQ\  @@ \ 
                  eq)\); 
merge\  = \ 
          Position[\(Count[#, \ \(z[n, _]\)[__]] &\)\  /@ \ rels, \ 
            4]; 
rels\  = \ {
Append[
                Delete[rels, merge], \ 
                EQ\  @@ \ rels[\([Join\  @@ \ 
                        merge]\)]], 
zl~
                  EQ~\((EQ\  @@ \ 
                      Product[
\(z[k, 4]\)[0]*\(z[k, 4]\)[
                            1]*\(z[k, 3]\)[0]*\(z[k, 3]\)[
                            2]*
If[
                            pd[\([k, 4]\)] > pd[\([k, 2]\)]\  || \ 
                              pd[\([k, 2]\)] - pd[\([k, 4]\)] > 
                                1, 
\(z[k, 1]\)[
                                1]*\(z[k, 4]\)[
                                2], 
\(z[k, 3]\)[
                                1]*\(z[k, 2]\)[
                                2]
], 
{k, n}
])\)\  //. \ {EQ[z_[0], z_[1], 
                      z_[2]]\  -> \ EQ[\(-1\)], \ 
                  EQ[\(-1\), \(-1\)]\  -> \ 
                    EQ[]}
}\  /. \ \(z[i_, j_]\)[
                k_]\  :> \ 
              m[ToExpression["z" <> ToString[i] <> ToString[j]], \ 
                k]
]\)

rels = Rels[Knot[4, 1]]

$Log  = 20;
RandomPermutation[l_List]  :=  
        Last  /@  Sort[{Random[], #} &  /@  
              l];
m[z_, k1_][k2_]  :=  
        m[z,  (k1 + k2)~Mod~3];
PachnerRules[z1_,  z2_,  z3_,  w1_,  
          w2_]  :=  {
z3  ->  EQ[w1, w2],  
          z2  ->  EQ[w1[2], w2[1]],  
          z1  ->  
            EQ[w1[1],  w2[2]], 
EQ[z1[2],  
              z2[1]]  ->  w1,  EQ[z1[1],  z2[2]]  ->  w2,  
          EQ[z2[2],  z3[1]]  ->  
            w1[1], 
EQ[z3[1],  z1[2]]  ->  w2[1],  
          EQ[z3[2],  z1[1]]  ->  w1[2],  
          EQ[z2[1],  z3[2]]  ->  
            w2[2]
};
\(\(RemoveBigon[s_, \ 
          k_Integer]\  := \((Replace[
            List @@ s[\([1, k]\)], \ {z1_m, 
                z2_m}\  :> \ \((
\(($Log > 
                      50)\) && Print[
                    StringForm["Doing bigon: k=``; s[[1,k]]=``; s=``", \ 
                      k, s[\([1, k]\)], \ s]]; 
t\  = \ 
                  s\  //. \ {EQ[z1, z2]\  -> \ EQ[], \ 
                      EQ[z1[1], \ z2[2]] -> EQ[], \ 
                      EQ[z1[2], \ z2[1]]\  -> 
                        EQ[]}; 
at = \(Position[
                        t[\([1]\)], \ # // First] &\)\  /@ \ {z1, 
                      z2}; 
If[\((Length\  /@ \ 
                        at)\)\  === \ {0, 0}, \ \(($Log > 25)\) && 
                    Print["Annular face removed..."]; \ 
                  Return[t]]; 
If[\((Length\  /@ \ 
                        at)\)\  =!= \ {1, 
                      1}\ , 
\(($Log > 25)\) && 
                    Print[StringForm["Bad bigon: at=``; t=``", \ at, \ 
                        t]]; \ Return[Failed], 
at1 = 
                    at[\([1, 1, 1]\)]; \ 
                  at2 = at[\([2, 1, 
                        1]\)]; 
Return[{
Append[
                        Delete[t[\([1]\)], \ {{at1}, \ {at2}}], \ 
                        EQ[t[\([1, at1]\)], t[\([1, at2]\)]]\  /. \ 
                          EQ[z1, \ z2]\  -> \ 
                            EQ[]], 
\(t[\([2]\)]\  /. \ {z1\  -> \ DeleteCases[t[\([1, at2]\)], \ z2], \ 
                            z2\  -> \ 
                              DeleteCases[t[\([1, at1]\)], \ 
                                z1]}\)\  //. 
{EQ[m[z_, 0], 
                              m[z_, 1], m[z_, 2]]\  -> \ EQ[\(-1\)], \ 
                          EQ[\(-1\), \(-1\)]\  -> \ 
                            EQ[]}
}]
]
)\)]; 
Return[
            Failed]
)\);\)\)
RemoveTriangle[Failed,  _]  =  
        Failed;
\(\(RemoveTriangle[s_, \ 
          k_Integer]\  := \((Replace[
            List @@ s[\([1, k]\)], \ {z1_m, z2_m, \ 
                z3_m}\  :> \ \((
\(($Log > 
                      50)\) && Print[
                    StringForm["Doing triangle: k=``; s[[1,k]]=``; s=``; Lengths = ``", \ k, s[\([1, k]\)], \ s, \ 
                      Length\  /@ \ s[\([1]\)]]]; 
{w1, 
                    w2} = \(m[#, 0] &\)\  /@ \ 
                    Unique[{w, 
                        w}]; 
t\  = \ {
Delete[s[\([1]\)], \ 
                        k]\  //. \ \((rules = 
                          PachnerRules[z1, z2, z3, w1, 
                            w2])\), 
s[\([2]\)]\  //. \ 
                      Join[rules, \ {
EQ[z3[1], \ 
                              z3[2]]\  -> \ 
                            EQ[\(-1\), \ w1[2], w2[1], w1[1], 
                              w2[2]], 
EQ[z2[1], \ 
                              z2[2]]\  -> \ 
                            EQ[\(-1\), \ w1, w2, w1[1], 
                              w2[2]], 
EQ[z1[1], \ 
                              z1[2]]\  -> \ 
                            EQ[\(-1\), \ w1, w2, w1[2], 
                              w2[1]], 
EQ[m[z_, 0], 
                              m[z_, 1], m[z_, 2]]\  -> \ EQ[\(-1\)], \ 
                          EQ[\(-1\), \(-1\)]\  -> \ 
                            EQ[]
}]
}; 
If[
                  FreeQ[t, \ 
                    z123 = First[z1] | First[z2] | 
                        First[z3]], 
Return[
                    t], 
\(($Log > 25)\) && 
                    Print[StringForm["Failed triangle; t=``; problems at ``", \ t, \ Select[t, \ \(! FreeQ[#, \ 
                                z123]\) &]]]
]
)\)]; 
Return[
            Failed]
)\);\)\)
\(\(ContractDomino[s_, \ k1_Integer, \ 
          k2_Integer]\  := \ \((
\(($Log > 70)\) && 
            Print[
              StringForm["{k1, k2} = {``, ``}; s[[1]][[{k1,k2}]]={``,``}", \ k1, 
                k2, \(s[\([1]\)]\)[\([k1]\)], \ \(s[\([1]\)]\)[\([k2]\)]]]; 
Replace[\(s[\([1]\)]\)[\([{k1, 
                  k2}]\)], 
{EQ[_m, \ _m, \ w1_m, w2_m], \ 
                EQ[wa_m, 
                  wb_m, \ _m, \ _m]}\  :> \ \((
\(($Log > 50)\) && 
                    Print[StringForm["Found domino at k1=``, k2=``", \ 
                        k1, \ k2]]; 
{z1, z2, 
                      z3} = \(m[#, 0] &\)\  /@ \ 
                      Unique[{z, z, z}]; 
t = 
                    Insert[s\  //. \ \((Reverse\  /@ \ 
                            PachnerRules[z1, z2, z3, w1, w2])\), \ 
                      EQ[z1, z2, z3], \ {1, \(-1\)}]; 
t = 
                    RemoveTriangle[RemoveTriangle[t, \ k2], \ 
                      k1]; 
\(($Log > 50)\) && \ 
                    Print["t=", \ t]; 
Return[
                    t]
)\) /; \ \((w1 == wa[1]\  && \ 
                    w2 == 
                      wb[2])\)
]; 
Return[
            Failed]
)\);\)\)
Red[rels_List]  :=  
      Module[
{s = rels,  ks,  k,  t,  at,  at1,  
          at2,  k1,  k2}, 
Label[Start];  
        s[[1]] = 
          RandomPermutation[
            DeleteCases[s[[1]],  
              EQ[]]]; 
If[(ev = 
                Complement[Union[Cases[s[[2]],  _Symbol,  Infinity]],  
                  Union[Cases[s[[1]],  _Symbol,  
                      Infinity]]])  =!=  {zl},  
          Print["Extra vars: ",  ev]]; 
Do[
          If[Failed  =!=  (t = RemoveBigon[s,  k]),  s = t;  
            Goto[Start]],  {k,  Length[s[[1]]]}]; 
Do[
          If[Failed  =!=  (t = RemoveTriangle[s,  k]),  s = t;  
            Goto[Start]],  {k,  
            Length[s[[1]]]}]; 
($Log > 50) && 
          Print["Contracting dominos; Lengths=",  
            Length /@  s[[1]],  " s=", 
            s]; 
Do[
If[{4, 
                  4} == (Length  /@  s[[1]][[{k1, 
                          k2}]])  &&  
              Failed  =!=  (t = ContractDomino[s,  k1,  k2]),  
            s = t;  Goto[Start]], 
{k2,  2,  
            Length[s[[1]]]},  {k1,  
            k2 - 1}]; 
s
]

{rels  = Rels[Knot[6, 1]];  s = Red[rels],  
      Length  /@  s[[1]],  Length  /@  s}

vars  =  
      Complement[Union[Cases[s,  _Symbol,  Infinity]],  {zl}]

eqns  =  
      ReplacePart[
(Plus  @@  # == 
                  2  Pi I) &  /@  
            Rest[Flatten[s]]   /.  {zl  ->  0,  
            m[w_, 0]  :>  Log[w],  
            m[w_, 1]  :>  Log[1 - 1/w],  
            m[w_, 2]  :>  
              Log[1/(1 - w)]}, 
4  Pi I,  {-1, 
          2}
]

{
sol  =  
        FindRoot  @@  
          Join[{eqns  /.  zl  -> 1},  ({#,  omega} &  /@  
                vars),  {MaxIterations  ->  
                1000}], 
Max[
        Abs[((First  /@  (eqns  /.  zl  ->  1))  /.  
                sol) - 
            1]], 
Plus  @@  (Vol  /@  (vars  /.  
                sol))
}

briefs  =  
      vars  /.  (sols = 
            NSolve[eqns  /.  zl  ->  1,  vars])

Sort[{Plus  @@  Vol  /@  #,  #} &  /@  
        briefs]

Jeff Weeks: Read about holonomies for CS/Vol computations within SnapPea in gluingequations.c; hyperbolic structure.c

leqns  =  Rest[Flatten[s]] /.  
          eq_EQ  :>  
            Plus  @@  Log  /@  eq   /.  {m[w_, 
              0]  :>  w,  
          m[w_, 1]  :>  1 - 1/w,  
          m[w_, 2]  :>  1/(1 - w),  
          zl  ->  1}

leqns  /.  sols

Im[leqns  /.  sols]  /.  
      arg_Real  :>  (Abs[arg - 2  Pi] < 1)

Position[
      Im[{2  Pi,  0}  /.  
          sols] ,  {(arg_?((Abs[arg - 2  Pi] < 1) &)) .. ,  
        arg_}]

Position[
      Im[leqns  /.  sols]  /.  
        arg_Real  :>  (Abs[arg - 2  Pi] < 
              1),  {True .. ,  False}]

Position[
      Apply[And,  leqns  /.  sols  /.  
          z_?NumberQ  :>  (Im[z] == 2  Pi),  {1}],  
      True]

Plus  @@  Vol  /@  (vars  /.  sols[[13]])