(* :Title: DGAlin, Version 1, 7 November 2003 *)

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

(* :Summary:
   This program calculates linearized homology groups for
   knot 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 `DGAlin.m', and then input the package 
   by typing `<<DGAlin.m'. For instructions on how to use the 
   program, type `?Instructions'.

   This program requires the use of the program
   IntegerSmithNormalForm.m, which can be downloaded from
   the Mathematica Information Center at
   <http://library.wolfram.com/infocenter/MathSource/682>.
*)

(* :Context: DGAlin` *)

(* :Mathematica Version: 4.0 *)



BeginPackage["DGAlin`"]

Instructions::usage =
"\!\(This\ program\ calculates\ linearized\ homology\ groups\ for\n
  knot\ contact\ homology, 
  via\ the\ commands\ HClin[braid]\ and\ HCabellin[braid] . \ 
\[IndentingNewLine]
      Write\ the\ knot\ as\ the\ closure\ of\ a\ braid, \
  and\ abbreviate\ the\ braid\ as\ a\ vector, \(where\ \[Sigma]\_k\ becomes\ \
k\ and\ \[Sigma]\_k\^\(-1\)\ becomes - k;\)\n
  hence, for\ example, 
  the\ braid\ \(\[Sigma]\_2\^\(-1\)\) \(\[Sigma]\_1\^\(-2\)\) \[Sigma]\_3\^3\ \
is\ written\ {\(-2\), \(-1\), \(-1\), 3, 3, 3} . \ 
\[IndentingNewLine]
Then\ HClin[braid]\ outputs\ the\ invariant\ factors\ of\ the\ abelian\ groups\ \
\(HC\_0\^lin\) \((K)\), \(HC\_1\^lin\) \((K)\), \(HC\_2\^lin\) \((K)\), 
  and\ HCabellin[
      braid]\ outputs\ the\ invariant\ factors\ of\ the\ abelian\ groups\ \
\(HC\_0\^\(ab, lin\)\) \((K)\), \(HC\_1\^\(ab, 
        lin\)\) \((K)\), \(HC\_2\^\(ab, lin\)\) \((K)\) . \ 
\[IndentingNewLine]For\ instance, \
  HClin[{1, 1, 1}]\ gives\ {{3}, {3, 3, 3, 0}, {0}}, 
  which\ means\ \(HC\_0\^lin\) \((trefoil)\) = \[DoubleStruckCapitalZ]\_3, \
\(HC\_1\^lin\) \((trefoil)\) = \[DoubleStruckCapitalZ]\_3\[CirclePlus]\
\[DoubleStruckCapitalZ]\_3\[CirclePlus]\[DoubleStruckCapitalZ]\_3\[CirclePlus]\
\[DoubleStruckCapitalZ], \(HC\_2\^lin\) \((trefoil)\) = \(\(\
\[DoubleStruckCapitalZ]\)\(.\)\)\)"


<< "IntegerSmithNormalForm.m"; 
Off[General::spell]; 
Off[General::spell1]; 
matrix12[list_] := Table[Flatten[
     Table[If[j == l, {}, {D[Flatten[list[[1]]][[i]], 
         a[j, l]]}], {j, 1, list[[2]]}, 
      {l, 1, list[[2]]}]], 
    {i, 1, Length[Flatten[list[[1]]]]}]; 
abelianize[n_] := Flatten[
    Table[If[i < j, {a[j, i] -> a[i, j]}, {}], {i, 1, n}, 
     {j, 1, n}]]; 
matrix12abel[list_] := 
   Table[Flatten[Table[Table[D[Flatten[list[[1]]][[i]] /. 
         abelianize[list[[2]]], a[j, l]], 
       {l, j + 1, list[[2]]}], {j, 1, list[[2]] - 1}]], 
    {i, 1, Length[Flatten[list[[1]]]]/2}]; 
HC0[list_] := Flatten[{{decomp[list]}, 
     cokern[matrix11[list]]}, 1]; 
matrix1[list_] := matrix12[matrix11[list]]; 
matrix1abel[list_] := matrix12abel[matrix11[list]]; 
strip[list_] := Flatten[Table[If[list[[i]] == 1, {}, 
      {list[[i]]}], {i, 1, Length[list]}]]; 
H0[list_] := coke[matrix1[list]]; 
matrix22[list_] := Table[Flatten[
     {Table[D[Flatten[list[[1]]][[i]], b[j, l]], 
       {j, 1, list[[2]]}, {l, 1, list[[2]]}], 
      Table[D[Flatten[list[[1]]][[i]], c[j, l]], 
       {j, 1, list[[2]]}, {l, 1, list[[2]]}]}], 
    {i, 1, Length[Flatten[list[[1]]]]}]; 
matrix2[list_] := matrix22[matrix21[list]]; 
abelianizebc[n_] := Flatten[Table[c[i, j] -> b[j, i], 
     {i, 1, n}, {j, 1, n}]]; 
upperdiag[list_] := 
   {Flatten[{Table[If[i < j, list[[1,1,i,j]], {}], 
        {i, 1, list[[2]]}, {j, 1, list[[2]]}], 
       Table[list[[1,2,i]], {i, 1, list[[2]]}]}], 
     list[[2]]} /. abelianizebc[list[[2]]]; 
matrix22abel[list_] := 
   Table[Flatten[{Table[D[list[[1,i]], b[j, l]], 
       {j, 1, list[[2]]}, {l, 1, list[[2]]}]}], 
    {i, 1, Length[list[[1]]]}]; 
matrix2abel[list_] := matrix22abel[
    upperdiag[matrix21[list]]]; 
smith2[array_] := ExtendedSmithForm[array][[2,2]]; 
nullity[diag_] := Sum[Product[If[diag[[i,j]] == 0, 1, 0], 
     {i, 1, Length[diag]}], {j, 1, Length[diag[[1]]]}]; 
rank[diag_] := Sum[Product[If[diag[[i,j]] == 0, 1, 0], 
     {j, 1, Length[diag[[1]]]}], {i, 1, Length[diag]}]; 
kern[smithfirst_, smithsecond_] := 
   Table[smithsecond[[i,j]], {i, 1, Length[smithsecond]}, 
    {j, Length[smithsecond[[1]]] - nullity[smithfirst] + 
      1, Length[smithsecond[[1]]]}]; 
kern[stuff_] := Table[stuff[[2,i,j]], 
    {i, 1, Length[stuff[[2]]]}, 
    {j, Length[stuff[[2,1]]] - nullity[stuff[[1]]] + 1, 
     Length[stuff[[2,1]]]}]; 
prelim[array_] := {SmithForm[array], 
    ExtendedSmithForm[array][[2,2]]}; 
newkernel[list_] := kern[prelim[matrix2[list]]]; 
homology[array_, newarray_] := 
   Table[LinearSolve[newarray, Table[array[[i,j]], 
      {i, 1, Length[array]}]], 
    {j, 1, Length[array[[1]]]}]; 
coke[array_] := Flatten[{strip[InvariantFactors[array]], 
     If[Length[array] < Length[array[[1]]], 
      Table[0, {i, 1, Length[array[[1]]] - 
         Length[array]}], {}]}]; 
Bur[n_, k_] = If[k > 0, Table[If[i == k, 
      If[j == k, 2, If[j == k + 1, -1, 0]], 
      If[i == k + 1, If[j == k, 1, 0], If[j == i, 1, 
        0]]], {i, 1, n}, {j, 1, n}], 
    Inverse[Bur[n, -k]]]; 
newprod[n_, {}] := IdentityMatrix[n]; 
newprod[n_, list_] := Bur[n, list[[1]]] . 
    newprod[n, trunc[list]]; 
new[n_, list_] := Inverse[Matr[n]] . newprod[n, list] . 
    Matr[n]; 
truncate[array_] := Table[array[[i,j]], 
    {i, 1, Length[array] - 1}, 
    {j, 1, Length[array[[1]]] - 1}]; 
Ale[list_] := newprod[size[list], list] - 
    IdentityMatrix[size[list]]; 
deter[list_] := Abs[Det[truncate[Ale[list]]]]; 
decomp[list_] := strip[InvariantFactors[
     truncate[Ale[list]]]]; 
invariant[list_] := {deter[list], decomp[list], H[list]}; 
qlin[i_, j_, m_, amatrix_] := 
   Simplify[If[i == m, If[j == m + 1, amatrix[[j,i]], 
      -amatrix[[i + 1,j]] + 2*amatrix[[i + 1,i]] + 
       2*amatrix[[i,j]]], If[i == m + 1, 
      If[j == m, amatrix[[j,i]], amatrix[[i - 1,j]]], 
      If[i == -m, If[j == -m + 1, amatrix[[j,i]], 
        amatrix[[i + 1,j]]], If[i == -m + 1, 
        If[j == -m, amatrix[[j,i]], -amatrix[[i - 1,j]] + 
          2*amatrix[[i - 1,i]] + 2*amatrix[[i,j]]], 
        If[j == m, -amatrix[[i,j + 1]] + 
          2*amatrix[[i,j]] + 2*amatrix[[j,j + 1]], 
         If[j == m + 1, amatrix[[i,j - 1]], 
          If[j == -m, amatrix[[i,j + 1]], If[j == -m + 1, 
            -amatrix[[i,j - 1]] + 2*amatrix[[i,j]] + 
             2*amatrix[[j,j - 1]], amatrix[[i,
             j]]]]]]]]]]]; 
slin[m_, amatrix_, n_, list_] := 
   Table[If[i == j, 0, qlin[i, j, list[[m]], amatrix]], 
    {i, 1, n}, {j, 1, n}]; 
pmatrixmodlin[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 = slin[k, amatrix, n, list]]; amatrix]; 
size[list_] := Max[Table[Abs[list[[i]]], 
      {i, 1, Length[list]}]] + 1; 
trunc[list_] := Table[list[[i]], {i, 2, Length[list]}]; 
trunclast[list_] := Table[list[[i]], 
    {i, 1, Length[list] - 1}]; 
rev[list_] := Table[list[[Length[list] + 1 - i]], 
    {i, 1, Length[list]}]; 
newarrayR[n_, list_] := If[list[[Length[list]]] > 0, 
    Table[If[i == list[[Length[list]]], 
      If[j == list[[Length[list]]], 
       pmatrixmodlin[n, trunclast[list]][[
        list[[Length[list]]],list[[Length[list]]] + 1]], 
       0], 0], {i, 1, n}, {j, 1, n}], 
    Table[If[i == Abs[list[[Length[list]]]] + 1, 
      If[j == Abs[list[[Length[list]]]] + 1, 
       pmatrixmodlin[n, trunclast[list]][[
        Abs[list[[Length[list]]]] + 1,
        Abs[list[[Length[list]]]]]], 0], 0], {i, 1, n}, 
     {j, 1, n}]]; 
newarrayL[n_, list_] := If[list[[Length[list]]] > 0, 
    Table[If[i == list[[Length[list]]], 
      If[j == list[[Length[list]]], 
       pmatrixmodlin[n, trunclast[list]][[
        list[[Length[list]]] + 1,list[[Length[list]]]]], 
       0], 0], {i, 1, n}, {j, 1, n}], 
    Table[If[i == Abs[list[[Length[list]]]] + 1, 
      If[j == Abs[list[[Length[list]]]] + 1, 
       pmatrixmodlin[n, trunclast[list]][[
        Abs[list[[Length[list]]]],
        Abs[list[[Length[list]]]] + 1]], 0], 0], 
     {i, 1, n}, {j, 1, n}]]; 
PhiRlinear[n_, {}] := Table[0, {i, 1, n}, {j, 1, n}]; 
PhiRlinear[n_, list_] := PhiRlinear[n, trunclast[list]] . 
     Transpose[Bur[n, list[[Length[list]]]]] - 
    Transpose[newprod[n, rev[trunclast[list]]]] . 
     newarrayR[n, list]; 
PhiLlinear[n_, {}] := Table[0, {i, 1, n}, {j, 1, n}]; 
PhiLlinear[n_, list_] := Bur[n, list[[Length[list]]]] . 
     PhiLlinear[n, trunclast[list]] - 
    newarrayL[n, list] . newprod[n, 
      rev[trunclast[list]]]; 
thetamatrix[n_] := Table[1, {i, 1, n}, {j, 1, n}]; 
newAmatrix[n_] := Table[If[i == j, 0, a[i, j]], 
    {i, 1, n}, {j, 1, n}]; 
dblinear[n_, list_] := newAmatrix[n] - 
    newprod[n, rev[list]] . newAmatrix[n] + 
    2*PhiLlinear[n, list] . thetamatrix[n]; 
dclinear[n_, list_] := newAmatrix[n] - 
    newAmatrix[n] . Transpose[newprod[n, rev[list]]] + 
    2*thetamatrix[n] . PhiRlinear[n, list]; 
ddlinear[n_, list_] := 
   Table[b[i, j], {i, 1, n}, {j, 1, n}] . 
     (IdentityMatrix[n] - Transpose[newprod[n, 
        rev[list]]]) - (IdentityMatrix[n] - 
      newprod[n, rev[list]]) . Table[c[i, j], {i, 1, n}, 
      {j, 1, n}]; 
delinear[n_, list_] := 
   Table[b[i, i] + Sum[newprod[n, rev[list]][[i,j]]*
       c[j, i], {j, 1, n}], {i, 1, n}]; 
newprod[n_, {}] := IdentityMatrix[n]; 
newprod[n_, list_] := Bur[n, list[[1]]] . 
    newprod[n, trunc[list]]; 
matrix11[list_] := {{dblinear[size[list], list], 
     dclinear[size[list], list]}, size[list]}; 
matrix21[list_] := {{ddlinear[size[list], list], 
     delinear[size[list], list]}, size[list]}; 
newHC0[list_] := Flatten[{{decomp[list]}, 
     cokern[newmatrix11[list]]}, 1]; 
invar[list_] := {decomp[list], Hnew[list]}; 
almostnew[arrayset_] := {coke[arrayset[[1]]], 
    coke[homology[Transpose[arrayset[[2]]], 
      kern[prelim[Transpose[arrayset[[1]]]]]]], 
    Table[0, {i, 1, rank[SmithForm[arrayset[[2]]]]}]}; 
Hjumblednew[list_] := {almostnew[{matrix1abel[list], 
      matrix2abel[list]}], almostnew[{matrix1[list], 
      matrix2[list]}]}; 
Hnew[list_] := Transpose[Hjumblednew[list]]; 
HClin[list_] := almostnew[{matrix1[list], 
     matrix2[list]}]; 
HCabellin[list_] := almostnew[{matrix1abel[list], 
     matrix2abel[list]}]; 

EndPackage[ ]