[Math] Random knot on six vertices

knot-theorypr.probabilitystick-knots

This question is inspired by Joseph O'Rourke's beautiful question on random knots. Choose an random ordered 6-tuple of points on the unit sphere in $\mathbf{R}^3$, and form a knot by connecting successive pairs of points in the 6-tuple by sticks (see the picture at Joseph's question). By known results on stick numbers, the resulting knot will either be the unknot or the trefoil knot. What is the probability of producing one or the other?

Best Answer

I wrote a program in Mathematica to sample knots from this distribution and test what proportion are the trefoil knot.

In order to tell if a given knot is the unknot or the trefoil, the program first checks the total curvature of the knot and applies the Fary-Milnor theorem: if the curvature is less than $4 \pi$, then it's the unknot. Half the time, this test identifies the unknot. I think it should be possible to compute the exact probability of the curvature being too small.

Next, the program projects the knot onto 100 random planes. If any of these projections has less than 3 crossings, then we are again considering the unknot. This test eliminates all but ~1% of cases.

Finally, if we're still not done, the program takes the projection with the least number of crossings and checks if the resulting knot diagram is tricolorable. Usually this diagram has three crossings and this test might be a bit of a sledgehammer, but this test completely distinguishes the unknot from the trefoil. (I don't use this test first because my implementation is very slow.)

In a test run of 10,000 random knots, 68 knots were determined to be the trefoil. The computation took about 12 minutes. Here's one of the trefoils it found:

An HSV-colored trefoil

The code follows. As usual, beware of bugs.

(* Random points, projections, those sorts of things *)
randsph[] := Normalize@Table[RandomVariate@NormalDistribution[], {3}]
randknot[] := Table[randsph[], {6}]
close[x_] := Join[x, {First[x]}]
project[ x_, frame_ ] := Flatten[frame[[2 ;; 3]] . Transpose[ {x} ]]
framify[x_] := Orthogonalize@{x, randsph[], randsph[]}
rotate[{x_, y_}] := {-y, x}
halfintersecthelper[a_, b_, c_, 
  d_] := (a - c) . rotate[b - a] / ((d - c) . rotate[b - a])
halfintersect[a_, b_, c_, d_] := 
 0 <= halfintersecthelper[a, b, c, d] <= 1
intersect[a_, b_, c_, d_] := 
 halfintersect[a, b, c, d] && halfintersect[c, d, a, b]
nintshelper[cknot3_, frame_] := 
 Module[{cknot2 = (project[#1, frame] &) /@ cknot3}, 
  Table[If[Abs[i - j] > 1 && Abs[i - j] != 5 && 
     intersect[cknot2[[i]], cknot2[[i + 1]], cknot2[[j]], 
      cknot2[[j + 1]]], {i, 
     halfintersecthelper[cknot2[[j]], cknot2[[j + 1]], cknot2[[i]], 
      cknot2[[i + 1]]], 
     If[over[cknot3[[i]], cknot3[[i + 1]], cknot3[[j]], 
       cknot3[[j + 1]], frame], +1, -1], {Min[i, j], Max[i, j]}}, {0, 
     0, 0, 0}], {i, 1, 6}, {j, 1, 6}]]
nints[cknot3_, frame_] := (#1[[3 ;; 4]] &) /@ 
  Union[Select[Flatten[nintshelper[cknot3, frame], 1], #1[[3]] != 0 &]]
curvature[cknot3_] := 
 Total@Table[
   VectorAngle[cknot3[[i + 1]] - cknot3[[i]], 
    cknot3[[1 + Mod[i + 1, 6]]] - cknot3[[i + 1]]], {i, 1, 6}]
overhelper[a_, b_, c_, d_] := (b - a)\[Cross](d - c)
over[a_, b_, c_, d_, frame_] := 
 overhelper[a, b, c, d].(c - a) overhelper[a, b, c, d].frame[[1]] > 0

(* Can this knot be tricolored? *)
vars[seq_] := x /@ Range@Length@seq
domains[xs_] := And @@ (#1 == 0 || #1 == 1 || #1 == 2 &) /@ xs
nonconstant[seq_] := ! 
  And @@ Table[x[i] == x[i + 1], {i, 1, Length[seq] - 1}]
overs[seq_] := 
 And @@ Module[{n = Length[seq]}, 
   Table[If[seq[[i, 1]] == +1, x[i] == x[1 + Mod[i, n]], True], {i, 1,
      n}]]
names[seq_] := Union[(#1[[2]] &) /@ seq]
overname[seq_, n_] := 
 x@First@Flatten[Position[seq, {+1, n}, {1}, Heads -> False]]
undername1[seq_, n_] := 
 x@First@Flatten[Position[seq, {-1, n}, {1}, Heads -> False]]
undername2[seq_, n_] := 
 x[1 + Mod[First@Flatten[Position[seq, {-1, n}, {1}, Heads -> False]],
     Length[seq]]]
overunder[seq_, n_] := 
 Mod[overname[seq, n] + undername1[seq, n] + undername2[seq, n], 
   3] == 0
overunders[seq_] := And @@ (overunder[seq, #1] &) /@ names@seq
conditions[seq_] := 
 domains[vars@seq] && overs@seq && overunders@seq && nonconstant@seq
tricolor[seq_] := FindInstance[conditions@seq, vars@seq]

(* Init *)
overalltrials = 0;
overallcount = 0;

(* Random trials! *)
First@
 Timing@Module[{trials = 10000, nframes = 100, count = 0, frames, i, 
    j, k, crossings, ncrossings, pgood, projn, projj},
   frames = framify /@ Table[randsph[], {nframes}];
   For[i = 1, i <= trials, i++,
    k = close[randknot[]];
    (* Angles *)
    pgood = If[curvature[k] >= 4 Pi, 0, -1];
    (* Projections *)
    projn = 20;
    projj = 0;
    For[j = 1, j <= nframes && pgood == 0, j++,
     crossings = nints[k, frames[[j]] ];
     ncrossings = Length@crossings/2;
     If[ncrossings < 3, pgood = -1];
     If[ncrossings < projn, projn = ncrossings; projj = j];
     ];
    If[pgood == 0, crossings = nints[k, frames[[projj]]];
     pgood = If[tricolor@crossings != {}, +1, -1];];
    (* Record *)
    If[pgood == +1 && count == 0, testk = k; 
     testf = frames[[projj]]];
    If[pgood == +1, count++];
    ];
   overalltrials += trials;
   overallcount += count;
   ]
overallcount
overalltrials
overallcount / overalltrials * 100. 

(* Draw a trefoil knot found by the random trials *)
nints[testk, testf] // MatrixForm
pk = Map[project[#, testf] &, testk ]
Graphics3D[{Thickness[0.02], Opacity[1], Specularity[White, 50], 
 Line[testk, VertexColors -> {Red, Yellow, Green, Cyan, Blue, Purple, 
 Red}]}, Axes -> False, PlotRange -> All, Boxed -> False]
Related Question