(* :Title: DGAabelian, Version 2, 7 November 2003 *)

(* :Author: Lenhard Ng, American Institute of Mathematics
   and Stanford University *)
  

(* :Summary:
   This is the commutative version of the Mathematica package 
   to calculate braid and knot invariants from contact homology. 
   This package accompanies `Knot and braid invariants
   from contact homology I', arXiv:math.GT/0302099, 
   and can be downloaded from the author's
   homepage at <http://alum.mit.edu/www/ng>.

   Save this file as `DGAabelian.m', and then input the package 
   by typing `<<DGAabelian.m'. For instructions on how to use the 
   program, type `?Instructions'.
*)

(* :Context: DGAabelian` *)

(* :Mathematica Version: 4.0 *)



BeginPackage["DGAabelian`"]


Help::usage = 
"Type `?Instructions' for help."

Instructions::usage = 
"This program calculates knot and braid contact homology invariants in the \
commutative setting. Abbreviate a braid as a vector, where \!\(\[Sigma]\_k\) \
becomes k and \!\(\[Sigma]\_k\^\(-1\)\) becomes \
-k; hence, for example, the \
braid \!\(\(\[Sigma]\_2\^\(-1\)\) \(\[Sigma]\_1\^\(-2\)\) \\[Sigma]\_3\^3\) \
is written {-2,-1,-1,3,3,3}. Type ?Braid for further instructions on the \
invariants for braids, ?Knot for further instructions on the invariants for \
knots, or ?Augment for further instructions on augmentation numbers."

Braid::usage = "Here are braid-related commands. Replace `braid' with the \
vector \
representing the desired braid in the following commands:
\[IndentingNewLine]phimatrix[braid] = (\!\(\[Phi]\_B\)(A));
\[IndentingNewLine]dBbraid[braid], dBbraidabel[braid] = \[PartialD]B for the (abelian) braid \
\
DGA;
\[IndentingNewLine]dbraid[braid], dbraidabel[braid] = the entire differential for the (abelian) \
\
braid DGA, written out;
\[IndentingNewLine]HC0braid[braid], HC0braidabel[braid] = the relations defining \!\(HC\_0\)(B), \
\
\!\(HC\_0\^ab\)(B);
\[IndentingNewLine]HC0braidGB[braid], HC0braidabelGB[braid] = a Gr\[ODoubleDot]bner basis for \
the ideal \
defining \!\(HC\_0\)(B), \!\(HC\_0\^ab\)(B)."

\!\(Knot::usage = 
\[IndentingNewLine]\*"\"\<Here are knot-related commands. \
Replace `braid' with the vector representing the desired braid in the \
following commands:
\[IndentingNewLine]phimatrix[braid] = (\\!\\(\[Phi]\\_B\\)(A));
\[IndentingNewLine]PhiL[braid] = \!\(\[CapitalPhi]\_B\^L\)(A);
\[IndentingNewLine]PhiR[braid] = \
\!\(\[CapitalPhi]\_B\^R\)(A);
\[IndentingNewLine]dB[braid], dBabel[braid] = \[PartialD]B for \
the (abelian) knot DGA;
\[IndentingNewLine]dC[braid], dD[braid], dE[braid] = \[PartialD]C, \
\[PartialD]D, {\[PartialD]\!\(e\_i\)} for the knot DGA;
\[IndentingNewLine]d[braid], \
dabel[braid] = the entire differential for the (abelian) knot DGA, written \
out;
\[IndentingNewLine]HC0[braid], HC0abel[braid] = the relations defining \!\(HC\_0\)(K), \
\!\(HC\_0\^ab\)(K);
\[IndentingNewLine]HC0GB[braid], HC0abelGB[braid] = a Gr\[ODoubleDot]bner \
basis for the ideal defining \!\(HC\_0\)(K), \!\(HC\_0\^ab\)(K).\[NewLine]
\[IndentingNewLine]For the \
alternate knot invariants, append `alt' to the end of the relevant command: \
phimatrixalt[braid], HC0GBalt[braid], etc.\>\""\)

\!\(Augment::usage = 
\[IndentingNewLine]\*"\"\<To calculate augmentation \
numbers for a braid B or a knot K, proceed as follows (here n is the number \
of strands in the braid):
\[IndentingNewLine]Aug(B,d): first enter \
`dBbraid[braid]', then enter \
`Aug[%,n,d]';
\[IndentingNewLine]\!\(Aug\^ab\)(B,d): first enter \
`dBbraidabel[braid]', then enter \
`Augabel[%,n,d]';
\[IndentingNewLine]Aug(K,d): first enter `dBC[braid]', then \
enter `Aug[%,n,d]';
\[IndentingNewLine]\!\(Aug\^ab\)(K,d): first enter \
`dBabel[braid]', then enter `Augabel[%,n,d]'.\[NewLine]
\[IndentingNewLine]For the \
alternate knot invariants, append `alt' to the end of the relevant command: \
dBCalt[braid], dBabelalt[braid].\>\""\)



Off[General::spell]; 
Off[General::spell1]; 
trunc[list_] := Table[list[[i]], {i, 2, Length[list]}]; 
size[list_] := Max[Table[Abs[list[[i]]], 
      {i, 1, Length[list]}]] + 1; 
phi[i_, j_, k_, amatrix_] := 
   Simplify[If[i == k, If[j == k + 1, amatrix[[j,i]], 
      -amatrix[[i + 1,j]] - amatrix[[i + 1,i]]*
        amatrix[[i,j]]], If[i == k + 1, 
      If[j == k, amatrix[[j,i]], amatrix[[i - 1,j]]], 
      If[i == -k, If[j == -k + 1, amatrix[[j,i]], 
        amatrix[[i + 1,j]]], If[i == -k + 1, 
        If[j == -k, amatrix[[j,i]], 
         -amatrix[[i - 1,j]] - amatrix[[i - 1,i]]*
           amatrix[[i,j]]], If[j == k, 
         -amatrix[[i,j + 1]] - amatrix[[i,j]]*
           amatrix[[j,j + 1]], If[j == k + 1, 
          amatrix[[i,j - 1]], If[j == -k, amatrix[[i,
            j + 1]], If[j == -k + 1, 
            -amatrix[[i,j - 1]] - amatrix[[i,j]]*
              amatrix[[j,j - 1]], amatrix[[i,
             j]]]]]]]]]]]; 
s[m_, amatrix_, n_, list_] := 
   Table[If[i == j, -2, phi[i, j, list[[m]], amatrix]], 
    {i, 1, n}, {j, 1, n}]; 
phimatrixmod[n_, list_] := 
   Module[{amatrix = Table[If[i == j, -2, a[i, j]], 
       {i, 1, n}, {j, 1, n}]}, 
    For[k = 1, k < Length[list] + 1, k++, 
      amatrix = s[k, amatrix, n, list]]; amatrix]; 
phimatrix[list_] := phimatrixmod[size[list], list]; 
PhiLmod[n_, list_] := Table[phimatrixmod[n + 1, list][[i,
      n + 1]] /. Table[a[m, n + 1] -> If[m == j, 1, 0], 
      {m, 1, n}], {i, 1, n}, {j, 1, n}]; 
PhiRmod[n_, list_] := Table[phimatrixmod[n + 1, list][[
      n + 1,j]] /. Table[a[n + 1, m] -> 
       If[m == i, 1, 0], {m, 1, n}], {i, 1, n}, 
    {j, 1, n}]; 
PhiL[list_] := PhiLmod[size[list], list]; 
PhiR[list_] := PhiRmod[size[list], list]; 
trivialA[n_] := Table[If[i == j, -2, a[i, j]], 
    {i, 1, n}, {j, 1, n}]; 
trivialB[n_] := Table[b[i, j], {i, 1, n}, {j, 1, n}]; 
trivialC[n_] := Table[c[i, j], {i, 1, n}, {j, 1, n}]; 

dBbraid[list_] := trivialA[size[list]] - 
    phimatrix[list]; 
dBbraidabel[list_] := dBbraid[list] /. 
    setabel[size[list]]; 
dbraid[list_] := Module[{}, 
    Print["The differential for the braid DGA of ", 
      list, " is given by:"]; For[i = 1, 
      i < size[list] + 1, i++, For[j = 1, j < i, j++, 
        Print["\[PartialD] b[", i, ",", j, "] = ", 
         Expand[dBbraid[list][[i,j]]]]]; 
       For[j = i + 1, j < size[list] + 1, j++, 
        Print["\[PartialD] b[", i, ",", j, "] = ", 
         Expand[dBbraid[list][[i,j]]]]]]]; 
dbraidabel[list_] := Module[{}, 
    Print["The differential for the abelian braid DGA of ", 
      list, " is given by:"]; For[i = 1, 
      i < size[list] + 1, i++, For[j = i + 1, j < size[list] + 1, j++,
        Print["\[PartialD] b[", i, ",", j, "] = ", 
         Expand[dBbraidabel[list][[i,j]]]]]]]; 

dB[list_] := (IdentityMatrix[size[list]] - PhiL[list]) . 
    trivialA[size[list]]; 
dBabel[list_] := dB[list] /. setabel[size[list]]; 
dC[list_] := trivialA[size[list]] . 
    (IdentityMatrix[size[list]] - PhiR[list]); 
dBC[list_] := {dB[list], dC[list]}; 
dD[list_] := -(IdentityMatrix[size[list]] - 
       PhiL[list]) . trivialC[size[list]] + 
    trivialB[size[list]] . (IdentityMatrix[size[list]] - 
      PhiR[list]); 
dDabel[list_]:= dD[list] /. setabel[size[list]] /.
       setabelbc[size[list]];
dE[list_] := Table[b[i, i] + 
     (PhiL[list] . trivialC[size[list]])[[i,i]], 
    {i, 1, size[list]}]; 
dEabel[list_]:= dE[list] /. setabel[size[list]] /.
       setabelbc[size[list]];
d[list_] := Module[{}, 
    Print["The differential for the knot DGA of ", list, 
      " is given by:"]; For[i = 1, i < size[list] + 1, 
      i++, For[j = 1, j < size[list] + 1, j++, 
       Print["\[PartialD] b[", i, ",", j, "] = ", 
        Expand[dB[list][[i,j]]]]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < size[list] + 1, j++, 
       Print["\[PartialD] c[", i, ",", j, "] = ", 
        Expand[dC[list][[i,j]]]]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < size[list] + 1, j++, 
       Print["\[PartialD] d[", i, ",", j, "] = ", 
        Expand[dD[list][[i,j]]]]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      Print["\[PartialD] e[", i, "] = ", 
       Expand[dE[list][[i]]]]]]; 
dabel[list_] := Module[{}, 
    Print["The differential for the abelian knot DGA of ", list, 
      " is given by:"]; For[i = 1, i < size[list] + 1, 
      i++, For[j = 1, j < size[list] + 1, j++, 
       Print["\[PartialD] b[", i, ",", j, "] = ", 
        Expand[dBabel[list][[i,j]]]]]]; 
     For[i = 1, i < size[list], i++, 
      For[j = i+1, j < size[list] + 1, j++, 
       Print["\[PartialD] d[", i, ",", j, "] = ", 
        Expand[dDabel[list][[i,j]]]]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      Print["\[PartialD] e[", i, "] = ", 
       Expand[dEabel[list][[i]]]]]]; 
HC0braid[list_] := Module[{}, 
    Print["The relations defining \!\(HC\_0\)(B), where \
B is the braid ", list, ", are:"]; 
     For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < i, j++, Print[Expand[
          dBbraid[list][[i,j]]], " = 0"]]; 
       For[j = i + 1, j < size[list] + 1, j++, 
        Print[Expand[dBbraid[list][[i,j]]], " = 0"]]]]; 
HC0[list_] := Module[{}, Print["The relations defining \
\!\(HC\_0\)(K), where K is the closure of ", list, 
      ", are:"]; For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < size[list] + 1, j++, 
       Print[Expand[dB[list][[i,j]]], " = 0"]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < size[list] + 1, j++, 
       Print[Expand[dC[list][[i,j]]], " = 0"]]]]; 
HC0braidGB[list_] := GroebnerBasis[dBbraid[list]]; 
HC0GB[list_] := GroebnerBasis[{dB[list], dC[list]}]; 
setabel[n_] := Flatten[Table[Table[a[j, i] -> a[i, j], 
      {j, i + 1, n}], {i, 1, n - 1}]]; 
setabelbc[n_] := Flatten[Table[c[j,i] -> b[i,j],
      {i, 1, n}, {j, 1, n}]];
HC0braidabel[list_] := Module[{}, 
    Print["The relations defining \!\(HC\_0\^ab\)(B), \
where B is the braid ", list, ", are:"]; 
     For[i = 1, i < size[list] + 1, i++, 
      For[j = i + 1, j < size[list] + 1, j++, 
       Print[Expand[dBbraid[list][[i,j]] /. 
          setabel[size[list]]], " = 0"]]]]; 
HC0abel[list_] := Module[{}, 
    Print["The relations defining \!\(HC\_0\^ab\)(K), \
where K is the closure of ", list, ", are:"]; 
     For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < size[list] + 1, j++, 
       Print[Expand[dB[list][[i,j]] /. 
          setabel[size[list]]], " = 0"]]]]; 
HC0braidabelGB[list_] := GroebnerBasis[
    dBbraid[list] /. setabel[size[list]]]; 
HC0abelGB[list_] := GroebnerBasis[dBabel[list]]; 

alistabel[n_] := Flatten[Table[Table[a[i, j], 
      {j, i + 1, n}], {i, 1, n - 1}]]; 
alist[n_] := Flatten[Table[Table[{a[i, j], a[j, i]}, 
      {j, i + 1, n}], {i, 1, n - 1}]]; 
base[n_, d_] := Table[IntegerDigits[i, d, n], 
    {i, 0, d^n - 1}]; 
shift[coll_, d_] := Table[Table[coll[[i,j]] + 
      Floor[d/2] - d + 1, {j, 1, Length[coll[[1]]]}], 
    {i, 1, Length[coll]}]; 
full[n_, d_] := shift[base[n, d], d]; 
eval[expr_, set_, list_] := 
   expr /. Flatten[Table[{list[[i]] -> set[[i]]}, 
      {i, 1, Length[list]}]]; 
cullabel[n_, expr_, coll_, d_] := 
   Flatten[Table[If[Mod[eval[expr, coll[[i]], 
         alistabel[n]], d] == 0, {coll[[i]]}, {}], 
     {i, 1, Length[coll]}], 1]; 
cull[n_, expr_, coll_, d_] := 
   Flatten[Table[If[Mod[eval[expr, coll[[i]], alist[n]], 
        d] == 0, {coll[[i]]}, {}], 
     {i, 1, Length[coll]}], 1]; 
Augabellist[{}, n_, d_] := full[n*(n - 1)/2, d]; 
Augabellist[list_, n_, d_] := cullabel[n, list[[1]], 
    Augabellist[trunc[list], n, d], d]; 
Auglist[{}, n_, d_] := full[n*(n - 1), d]; 
Auglist[list_, n_, d_] := cull[n, list[[1]], 
    Auglist[trunc[list], n, d], d]; 
Augabel[list_, n_, d_] := 
   Length[Augabellist[Flatten[list], n, d]]; 
Aug[list_, n_, d_] := Length[Auglist[Flatten[list], n, 
     d]]; 


phialt[i_, j_, k_, amatrix_] := 
   Simplify[If[i == k, If[j == k + 1, -amatrix[[j,i]], 
      amatrix[[i + 1,j]] + amatrix[[i + 1,i]]*
        amatrix[[i,j]]], If[i == k + 1, 
      If[j == k, -amatrix[[j,i]], -amatrix[[i - 1,j]]], 
      If[i == -k, If[j == -k + 1, -amatrix[[j,i]], 
        -amatrix[[i + 1,j]]], If[i == -k + 1, 
        If[j == -k, -amatrix[[j,i]], 
         amatrix[[i - 1,j]] - amatrix[[i - 1,i]]*
           amatrix[[i,j]]], If[j == k, 
         amatrix[[i,j + 1]] - amatrix[[i,j]]*
           amatrix[[j,j + 1]], If[j == k + 1, 
          -amatrix[[i,j - 1]], If[j == -k, 
           -amatrix[[i,j + 1]], If[j == -k + 1, 
            amatrix[[i,j - 1]] + amatrix[[i,j]]*
              amatrix[[j,j - 1]], amatrix[[i,
             j]]]]]]]]]]]; 
salt[m_, amatrix_, n_, list_] := 
   Table[If[i == j, 0, phialt[i, j, list[[m]], 
      amatrix]], {i, 1, n}, {j, 1, n}]; 
phimatrixmodalt[n_, list_] := 
   Module[{amatrix = Table[If[i == j, 0, a[i, j]], 
       {i, 1, n}, {j, 1, n}]}, 
    For[k = 1, k < Length[list] + 1, k++, 
      amatrix = salt[k, amatrix, n, list]]; amatrix]; 
phimatrixalt[list_] := phimatrixmodalt[size[list], 
    list]; 
PhiLmodalt[n_, list_] := 
   Table[phimatrixmodalt[n + 1, list][[i,n + 1]] /. 
     Table[a[m, n + 1] -> If[m == j, 1, 0], {m, 1, n}], 
    {i, 1, n}, {j, 1, n}]; 
PhiRmodalt[n_, list_] := 
   Table[phimatrixmodalt[n + 1, list][[n + 1,j]] /. 
     Table[a[n + 1, m] -> If[m == i, 1, 0], {m, 1, n}], 
    {i, 1, n}, {j, 1, n}]; 
PhiLalt[list_] := PhiLmodalt[size[list], list]; 
PhiRalt[list_] := PhiRmodalt[size[list], list]; 
trivialAalt[n_] := Table[If[i == j, 0, a[i, j]], 
    {i, 1, n}, {j, 1, n}]; 

dBalt[list_] := (IdentityMatrix[size[list]] + 
     PhiLalt[list]) . trivialAalt[size[list]]; 
dBabelalt[list_] := dBalt[list] /. 
    setabelalt[size[list]]; 
dCalt[list_] := trivialAalt[size[list]] . 
    (IdentityMatrix[size[list]] + PhiRalt[list]); 
dBCalt[list_] := {dBalt[list], dCalt[list]}; 
dDalt[list_] := -(IdentityMatrix[size[list]] + 
       PhiLalt[list]) . trivialC[size[list]] + 
    trivialB[size[list]] . (IdentityMatrix[size[list]] + 
      PhiRalt[list]); 
dDabelalt[list_]:= dDalt[list] /. setabelalt[size[list]] /.
       setabelbcalt[size[list]];
dEalt[list_] := Table[b[i, i] - 
     (PhiLalt[list] . trivialC[size[list]])[[i,i]], 
    {i, 1, size[list]}]; 
dEabelalt[list_]:= dEalt[list] /. setabelalt[size[list]] /.
       setabelbcalt[size[list]];
dalt[list_] := Module[{}, 
    Print[
      "The differential for the alternate knot DGA of ", 
      list, " is given by:"]; For[i = 1, 
      i < size[list] + 1, i++, For[j = 1, 
       j < size[list] + 1, j++, Print["\[PartialD] b[", i, ",", 
        j, "] = ", Expand[dBalt[list][[i,j]]]]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < size[list] + 1, j++, 
       Print["\[PartialD] c[", i, ",", j, "] = ", 
        Expand[dCalt[list][[i,j]]]]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < size[list] + 1, j++, 
       Print["\[PartialD] d[", i, ",", j, "] = ", 
        Expand[dDalt[list][[i,j]]]]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      Print["\[PartialD] e[", i, "] = ", 
       Expand[dEalt[list][[i]]]]]]; 
dabelalt[list_] := Module[{}, 
    Print["The differential for the alternate abelian knot DGA of ", list, 
      " is given by:"]; For[i = 1, i < size[list] + 1, 
      i++, For[j = 1, j < size[list] + 1, j++, 
       Print["\[PartialD] b[", i, ",", j, "] = ", 
        Expand[dBabelalt[list][[i,j]]]]]]; 
     For[i = 1, i < size[list], i++, 
      For[j = i+1, j < size[list] + 1, j++, 
       Print["\[PartialD] d[", i, ",", j, "] = ", 
        Expand[dDabelalt[list][[i,j]]]]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      Print["\[PartialD] e[", i, "] = ", 
       Expand[dEabelalt[list][[i]]]]]]; 
HC0alt[list_] := Module[{}, Print["The relations \
defining alternate \!\(HC\_0\)(K), where K is the \
closure of ", list, ", are:"]; For[i = 1, 
      i < size[list] + 1, i++, For[j = 1, 
       j < size[list] + 1, j++, 
       Print[Expand[dBalt[list][[i,j]]], " = 0"]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < size[list] + 1, j++, 
       Print[Expand[dCalt[list][[i,j]]], " = 0"]]]]; 
HC0GBalt[list_] := GroebnerBasis[{dBalt[list], 
     dCalt[list]}]; 
setabelalt[n_] := Flatten[
    Table[Table[a[j, i] -> -a[i, j], {j, i + 1, n}], 
     {i, 1, n - 1}]]; 
setabelbcalt[n_] := Flatten[Table[c[j,i] -> -b[i,j],
      {i, 1, n}, {j, 1, n}]];
HC0abelalt[list_] := Module[{}, 
    Print["The relations defining alternate \
\!\(HC\_0\^ab\)(K), where K is the closure of ", list, 
      ", are:"]; For[i = 1, i < size[list] + 1, i++, 
      For[j = 1, j < size[list] + 1, j++, 
       Print[Expand[dBabelalt[list][[i,j]]], " = 0"]]]]; 
HC0abelGBalt[list_] := GroebnerBasis[
    dBabelalt[list]]; 



EndPackage[ ]