The Implicitization of a Trefoil Knot

Mathematica is used to find an implicit description of the trefoil knot.

by Michael Trott

TubeFunctionalShort[
       curve_List       (* parametric representation of
                           the space curve               *),
       cp_              (* curve parameter               *),
       rr_              (* tube radius                   *),
       cpInt_List       (* domain of the curve parameter *),
       ppq_             (* PlotPoints cross section      *),
       ppl_             (* PlotPoints transversal        *),
       opts___          (* possible options for the plot *)] :=
Show[Graphics3D[
Drop[MapThread[Polygon[{#1, #2, #3, #4}]&,
 {#[[1]], #[[2]], RotateRight[#[[3]]], #[[3]]}&[
  {#, RotateRight[#], Transpose[RotateRight[
   Transpose[#]]]}&[
 Function[u1n,
 MapThread[Map[Function[temp, temp + #1], #2]&,
   {Function[cp, Evaluate[N[curve]]] /@ u1n,
    Partition[(Plus @@ #)& /@ Distribute[{
     Transpose[Function[{d1n, d2n},
     {(#/Sqrt[#.#])& /@ MapThread[(#2/(#1.#1) -
      (#1.#2) #1/(#1.#1)^2)&, {d1n, d2n}],
 MapThread[(Function[{t1, t2, t3, k1, k2, k3},
         {t2 k3 - t3 k2, t3 k1 - t1 k3, t1 k2 - t2 k1}] @@
          Flatten[{##}]&), {(#/Sqrt[#.#])& /@
    d1n, (#/Sqrt[#.#])& /@ MapThread[(#2/(#1.#1) -
     (#1.#2) #1/(#1.#1)^2)&, {d1n, d2n}]}]}][
  Function[cp, Evaluate[N[D[curve, cp]]]] /@ u1n,
  Function[cp, Evaluate[N[D[curve, {cp, 2}]]]] /@ u1n]],
    rr N[{Cos[#], Sin[#]}& /@
   N[Range[0, 2Pi(1 - 1/ppq), 2Pi/ppq]]]},
             List, List, List, Times], {ppq}]}]][
       N[Range[cpInt[[1]], cpInt[[2]],
       -Subtract @@ cpInt/ppl]]]]], 2], 1]], opts]

TubeFunctionalShort[
      {Sin[t] + 2 Sin[2 t], Cos[t] - 2 Cos[2 t], Sin[3 t]}
                     (* parametric representation     *),
       t             (* curve parameter               *),
       1/2           (* tube radius                   *),
       {0, 2Pi}      (* domain of the curve parameter *),
       12            (* PlotPoints cross section      *),
       100           (* PlotPoints transversal        *)];

With[{c = {R Sin[t], R Cos[t], 0} 
  (* the center line of a torus *)},
Factor @ (Subtract @@ (Eliminate[# == 0& /@
(({#.#&[{x, y, z} - c] - r^2, 
        D[c, t].({x, y, z} - c)} /.
    {Cos[t] -> (1 - s^2)/(1 + s^2), 
     Sin[t] -> 2 s/(1 + s^2)})), {s}]))]

c = {Sin[t] + 2 Sin[2t], 
     Cos[t] - 2 Cos[2t], Sin[3t]}

{#.#&[{x, y, z} - c] - (1/2)^2, 
         D[c, t].({x, y, z} - c)}

% /. {Cos[3 t] -> Cos[t]^3 - 3 Cos[t] Sin[t]^2,
      Sin[3 t] -> 3 Cos[t]^2 Sin[t] - Sin[t]^3,
      Cos[2 t] -> Cos[t]^2 - Sin[t]^2,
      Sin[2 t] -> 2 Cos[t] Sin[t]}

% /. {Cos[t] -> (1 - s^2)/(1 + s^2), 
      Sin[t] -> 2 s/(1 + s^2)}

{p1, p2} = Numerator /@ Together /@ %

m = p1 (p2 /. s -> alpha) - (p1 /. s -> alpha) p2;

mat1 = CoefficientList[#, alpha]& /@
        CoefficientList[Cancel[m/(s - alpha)], s];

Short[mat1, 6]

mat2 = Map[Replace[#, _Integer rest_ ->  rest]&, 
                                       mat1, {2}];

allElems = Union[Flatten[mat2]];

Length[allElems]

vars = Array[ToExpression[
         "a" <> If[#<10, "0", ""] <> 
               ToString[#]]&, 78];

repRules = Apply[Rule, 
                 Transpose[{allElems, vars}], {1}];

aRules = Dispatch[Reverse /@ repRules];

Short[First[aRules], 8]

mat3 = mat2 //. repRules

minor[mat_, {column_, raw_}] := 
      Delete[#, column]& /@ Delete[mat,raw]

detRule = det[mat_] :>
 Sum[(-1)^(1 + raw) mat[[raw, 
       1]] det[minor[mat, {1, raw}]],
     {raw, 1, Length[mat]}];

LaplaceStep[expr_, j_] :=
 Module[{oldCoeffs, a0, a1, a2, a3, 
         allDets, tempVars, repRules,
         repRulesDis, expa, res, 
         collected, RevRepRules},
 a0 = List @@ expr;
 (* the coefficients from the input det sum *)
 oldCoeffs = a0 /. {_det -> 1};
 (* every coefficient becomes a new 
       variable of the form B[j][k] *)
 oldCoeffsRepRules[j] = 
    Apply[Rule, 
          Transpose[{oldCoeffs, ar = Array[B[j], 
                     Length[oldCoeffs]]}], {1}];
 (* use the abbreviated coefficients *)
 a1 = a0/oldCoeffs;
 (* reduce every det *)
 a2 = a1 /. detRule;
 (* all determinants of the input expression *)
 allDets = Union[Cases[a2, _det, {0, Infinity}]];
 (* form auxiliary variables (in form of powers) *)
 tempVars = Array[C^(#+2)&, Length[allDets]];
 (* for faster collecting we use the 
             C^i powers, not the dets *)
 repRules = 
   Apply[Rule, Transpose[{allDets, tempVars}], {1}];
 (* after a few iterations we have many replacements 
    here, so we dispatch the resulting list *)
 repRulesDis = 
  Dispatch[Apply[Rule, 
           Transpose[{allDets, tempVars}], {1}]];
 (* we collect with respect to the new dets *)
 a3 = a2 /. repRulesDis;
 a4 = a3.ar;
 expa = Expand[a4];
 res = Table[Coefficient[expa, C, i], 
             {i, 0, Exponent[a4, C]}];
 collected = Drop[res, 3].tempVars;
 RevRepRules = Dispatch[Reverse /@ repRules];
 collected /. RevRepRules]

j = 0;
s1 = NestList[(j=j+1;LaplaceStep[#,j])&, 
              det[mat3] /. detRule, 3];

Length[Last[s1]]

Reverse /@ oldCoeffsRepRules[1]

Short[Reverse /@ oldCoeffsRepRules[2], 6]

Short[Reverse /@ oldCoeffsRepRules[3], 6]

Length[Union[Cases[#, _det, 
                   {0,Infinity}]]]& /@ s1

Take[Last[s1], {1, 2}]

j = 0;
s1 = NestList[(j=j+1;LaplaceStep[#,j])&, 
     det[mat3] /. detRule, 11];

j = 0;
s1 = NestList[(j = j + 1;LaplaceStep[#, j])&, 
       det[mat3] /. detRule, 7];

Do[deterSum=s1[[lei+1]];
   (* the prefactor definitions *)
   BRules = Dispatch[Join @@ Table[Reverse /@ 
             oldCoeffsRepRules[i], {i,1,lei}]];
   (* the first determinant *)
   deter = Cases[deterSum[[1]], _det][[1]];
   (* the prefactor itself *)
   preFactor = (deterSum[[1]]/deter);
   (* expand the determinant, result 
      is a sum of products of the as*)
   timeDeter = Timing[(expDeter = 
                  Expand[deter /. det -> Det]);];
   (* replace explicit x, y, z dependent 
      expressions for the as and expand *)
   timeDeterElemExpand=Timing[(d = 
          (Expand[Take[expDeter,1] //. 
                 BRules //. aRules]));];
   (* expand the prefactor *)
   timePrefactorExpand=Timing[
    (p=Expand[preFactor //. BRules //. aRules]);];
   timeProd=Timing[(prod = Expand[d p]);];  
   (* estimating the time it would take to 
      calculate the whole expression *)
   Print[{lei, Length[deterSum] *(
    (r1 = (timeProd[[1]] + 
            timePrefactorExpand[[1]] + 
            timeDeter[[1]])) +
    (r2 = 
     (timeDeterElemExpand[[1]] Length[expDeter])))/
          3600 hours/Second, r2/r1}];, {lei, 3, 7}]

(* all the replacement rules we need here *)
BRules = Dispatch[Join @@ Table[Reverse /@ 
           oldCoeffsRepRules[i], {i, 1, 6}]];
                              
(* the expression to calculate *)
symbDeter = s1[[7]];                            
             
(* the loop for each term of the sum symbDeter *)                            
 
 Do[(* to see some progress *) Print[i];
   Module[{elem, deter, preFactor, 
           expaDeter, d, p, prod},
    elem = symbDeter[[i]];
    deter = Cases[elem, _det][[1]];
    preFactor = elem/deter;
    expaDeter = Expand[deter /. det -> Det];
    d = Expand[expaDeter //. BRules //. aRules]; 
    p = Expand[preFactor //. BRules //. aRules];
    prod = Expand[d p]; 
    Save["directoryAndFileName"<>ToString[i], prod];
    Share[]],
 {i, 1, 792}]

su=0;
Do[Get["directoryAndFileName"<>ToString[i]];
   su = su + prod,
   {i, 1, 792}];
su;

treFoilKnotPoly[x_, y_, z_] = Factor[su];

Block[{x = -1,  y = 2,  z = 3},
      {Resultant[##, s]& @@ {p1, p2}, 
           treFoilKnotPoly[x, y, z]}]

Needs["Graphics`ContourPlot3D`"]

ContourPlot3D[
   Evaluate[treFoilKnotPoly[x, y, z]],
     {x, -4, 4}, {y, -4, 4}, {z, -3/2, 3/2},
     MaxRecursion -> 1, Boxed -> False,
     PlotPoints -> {{21, 4}, {21, 4}, {13, 4}}];

Needs["Graphics`ContourPlot3D`"]

Show[Graphics3D[
(* convert back from polar coordinates 
   to Cartesian coordinates *)
Apply[{#1 Cos[#2], #1 Sin[#2], #3}&,
Cases[ContourPlot3D[Evaluate[
768 x^4 - 1024 x^5 - 128 x^6 + 512 x^7 - 
80 x^8 - 64 x^9 + 16 x^10 + 144 x^2 y^2 - 
768 x^3 y^2 - 136 x^4 y^2 + 896 x^5 y^2 - 
183 x^6 y^2 - 176 x^7 y^2 + 52 x^8 y^2 + 
400 y^4 + 256 x y^4 - 912 x^2 y^4 + 
256 x^3 y^4 + 315 x^4 y^4 - 144 x^5 y^4 - 
16 x^6 y^4 + 4 x^8 y^4 - 904 y^6 - 
128 x y^6 + 859 x^2 y^6 - 16 x^3 y^6 - 
200 x^4 y^6 + 16 x^6 y^6 + 441 y^8 + 
16 x y^8 - 224 x^2 y^8 + 24 x^4 y^8 - 
76 y^10 + 16 x^2 y^10 + 4 y^12 - 
2784 x^3 y z + 4112 x^4 y z - 
968 x^5 y z - 836 x^6 y z + 416 x^7 y z - 
48 x^8 y z + 1312 x y^3 z + 2976 x^2 y^3 z -
5008 x^3 y^3 z - 12 x^4 y^3 z + 
2016 x^5 y^3 z - 616 x^6 y^3 z - 
64 x^7 y^3 z + 32 x^8 y^3 z - 1136 y^5 z - 
4040 x y^5 z + 2484 x^2 y^5 z + 
2784 x^3 y^5 z - 1560 x^4 y^5 z - 
192 x^5 y^5 z + 128 x^6 y^5 z + 1660 y^7 z + 
1184 x y^7 z - 1464 x^2 y^7 z - 
192 x^3 y^7 z + 192 x^4 y^7 z - 472 y^9 z - 
64 x y^9 z + 128 x^2 y^9 z + 32 y^11 z - 
752 x^4 z^2 + 1808 x^5 z^2 - 1468 x^6 z^2 + 
512 x^7 z^2 - 64 x^8 z^2 + 6280 x^2 y^2 z^2 - 
5728 x^3 y^2 z^2 - 4066 x^4 y^2 z^2 +
5088 x^5 y^2 z^2 - 820 x^6 y^2 z^2 - 
384 x^7 y^2 z^2 + 96 x^8 y^2 z^2 -
136 y^4 z^2 - 7536 x y^4 z^2 + 
112 x^2 y^4 z^2 + 8640 x^3 y^4 z^2 -
2652 x^4 y^4 z^2 - 1152 x^5 y^4 z^2 + 
400 x^6 y^4 z^2 + 2710 y^6 z^2 +
4064 x y^6 z^2 - 3100 x^2 y^6 z^2 - 
1152 x^3 y^6 z^2 + 624 x^4 y^6 z^2 -
1204 y^8 z^2 - 384 x y^8 z^2 + 
432 x^2 y^8 z^2 + 112 y^10 z^2 +
3896 x^3 y z^3 - 7108 x^4 y z^3 + 
3072 x^5 y z^3 + 768 x^6 y z^3 -
768 x^7 y z^3 + 128 x^8 y z^3 - 3272 x y^3 z^3 - 
4936 x^2 y^3 z^3 + 8704 x^3 y^3 z^3 - 
80 x^4 y^3 z^3 - 2496 x^5 y^3 z^3 + 
608 x^6 y^3 z^3 + 2172 y^5 z^3 + 
5632 x y^5 z^3 - 2464 x^2 y^5 z^3 - 
2688 x^3 y^5 z^3 + 1056 x^4 y^5 z^3 - 
1616 y^7 z^3 - 960 x y^7 z^3 + 800 x^2 y^7 z^3 +
224 y^9 z^3 + 752 x^4 z^4 - 1792 x^5 z^4 + 
1472 x^6 z^4 - 512 x^7 z^4 + 64 x^8 z^4 - 
3031 x^2 y^2 z^4 + 1936 x^3 y^2 z^4 + 
2700 x^4 y^2 z^4 - 2304 x^5 y^2 z^4 + 
448 x^6 y^2 z^4 + 697 y^4 z^4 + 3728 x y^4 z^4 +
24 x^2 y^4 z^4 - 3072 x^3 y^4 z^4 + 
984 x^4 y^4 z^4 - 1204 y^6 z^4 - 
1280 x y^6 z^4 + 880 x^2 y^6 z^4 + 280 y^8 z^4 - 
800 x^3 y z^5 + 1488 x^4 y z^5 - 768 x^5 y z^5 + 
128 x^6 y z^5 + 992 x y^3 z^5 + 
1016 x^2 y^3 z^5 - 1728 x^3 y^3 z^5 + 
480 x^4 y^3 z^5 - 472 y^5 z^5 - 960 x y^5 z^5 + 
576 x^2 y^5 z^5 + 224 y^7 z^5 + 16 x^4 z^6 +
388 x^2 y^2 z^6 - 384 x^3 y^2 z^6 + 
96 x^4 y^2 z^6 - 76 y^4 z^6 - 384 x y^4 z^6 + 
208 x^2 y^4 z^6 + 112 y^6 z^6 - 64 x y^3 z^7 +
32 x^2 y^3 z^7 + 32 y^5 z^7 + 4 y^4 z^8 /.
 (* to polar coordinates *)  
{x -> r Cos[p], y -> r Sin[p]}],
{r, 0.6, 3.3}, {p, 0, 2Pi}, {z, -1.3, 1.3}, 
PlotPoints -> {12, 32, 16},
MaxRecursion -> 0, DisplayFunction -> Identity],
                   _Polygon, Infinity], {-2}]]];