(* :Title: framedDGA, Version 1, 7 July 2004 *)

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

(* :Summary:
   This is the noncommutative version of the Mathematica package 
   to calculate framed knot contact homology. 
   This package accompanies `Framed knot contact homology',
   arXiv:math.GT/0407071,
   and can be downloaded from the author's
   homepage at <http://alum.mit.edu/www/ng>.
   It requires installing the NCAlgebra package from
   <http://www.math.ucsd.edu/~ncalg/>.

   Save this file as `framedDGA.m', and then input the package 
   by typing `<<framedDGA.m'.
*)

(* :Context: framedDGA` *)

(* :Mathematica Version: 4.0 *)


Instructions::usage = 
"This program calculates framed knot contact homology invariants. \
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 ?Main for further instructions on the DGA \
and cord algebra invariants, and ?Augment for further instructions on \
augmentation numbers."

Help::usage = 
"Type ?Instructions for instructions on using this package."

Main::usage = 
"\!\(In\ the\ commands\ below, \ \ \(`braid'\ is\ the\ vector\ representing\ a\
\ braid\ whose\ closure\ is\ the\ desired\ knot . \ \ \
\[IndentingNewLine]phimatrix[braid]\  = \ \((\(\[Phi]\_B\) \((A)\))\);\)\ \
\[IndentingNewLine]
  \t\(PhiL[braid]\  = \ \(\[CapitalPhi]\_B\^L\) \((A)\);\)\ \[IndentingNewLine]
  \t\(PhiR[braid]\  = \ \(\[CapitalPhi]\_B\^R\) \((A)\);\)\ \[IndentingNewLine]
  \tdB[braid], dC[braid], \ dD[
        braid], \ dE[braid]\  = \ \ \[PartialD]B, \ \[PartialD]C, \ \[PartialD]D, \ \({\[PartialD]e\_i}\ for\ 
    the\ 0 - framed\ knot\ DGA;\)\ \[IndentingNewLine]
  \td[braid]\  = \ the\ 
          entire\ differential\ for\ the\ 0 - framed\ knot\ DGA, \ written\ \
out; \[IndentingNewLine]\tcordalg[
      braid]\  = \({\[PartialD]B, \[PartialD]C}\  = \ 
          the\ relations\ defining\ the\ cord\ algebra\ of\ the\ \
\(\(knot\)\(.\)\(\)\)\)\)"

Augment::usage = 
"\!\(To\ calculate\ augmentation\ numbers\ for\ a\ knot\ K\ which\ is\ the\ \
closure\ of\ some\ braid, \ first\ enter\ `cordalg[braid]', \ then\ enter\ \
`Aug[%, p, \[Lambda], \[Mu]]' . \ The\ result\ is\ \(\(Aug[K, \[DoubleStruckCapitalZ]\_p, \
\[Lambda], \[Mu]]\)\(.\)\)\)"


<< "C:\\NC\\SetNCPath.m"; 
<< "NCGB.m"; 
<< "IntegerSmithNormalForm.m"; 
Off[General::spell]; 
Off[General::spell1]; 
Print[ "The package framedDGA.m calculates framed knot contact homology. \
For instructions, type ?Instructions." ];
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]; 
phimatrixold[list_] := phimatrixmod[size[list], list]; 
muize[list_] := Table[If[i == j, -1 - \[Mu],
      If[i < j, \[Mu] list[[i,j]], list[[i,j]]]],
       {i, 1, Length[list]}, {j, 1, Length[list]}];
phimatrix[list_] := muize[phimatrixold[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, -1 - \[Mu],
        If[i > j, a[i, j], a[i, j] \[Mu]]],
      {i, 1, n}, {j, 1, n}];
vect[list_, ll_] := Table[If[i == 1, ll, 1], {i, 1, size[list]}];
diag[list_] := Table[If[i == j, list[[i]], 0], {i, 1, 
            Length[list]}, {j, 1, Length[list]}];
phinew[list_] := MatMult[diag[vect[list, \[Lambda]]], PhiL[
        list], trivialA[size[list]], PhiR[list], diag[vect[list, 1/\[Lambda]]]];
dBframed[list_] := trivialA[size[list]] -
       MatMult[diag[vect[list, \[Lambda]]], PhiL[list], trivialA[size[list]]];
dCframed[list_] := 
    trivialA[size[list]] - MatMult[trivialA[size[list]], PhiR[list], \
diag[vect[list, 1/\[Lambda]]]];
trivialB[n_] := Table[b[i, j], {i, 1, n}, {j, 1, n}];
trivialC[n_] := Table[c[i, j], {i, 1, n}, {j, 1, n}];
dDframed[list_] := -MatMult[IdentityMatrix[size[
    list]] - MatMult[diag[vect[list, \[Lambda]]], PhiL[list]], trivialC[size[list]]] \
+ MatMult[trivialB[size[list]], IdentityMatrix[size[list]] - MatMult[
             PhiR[list], diag[vect[list, 1/\[Lambda]]]]];
dEframed[list_] :=
       Table[b[i, i] + MatMult[diag[vect[list, \[Lambda]]], PhiL[list], 
        trivialC[size[list]]][[i, i]], {i, 1, size[list]}];
writhe[list_] := Sum[If[list[[i]] > 0, 1, -1], {i, 1, Length[list]}];
dB[list_] := dBframed[list] /. {\[Lambda] -> \[Lambda] \[Mu]^(-writhe[list])};
dC[list_] := dCframed[list] /. {\[Lambda] -> \[Lambda] \[Mu]^(-writhe[list])};
dD[list_] := dDframed[list] /. {\[Lambda] -> \[Lambda] \[Mu]^(-writhe[list])};
dE[list_] := dEframed[list] /. {\[Lambda] -> \[Lambda] \[Mu]^(-writhe[list])};
cordalg[list_] := {dB[list],dC[list]};
d[list_] := Module[{}, 
    Print["The differential for the 0-framed 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, "] = ", 
        NCExpand[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, "] = ", 
        NCExpand[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, "] = ", 
        NCExpand[dD[list][[i,j]]]]]]; 
     For[i = 1, i < size[list] + 1, i++, 
      Print["\[PartialD] e[", i, "] = ", NCExpand[
        dE[list][[i]]]]]]; 

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_, l_, m_] := 
   expr /. Flatten[{Table[{list[[i]] -> set[[i]]}, 
      {i, 1, Length[list]}], \[Lambda]->l, \[Mu]->m}]; 
cull[n_, expr_, coll_, d_, l_, m_] := 
   Flatten[Table[If[Mod[Numerator[eval[expr, coll[[i]], alist[n], l, m]], 
        d] == 0, {coll[[i]]}, {}], {i, 1, Length[coll]}], 
    1]; 
Auglist[{}, n_, d_, l_, m_] := full[n*(n - 1), d]; 
Auglist[list_, n_, d_, l_, m_] := cull[n, list[[1]], 
    Auglist[trunc[list], n, d, l, m], d, l, m]; 
Augmod[list_, n_, d_, l_, m_] := If[
      {GCD[d, l], GCD[d, m]} == {1, 1},
      Length[Auglist[Flatten[list], n, d, l, m]],
      Print[ "The values for \[Lambda] and \[Mu] 
      must be relatively prime to the modulus." ]];
Aug[list_, p_, l_, m_] := Augmod[list, Length[list[[1]]], p, l, m];