BeginPackage["HFhat`"]  

Print["Version: 2006.09.14. Loading the package HFhat`..."]

(*** variables ***)
OriginDiagram; NiceDiagram;
CFhat;

(*** functions ***)
ComputeChainComplex;
GetHFhat;   GetHFKhat;

Begin["`Private`"]

(*******************************************************
	The main function
*******************************************************)

ComputeChainComplex[hd_]:= Module[
	{},
	Clear[HeeGenus, OriginDiagram, NiceDiagram];
	Clear[CFhat, CFhatGraph, HFhat, HFhatGraph];
	Clear[CFKhat, CFKhatGraph, HFKhat, HFKhatGraph];
	Clear[KnotGenus, KnotFibred];
	Clear[DS, RS, VS];
	
	Get["LinearAlgebra`MatrixManipulation`"];		(* we need this package for matrix operation *)

	GetInitData[hd];								(* Initialize the data *)
	GetNiceDiagram[];								(* get a nice Heegaard diagram *)
	GetCFhat[];										(* get the chain complex *)

	Print["The chain complex is computed successfully"];
]

GetHFhat[spc_] := Module[
	{RHS, FGS, Diff, GS, HoloDisk, Homology, SP, FGS1, Diff1, parents},

	{RHS, FGS, Diff, GS, HoloDisk, Homology, SP} = CFhat;
	
	If [Length[Select[FGS, (#[[1]]==spc)&]]==0, Print["No homology in this Spin^c structure"]; Return[]];
	Diff = Select[Diff, (FGS[[ #[[1]], 1]]==spc && FGS[[ #[[2]],1 ]]==spc)&];
	
	FGS = Transpose[ Append[Transpose[FGS], Table[1,{Length[FGS]}]]];
	
	While[True,
		If[Length[Diff]==0, Break[]];

		(*Print["Kill", Diff[[1]]];*)
		
		{father, child} = Diff[[1]];
		FGS[[ father, 5]] = 0;	
		FGS[[ child,  5]] = 0;
		Diff = Rest[Diff];
		
		NewFathers  = Select[Diff, (#[[2]]==child )&];
		NewChildren = Select[Diff, (#[[1]]==father)&];
		Diff = Select[Diff, (#[[1]]!=father && #[[1]]!=child && #[[2]]!=father && #[[2]]!=child)&];
		If [Length[NewFathers]>0 && Length[NewChildren]>0,
			NewFathers  = Transpose[NewFathers] [[1]];
			NewChildren = Transpose[NewChildren][[2]];
			NewDiff = CartesianProduct[NewFathers, NewChildren];
			For[i=1, i<=Length[NewDiff], i++,
				If[ MemberQ[Diff, NewDiff[[i]]],
					Diff=Drop[Diff,{Position[Diff, NewDiff[[i]]][[1,1]]}]
					,
					Diff=Append[Diff, NewDiff[[i]]];
				];
			];
		];
		
		Survived = Select[FGS, (#[[1]]==spc && ##[[5]]==1)&];
		(*Print[Survived, Diff];
		Print["Remaining", {Length[Survived], Length[Diff]}];*)
		
		(** check the validity **)
		For[i=1, i<=Length[Diff], i++,
			If[ FGS[[Diff[[i,1]],5]]==0, Print[ Diff[[i,1]], " not survived "]; Return[]];
		];
	];
	
	Survived = Select[FGS, (#[[1]]==spc && ##[[5]]==1)&];
	minGrading = Min[Transpose[Survived][[2]]];
	maxGrading = Max[Transpose[Survived][[2]]];
	
	HeeHomology = Table[ {minGrading+n-1, 0}, {n, 1, maxGrading-minGrading+1}];
	For[i=1, i<=Length[Survived], i++,
		HeeHomology[[ Survived[[i,2]]-minGrading+1, 2 ]]++;
	];
	Print[HeeHomology //MatrixForm];
]

GetHFKhat[spc_]:=Module[
	{RHS, FGS, Diff, GS, HoloDisk, Homology, SP, FGS1, Diff1, parents},

	{RHS, FGS, Diff, GS, HoloDisk, Homology, SP} = CFhat;
	
	If [Length[Select[FGS, (#[[1]]==spc)&]]==0, Print["No homology in this Spin^c structure"]; Return[]];
	Diff = Select[Diff, (FGS[[ #[[1]], 1]]==spc && FGS[[ #[[2]],1 ]]==spc && FGS[[#[[1]],3]]==FGS[[#[[2]],3]])&];
	
	FGS = Transpose[ Append[Transpose[FGS], Table[1,{Length[FGS]}]]];
	
	While[True,
		If[Length[Diff]==0, Break[]];

		(*Print["Kill", Diff[[1]]];*)
		
		{father, child} = Diff[[1]];
		FGS[[ father, 5]] = 0;	
		FGS[[ child,  5]] = 0;
		Diff = Rest[Diff];
		
		NewFathers  = Select[Diff, (#[[2]]==child )&];
		NewChildren = Select[Diff, (#[[1]]==father)&];
		Diff = Select[Diff, (#[[1]]!=father && #[[1]]!=child && #[[2]]!=father && #[[2]]!=child)&];
		If [Length[NewFathers]>0 && Length[NewChildren]>0,
			NewFathers  = Transpose[NewFathers] [[1]];
			NewChildren = Transpose[NewChildren][[2]];
			NewDiff = CartesianProduct[NewFathers, NewChildren];
			For[i=1, i<=Length[NewDiff], i++,
				If[ MemberQ[Diff, NewDiff[[i]]],
					Diff=Drop[Diff,{Position[Diff, NewDiff[[i]]][[1,1]]}]
					,
					Diff=Append[Diff, NewDiff[[i]]];
				];
			];
		];
		
		Survived = Select[FGS, (#[[1]]==spc && ##[[5]]==1)&];
		(*Print[Survived, Diff];
		Print["Remaining", {Length[Survived], Length[Diff]}];*)
		
		(** check the validity **)
		For[i=1, i<=Length[Diff], i++,
			If[ FGS[[Diff[[i,1]],5]]==0, Print[ Diff[[i,1]], " not survived "]; Return[]];
		];
	];
	
	Survived = Select[FGS, (#[[1]]==spc && ##[[5]]==1)&];
	minGrading = Min[Transpose[Survived][[2]]];
	maxGrading = Max[Transpose[Survived][[2]]];
	minFiltration = Min[Transpose[Survived][[3]]];
	maxFiltration = Max[Transpose[Survived][[3]]];
	
	HeeHomology = Table[ 0, {maxGrading-minGrading+1}, {maxFiltration-minFiltration+1}];
	For[i=1, i<=Length[Survived], i++,
		HeeHomology[[ (maxGrading-Survived[[i,2]]+1), (Survived[[i,3]]-minGrading+1) ]]++;
	];
	Print[HeeHomology //MatrixForm];
]

(******************************************************
	General functions
******************************************************)

(*** return the Cartesian product of two lists ***)
CartesianProduct[a_List, b_List] := Flatten[Outer[List, a, b, 1, 1], 1]

(*** give the mod function i Mod(n), if zero, return n instead of zero***)
ModUp[i_, n_] := If[Mod[i,n]==0,n,Mod[i,n]]

(** Swap two elements of two lists **)
SwapList[L_, i_, j_] := Module[
	{newL},
	newL = L;
	newL[[i]] = L[[j]];
	newL[[j]] = L[[i]];
	newL
]

(** Smith form for a matrix **)
GetSmithForm[oM_] := Module[
	{i, j, k, M, m, n, L, R},

	M = oM;

	(** check whether this is a matrix **)
	If[MatrixQ[M]==False, Print["Are you kidding? Sorry, I accept a matrix only!!"]; Abort[]];

	(** size of M **)
	{m, n} = Dimensions[M];
	sizeM = Min[m, n];

	(** check whether this is an integer matrix **)	
	If[ Count[ Flatten[M], _Integer] != m*n, Print["Do not fool me. Give me an integer matrix"]; Abort[]];

	(** Initialize L & R **)
	L = Table[If[i==j, 1, 0], {i, 1, m}, {j, 1, m}];
	R = Table[If[i==j, 1, 0], {i, 1, n}, {j, 1, n}];
	
	(** First make the matrix a diagonal one **)
	For[i=1, i<=sizeM, i++,
		(** check whether the lower right submatrix starting (i,i) is zero already. If so, we are done **)
		subM = LinearAlgebra`MatrixManipulation`TakeMatrix[M, {i,i}, {m,n}];
		If [ Count[subM, 0, 2] == (m-i+1)*(n-i+1), Break[] ];
		
		(** if M[[i,i]] is zero, try to find a nonzero term and switch that to [[i,i]] **)
		If [ M[[i,i]]==0,
			For[j=i, j<=m, j++,
				For[k=i, k<=n, k++,
					If[ M[[j,k]] != 0, Break[]];
				];
				If[ k<n+1, Break[] ];
			];
			M = SwapList[Transpose[M], i, k];				(** exchange i-th column and k-th column **)
			R = Transpose[SwapList[Transpose[R], i, k]];	(** do the same for R **)
			
			M = SwapList[Transpose[M], i, j];				(** exchange i-th row and j-th row **)
			L = SwapList[L, i, j];							(** L will exchange i-th column and j-th column**)
		];

		(** who like negative stuff? **)
		If [ M[[i,i]]<0, 
			M[[i]] = -M[[i]];					(** multiply i-th row by -1 **)
			L[[i]] = -L[[i]];					(** do the same for L **)
		];

		(** while there is any entry M[[j,i]] in the i-th column not divided by M[[i,i]] **)
		For [j=i+1, j<=m, j++,
			While [ True,
				r = Mod[M[[j,i]], M[[i,i]]];
				q = ( M[[j,i]] - r ) / M[[i,i]];
				M[[j]] = M[[j]] - q * M[[i]];					(** add (-q) x (i-th row) to j-th row **)
				L[[j]] = L[[j]] - q * L[[i]];					(** do the same for L **)
				
				If[r==0, Break[]];

				M = SwapList[M, i, j];							(** exchange i-th row and j-th row **)
				L = SwapList[L, i, j];							(** do the same for L **)
			];
		];
		
		(** while there is any entry M[[i, j]] in the i-th row not divided by M[[i,k]] **)
		M = Transpose[M]; R = Transpose[R];
		For [j=i+1, j<=n, j++,
			While [ True,
				r = Mod[M[[j,i]], M[[i,i]]];
				q = ( M[[j,i]] - r ) / M[[i,i]];
				M[[j]] = M[[j]] - q * M[[i]];					(** add (-q) x (i-th column) to j-th column **)
				R[[j]] = R[[j]] - q * R[[i]];					(** do the same for R **)
				
				If[r==0, Break[]];

				M = SwapList[M, i, j];							(** exchange the i-th column and the j-th column **)
				R = SwapList[R, i, j];							(** do the same to R **)
			];
		];
		M=Transpose[M]; R=Transpose[R];
	];

	(*Print[L //MatrixForm, oM//MatrixForm, R//MatrixForm, M // MatrixForm, L . oM . R //MatrixForm];*)
	M = Join[Table[M[[i,i]],{i,1,Min[m,n]}], If[n>m, Table[0, {n-m}], {}]];

	{L, R, M}
]

(*******************************************************
	Initialize Data: regions, vertices, 
*******************************************************)

GetInitData[hd_] := Module[
	{HeeGenus=0, DS={}, Dnum=0, ES={}, Enum=0, AS={}, BS={}, CS={}, VS={}, Vnum=0, HeeDist=0,
		i, j, k, r, s, tmp={}
	},

	DS = hd; Dnum = Length[DS];
	
	(** Check validity of the Heegaard diagram **)
	For[i=1, i<=Length[DS], i++,
		If[Mod[Length[DS[[i]]],2]!=0, Print["The disk ", i, " has an odd number of edges"]; Exit[]]
	];
	tmp = Flatten[DS];
	Enum = Length[tmp]/2;

	(** edges should be indexed from one to the number of edges **)
	If[Enum != Max[DS], Print["Please index the edges from 1 to the number of edges"]; Exit];

	(** each edge should appears twice, one postive and one negative **)
	For[i=1, i<=Enum, i++,
		j = Count[DS, +i, 2];
		k = Count[DS, -i, 2];
		If[j!=1, Print["The edge ", +i, " appears ", j, " times"]; Exit];
		If[k!=1, Print["The edge ", -i, " appears ", k, " times"]; Exit];
	];

	(** edges (left disk, right disk, type, index **)
	ES = Table[{0, 0, 0, 0},{Enum}];
	For[i=1, i<=Dnum, i++,
		For[j=1, j<=Length[DS[[i]]], j++,
			k = DS[[i,j]];
			If[k>0, ES[[+k, 1]]=i, ES[[-k, 2]]=i];
			ES[[Abs[k],3]] = ModUp[j,2];
		]
	];

	(** Find all alpha curves **)
	r=1;		(** j=index for curve, k=disk on the left **) 
	For[i=1, i<=Enum, i++,
		If[ES[[i,3]]==2 || ES[[i,4]]>0, Continue[]];
		ES[[i,4]]=r;
		AS=Append[AS, {i}];
		j=i;
		s = ES[[j, 1]];
		k = ModUp[Position[DS[[s]], j][[1,1]]+1, Length[DS[[s]]]];
		j = DS[[s, k]];
		s = If[j>0, ES[[j,2]], ES[[Abs[j],1]]];
		k = ModUp[Position[DS[[s]], -j][[1,1]]+1, Length[DS[[s]]]];
		j = DS[[s, k]];
		While[j!=i,
			AS[[r]]=Append[AS[[r]], j];
			ES[[j, 4]]=r;
			s = ES[[j, 1]];
			k = ModUp[Position[DS[[s]], j][[1,1]]+1, Length[DS[[s]]]];
			j = DS[[s, k]];
			s = If[j>0, ES[[j,2]], ES[[Abs[j],1]]];
			k = ModUp[Position[DS[[s]], -j][[1,1]]+1, Length[DS[[s]]]];
			j = DS[[s, k]]
		];
		r++
	];

	HeeGenus=r-1;

	(** Find all beta curves **)
	r=1;		(** j=index for curve, k=disk on the left **) 
	For[i=1, i<=Enum, i++,
		If[ES[[i,3]]==1 || ES[[i,4]]>0, Continue[]];
		ES[[i,4]]=r;
		BS=Append[BS, {i}];
		j=i;
		s = ES[[j, 1]];
		k = ModUp[Position[DS[[s]], j][[1,1]]+1, Length[DS[[s]]]];
		j = DS[[s, k]];
		s = If[j>0, ES[[j,2]], ES[[Abs[j],1]]];
		k = ModUp[Position[DS[[s]], -j][[1,1]]+1, Length[DS[[s]]]];
		j = DS[[s, k]];
		While[j!=i,
			BS[[r]]=Append[BS[[r]], j];
			ES[[j, 4]]=r;
			s = ES[[j, 1]];
			k = ModUp[Position[DS[[s]], j][[1,1]]+1, Length[DS[[s]]]];
			j = DS[[s, k]];
			s = If[j>0, ES[[j,2]], ES[[Abs[j],1]]];
			k = ModUp[Position[DS[[s]], -j][[1,1]]+1, Length[DS[[s]]]];
			j = DS[[s, k]]
		];
		r++
	];

	r--;

	(* Print["Genus=", HeeGenus, "    AS=", AS//MatrixForm, "    BS=", BS//MatrixForm];*)

	(** The number of alpha curves and the number of beta curves must be the same **)
	If[r!=HeeGenus, Print["There are ", HeeGenus, " alpha curves and ", r, " beta curves"]; Exit[]];

	(** 
		check whether the Heegaard diagram represents a rational homology sphere 
		But we are able to deal with any (orientable) three-manifold.
	**)
	BA=Table[0, {HeeGenus}, {HeeGenus}];
	For[i=1, i<=HeeGenus, i++,
		For[j=1, j<=Length[AS[[i]]], j++,
			s = ES[[ AS[[i,j]], 1]];
			k = DS[[ s, ModUp[ Position[DS[[s]], AS[[i,j]]][[1,1]]+1, Length[DS[[s]]] ] ]];
			BA[[ i, ES[[ Abs[k], 4]] ]] += Sign[k];
		];
	];
	(* If[Det[BA]==0, Print["Not a rational homology sphere"]; Abort[]]; *)
	(* If[Abs[Det[BA]]==1, (* integer homology sphere *) Print["integer homology sphere"]]; *)
	
	(** add length and badness disks (distance, badness, list of edges) **)
	For[i=1, i<=Dnum, i++, DS[[i]] = {-1, Length[DS[[i]]]/2, DS[[i]]} ];

	(** Compute the distance of each disk **)
	For[DS[[1,1]]=0; HeeDist=0, Count[Table[DS[[i,1]],{i,1,Dnum}], -1]>0, HeeDist++,
		For[i=1, i<=Enum, i++,    
			If[ES[[i,3]]==1, Continue[]];				(** an alpha curve, continue **)
			If[DS[[ES[[i,1]],1]]==HeeDist && DS[[ES[[i,2]],1]]==-1,
				DS[[ES[[i,2]],1]]=HeeDist+1
			];
			If[DS[[ES[[i,2]],1]]==HeeDist && DS[[ES[[i,1]],1]]==-1,
				DS[[ES[[i,1]],1]]=HeeDist+1
			];
		];
	];
	
	OriginDiagram = {HeeGenus, Dnum, Enum, Vnum, AS, BS, CS, DS, ES, VS};
	
	Print["All datas initialized"];
]
		
(*******************************************************
	Make the Heegaard Diagram nice
*******************************************************)

GetNiceDiagram[] := Module[
	{HeeGenus, DS, Dnum, ES={}, Enum=0, AS={}, BS={}, CS={}, VS={}, Vnum=0, HeeDist=0,
		betastar, diskm, diskstar, badm, edgehead, edgeheadpos, diskhead, finger, edgetail, edgetailpos,
		signstar,
		i, j, k, r, s, tmp={}
	},

	(** 
		Get the data from OriginDiagram,
		What we use: Dnum, Enum, HeeGenus, DS, ES.
		What we do not use: AS, BS, CS, VS, Vnum 
	**)
	{HeeGenus, Dnum, Enum, Vnum, AS, BS, CS, DS, ES, VS}=OriginDiagram;

	(** Get the largest distance of bad disks **)
	HeeDist = Max[ Transpose[Select[DS, (#[[2]]>2)&]] [[1]] ];

	(* Print["DS=", Transpose[Prepend[Transpose[DS], Range[Length[DS]]]]//MatrixForm, 
		"     ES=", Transpose[Prepend[Transpose[ES], Range[Length[ES]]]] //MatrixForm]; *)

	While[HeeDist>0,
		(** Find the distance HeeDist disk (diskm) to start the algorithm **)
		For[badm=Enum; i=1, i<=Length[DS], i++,
			If[ DS[[i,1]]==HeeDist && DS[[i,2]]<badm && DS[[i,2]]>2, badm=DS[[i,2]]; diskm=i];
		];
		
		(** Find betastar and diskstar **)
		For[i=2, i<=Length[DS[[diskm,3]]], i+=2,
		    betastar = DS[[diskm, 3, i]];
			diskstar = If[betastar>0, ES[[betastar, 2]], ES[[-betastar, 1]]];
			If[DS[[diskstar, 1]]==HeeDist-1, Break[];];
		];

		(** find our finger **)
		edgeheadpos = Position[DS[[diskstar, 3]], -betastar][[1,1]];
		finger = {{diskstar, -betastar, edgeheadpos}};
		edgeheadpos = ModUp[i+3, Length[DS[[diskm,3]]] ];
		edgehead = DS[[diskm, 3, edgeheadpos]];
		finger = Append[finger, {diskm, edgehead, edgeheadpos}];
		While[True,
			diskhead = If[edgehead>0, ES[[edgehead, 2]], ES[[-edgehead, 1]] ];
			
			(** if we reach a bigon or a region with smaller distance **)
			If[ DS[[diskhead,1]]<HeeDist || DS[[diskhead,2]]==1, Break[];];

			(** if we reach another bad region of distance Heedist **)
			If[ DS[[diskhead,1]]==HeeDist && DS[[diskhead, 2]]>2 && diskhead!=diskm, Break[]];

			(** If we come back to diskm **)
			If[ diskhead==diskm,
				edgeheadpos = Position[DS[[diskm, 3]], -edgehead][[1,1]];

				(** if we come back via the next (adjacent) edge, we are done **)
				If [ Mod[ edgeheadpos - finger[[2, 3]] - 2, Length[DS[[diskm,3]]]]==0, Break[]];
				
				(** if we come back via the previous (adjacent) edge, we are done, but we need to reverse the finger **)
				If [ Mod[ edgeheadpos - finger[[2, 3]] + 2, Length[DS[[diskm,3]]]]==0, 
					(** start to reverse the finger **)
					
					Break[];
				];
				
				(** 
					if we are not coming back via an adjacent edge, we start the finger over from 
					the beta edge next to DS[[diskm, 3, finger[[2,3]]]]
				**)
				finger = Take[finger, 2];
				edgeheadpos = ModUp[finger[[2, 3]]+2, Length[DS[[diskm, 3]]]];
				edgehead = DS[[diskm, 3, edgeheadpos]];
				finger[[2]]={diskm, edgehead, edgeheadpos};
				Continue[];
			];

			(** now we reach a square region of distance >= Heedist **)
			edgeheadpos = ModUp[
				Position[DS[[diskhead,3]], -edgehead][[1,1]]+2, 
				Length[ DS[[diskhead,3]] ] 
			];
			edgehead = DS[[diskhead, 3, edgeheadpos]];
			finger = Append[finger, {diskhead, edgehead, edgeheadpos}];
		];

		(* print our finger *)
		(*Print["diskm=", diskm, " badm=", badm, " diskstar=", diskstar, " betastar=", betastar]; 
		Print["finger=", TableForm[finger,TableDepth->1]]; *)

		If[ diskhead==diskm,
		(** handle slide **)
			Print["This part will be written recently, please semd me this example so that I could test this part. Thanks."];
			Abort[]
			,
		(** finger move **)
			signstar = Sign[finger[[1,2]]];  (* sign of betastar, we need to use a lot *)
			betastarcurve = ES[[Abs[betastar], 4]];

			(** modify disk (diskstar) **) 
			DS[[diskstar,2]] = DS[[diskstar,2]]+1;  			(* badness increase by one *)
			DS[[diskstar,3]] = Insert[ DS[[diskstar,3]], Sign[finger[[2,2]]]*(Enum+2), finger[[1, 3]] ];
			DS[[diskstar,3]] = Insert[ DS[[diskstar,3]], signstar*(Enum+1), finger[[1, 3]] ];
				
			(** add beta edges (Enum+1)**)
			ES = If [signstar>0,
				Append[ ES, {diskstar, Dnum+1, 2, betastarcurve} ], 
				Append[ ES, {Dnum+1, diskstar, 2, betastarcurve} ]
			];
				
			(** add alpha edges (Enum+2) and (Enum+3) **)
			edgehead = finger[[2, 2]];
			If [edgehead>0,
				ES = Append[ ES, {diskstar, Dnum+2, 1, ES[[Abs[edgehead], 4]]} ];
				ES = Append[ ES, {Dnum+1,   Dnum+3, 1, ES[[Abs[edgehead], 4]]} ]
				,
				ES = Append[ ES, {Dnum+2, diskstar, 1, ES[[Abs[edgehead], 4]]} ];
				ES = Append[ ES, {Dnum+3, Dnum+1,   1, ES[[Abs[edgehead], 4]]} ]
			];

			(** add disk (Dnum+1) **)
			DS[[diskm, 3]] = RotateLeft[ DS[[diskm,3]], Position[DS[[diskm,3]], betastar][[1,1]] ];
			s = Position[DS[[diskm,3]], finger[[2,2]]][[1,1]];
			DS = Append[DS, {HeeDist+1, (s+1)/2,
				Join[ {Sign[edgehead]*(Enum+3), Sign[betastar]*(Enum+1)}, Take[DS[[diskm,3]], s-1] ]
			}];
			
			(** modify the edges of the disk (Dnum+1) **)
			For[i=3, i<=s+1, i++,
				r = Last[DS][[3, i]];
				ES[[Abs[r], If[r>0, 1, 2] ]] =Dnum+1
			];

			(** modify disk (diskm) **)
			DS[[diskm]] = {HeeDist, (Length[DS[[diskm,3]]]-s+1)/2, Drop[DS[[diskm,3]],s-1]};
	
			For[i=3, i<=Length[finger], i++,
				(** add beta edges (Enum+4*i-8) and (Enum+4*i-7) **)
				If[ signstar>0,
					ES = Append[ES, {Dnum+2*i-4, finger[[i,1]], 2, betastarcurve}];
					ES = Append[ES, {Dnum+2*i-4, Dnum+2*i-3,    2, betastarcurve}]
					,
					ES = Append[ES, {finger[[i,1]], Dnum+2*i-4, 2, betastarcurve}];
					ES = Append[ES, {Dnum+2*i-3,    Dnum+2*i-4, 2, betastarcurve}]
				];
					
				(** add alpha edges (Enum+4*i-6) and (Enum+4*i-5) **)
				r = ES[[ Abs[ finger[[i,2]] ], 4]];
				If[ finger[[i,2]]>0,
					ES = Append[ES, {Dnum+2*i-4, Dnum+2*i-2, 1, r}];
					ES = Append[ES, {Dnum+2*i-3, Dnum+2*i-1, 1, r}]
					,
					ES = Append[ES, {Dnum+2*i-2, Dnum+2*i-4, 1, r}];
					ES = Append[ES, {Dnum+2*i-1, Dnum+2*i-3, 1, r}]
				];
					
				(** add disks (Dnum+2*i-4) and disks (Dnum+2*i-3) **)
				DS = Append[DS, {HeeDist, 2, {
					(Enum+4*i-10) * (-Sign[finger[[i-1,2]]]), 
					(Enum+4*i-7) * signstar, 
					(Enum+4*i-6) * Sign[finger[[i,2]]], 
					(Enum+4*i-8) * signstar
				}}];
				r = DS[[finger[[i,1]], 3, ModUp[finger[[i,3]]+3,4] ]]; 	(* the right edge to be modified*)
				DS = Append[DS, {HeeDist, 2, {
					(Enum+4*i-9) * (-Sign[finger[[i-1,2]]]), 
					r, 
					(Enum+4*i-5) * Sign[finger[[i,2]]], 
					(Enum+4*i-7) * (-signstar)
				}}];
				
				(** modify the edge (r) **)
				If[r>0, ES[[r,1]]=Dnum+2*i-3, ES[[-r,2]]=Dnum+2*i-3];

				(** modify disk (finger[[i,1]]) **)					
				DS[[ finger[[i,1]], 3, ModUp[finger[[i,3]]+3,4] ]] = (-signstar) * (Enum+4*i-8);
			];

			(** add disk (Dnum+2*i-2) **)
			i = Length[finger];
			DS = Append[DS, {DS[[diskhead,1]]+1, 1, {(-Sign[Last[finger][[2]]])*(Enum+4i-6), signstar * (Enum+4*i-4)}}];
				
			(** add beta edge (Enum+4*i-4) **)
			ES = If[ signstar>0,
				Append[ES, {Dnum+2*i-2, diskhead, 2, betastarcurve}],
				Append[ES, {diskhead, Dnum+2*i-2, 2, betastarcurve}]
			];
				
			(** modify alpha edge (Enum+4*i-5) **);
			ES[[Enum+4*i-5, If[finger[[i,2]]>0,2,1] ]] = diskhead;

			(** modify disk (diskhead) **)
			DS[[diskhead,2]] = DS[[diskhead,2]]+1;			(* badness increase by one *)
			DS[[diskhead,3]] = Insert[ DS[[diskhead,3]], 
				Unevaluated[Sequence@@{
					(Enum+4*i-4) * (-signstar), 
					(Enum+4*i-5) * (-Sign[finger[[i,2]]])
				}],
				Position[DS[[diskhead,3]], -finger[[i,2]]][[1,1]] +1
			];
			
			Dnum = Length[DS];		(* Dnum = Dnum + 2*(i-1) *)
			Enum = Length[ES];		(* Enum = Enum + 4*(i-1) *)
		];
			
		(** re-compute the distance and loop **)
		HeeDist = Max[Transpose[Select[DS, (#[[2]]>2)&]] [[1]]];
		
		(*Print["DS=", Transpose[Prepend[Transpose[DS], Range[Length[DS]]]]//MatrixForm, 
			"     ES=", Transpose[Prepend[Transpose[ES], Range[Length[ES]]]] //MatrixForm];*)
		
		(* check validity of resultin diagram *)
		For[i=1, i<=Length[DS], i++,
			For[j=1, j<=Length[DS[[i,3]]], j++,
				k = DS[[i,3,j]];
				If[k>0, If[ ES[[k,1]]!= i, 
					Print["The left disk of edge ", k, " should be ", i];
					If[Input["Input 1 to abort"]==1, Abort[]];
				]];
				If[k<0, If[ ES[[-k,2]]!= i, 
					Print["The left disk of edge ", k, " should be ", i];
					If[Input["Input 1 to abort"]==1, Abort[]];
				]];
				
			];
		];
		tmp = Flatten[Transpose[DS][[3]]];
		If[ Length[tmp] != 2 * Length[ES],
			Print["There are ", Length[ES], " edges, while the number of edges in DS is ", Length[tmp]]; Abort[]
		];
		For[i=1, i<=Length[ES], i++,
			If[ Count[tmp, +i] != 1, Print["The edge ", +i, " appears in DS ", Count[tmp, +i], " times"]; Abort[] ];
			If[ Count[tmp, -i] != 1, Print["The edge ", -i, " appears in DS ", Count[tmp, -i], " times"]; Abort[] ];
		];
	];
	
	(* Print["DS=", Transpose[Prepend[Transpose[DS], Range[Length[DS]]]]//MatrixForm, 
		"     ES=", Transpose[Prepend[Transpose[ES], Range[Length[ES]]]] //MatrixForm]; *)

	NiceDiagram = {HeeGenus, Dnum, Enum, Vnum, AS, BS, CS, DS, ES, VS};
	
	Print["Nice diagram achieved"];
]

(*******************************************************
	Read the chain complex CFhat
	{{grading, index, alphaindex, alphapoint, betaindex, betapoint},
	 {totalindex, totalindex}}
*******************************************************)
    
(*** Get the chain complex ***)
GetCFhat[] := Module[
	{HeeGenus, DS, Dnum, ES={}, Enum=0, AS={}, BS={}, CS={}, VS={}, Vnum=0, HeeDist=0,
		GS={},
		i, j, k, r, s, tDS={}, tES={}, finger={},
		SF={}, RHS, FGS, HoloDisk, Diff, IndexZeroDisk, Homology, SP
	},
	{HeeGenus, Dnum, Enum, Vnum, AS, BS, CS, DS, ES, VS}=NiceDiagram;
	AS=BS=CS=Table[{},{n,HeeGenus}];
	VS=GS={};
	
	(** Find all alpha curves **)
	For[i=1, i<=HeeGenus, i++,
		j = r = Position[Transpose[Drop[Transpose[ES],2]], {1,i}][[1,1]];	(* some edge belonging to \alpha_i *)
		AS[[i]] = Append[AS[[i]], r];										(* add to AS[[i]] *)
		While[True,
			s = ES[[j, 1]]; 																	(* (s) = left disk of edge (r) *)
			j = DS[[s, 3, ModUp[Position[DS[[s,3]], j][[1,1]]+1, Length[DS[[s,3]]]] ]];			(* next edge of (j) on disk (s) *)
			s = ES[[Abs[j], If[j>0, 2, 1] ]];													(* the disk on the other side of (j) *)
			j = DS[[s, 3, ModUp[Position[DS[[s,3]], -j][[1,1]]+1, Length[DS[[s,3]]]] ]];		(* the edge continuing \alhpa_i *)
			
			If[ j==r, Break[]];				(* if we got a loop, we are done with \alpha_i *)
			AS[[i]]=Append[AS[[i]], j];		(* if not, add to \alpha_i *)
		];
	];
	
	(** Find all beta curves **)
	For[i=1, i<=HeeGenus, i++,
		j = r = Position[Transpose[Drop[Transpose[ES],2]], {2,i}][[1,1]];		(* some edge (r) belonging to \beta_i *)
		BS[[i]] = Append[BS[[i]], r];
		While[True,
			s = ES[[j, 1]]; 													(* the left disk of edge (r) *)
			j = DS[[s, 3, ModUp[Position[DS[[s,3]], j][[1,1]]+1, Length[DS[[s,3]]]] ]];
			s = ES[[Abs[j], If[j>0, 2, 1] ]];
			j = DS[[s, 3, ModUp[Position[DS[[s,3]], -j][[1,1]]+1, Length[DS[[s,3]]]] ]];
			
			If[ j==r, Break[]];
			BS[[i]]=Append[BS[[i]], j];
		];
	];
	
	(** check whether the curves match what they are assigned in ES **)
	(*For[ i=1, i<=HeeGenus, i++,
		For[ j=1, j<=Length[AS], j++,
			If [ ES[[AS[[i,j]],3]] != 1 or ES[[AS[[i,j]],4]] != i, 
				Print["Somethin wrong with edge", AS[[i,j]]]; Abort[]
			];
		];
	];
	For[ i=1, i<=HeeGenus, i++,
		For[ j=1, j<=Length[BS], j++,
			If [ ES[[BS[[i,j]],3]] != 2 or ES[[BS[[i,j]],4]] != i, 
				Print["Somethin wrong with edge", BS[[i,j]]]; Abort[];
			];
		];
	];*)

	(** 
		Find our gamma curves-- CS[[i]] corresponds to AS[[i]]
		our stragey to find \gamma_i is: find two disks separated by an \alpha_i edge, now try to find
		a shortest path from one disk to the other in the complement of alpha and gamma curves,
		So what we do is to find the distance of the other disk from the fist one
	**)
	tES = Transpose[ Append[ Transpose[ES], Range[Length[ES]]]];
	For[tDS=DS; i=1, i<=HeeGenus, i++,
		alphaedge = AS[[i, 1]];							(* some edge on \alpha_i *)
		{disk1, disk2} = Take[tES[[alphaedge]], 2];		(* disk1=left disk, disk2=right disk *)

		(** define distances with respect to (disk1) **)
		tDS = Transpose[Join[{Table[-1, {Length[tDS]}]}, Drop[Transpose[tDS],1]]];		(* first set all distances to be -1 *)
		tDS[[disk1,1]]=0; 																(* set the distance of (disk1) to be 0 *)
		For[r=0, Count[Transpose[tDS][[1]], -1]>0 && tDS[[disk2,1]]==-1, r++,			(* r will record the distance *)
			For[j=1, j<=Length[tES], j++,
				If[tES[[j,3]]!=2, Continue[]];	(** If this is an alpha, continue **)
				If[tDS[[tES[[j,1]],1]]==r && tDS[[tES[[j,2]],1]]==-1, tDS[[tES[[j,2]],1]]=r+1 ];
				If[tDS[[tES[[j,2]],1]]==r && tDS[[tES[[j,1]],1]]==-1, tDS[[tES[[j,1]],1]]=r+1 ];
			];
		];
 		
		(** now find a shortest path from disk2 to disk1, as like a finger move**)
		finger = {{disk1, alphaedge, Position[tDS[[disk1,3]], alphaedge][[1,1]]}};
		CS[[i]] = {alphaedge};
		For[disk=disk2; edge=alphaedge, r>0, r--,
			(** find an adjacent disk of distance (r-1) **)
			tDS[[disk,3]] = RotateLeft[ tDS[[disk,3]], Position[ tDS[[disk,3]], -edge][[1,1]]-1];
			For[j=1, True, j++,
				edge = tDS[[disk, 3, j]];
				If[tES[[Abs[edge],3]]!=2, Continue[]];
				s = tES[[Abs[edge], If[edge>0, 2, 1]]];
				If[ tDS[[s,1]]==r-1, Break[]];
			];
			finger = Append[finger, {disk, edge, Position[ tDS[[disk, 3]], edge][[1,1]]} ];
			CS[[i]] = Append[CS[[i]], Sign[edge] * tES[[Abs[edge], 5]]];
			disk = s;
		];
		tDS[[disk,3]] = RotateLeft[ tDS[[disk,3]], Position[ tDS[[disk,3]], -edge][[1,1]]-1];
		finger[[1,3]] = Position[ tDS[[disk1, 3]], alphaedge][[1,1]];
		
		(** now we separate the disks using this gamma curve (to make sure the gamma curves are pairwise disjoint) **)
		tEnum = Length[tES];
		tDnum = Length[tDS];
		L     = Length[finger];
		For[j=1, j<=L, j++,
			k = ModUp[j-1,L];
			
			(** add the disk (tDnum + j) **)
			tDS = Append[tDS, {0, 0, 
				Join[{
					(tEnum + 2*k) * (-Sign[finger[[k, 2]]]),
					(tEnum + 2*j -1),
					(tEnum + 2*j) * Sign[finger[[j,2]]]
				}, Drop[tDS[[finger[[j,1]],3]], finger[[j,3]]]]
			}];
			
			(** modify edges of old edges of disk (tDnum+j) **)
			For[l=4, l<=Length[tDS[[tDnum+j,3]]], l++,
				edge = tDS[[tDnum+j, 3, l]];
				If[edge>0, tES[[edge,1]]=tDnum+j, tES[[-edge, 2]]=tDnum+j];
			];

			(** modify disk finger[[j,1]] **)
			tDS[[finger[[j,1]],3]] = Append[ Take[ tDS[[finger[[j,1]],3]], finger[[j,3]]], -(tEnum+2*j-1)];
			
			(** add edges (tEnum+2j-1) and (tEnum+2j) **)
			tES = Append[tES, {tDnum+j, finger[[j,1]], 3, i, 0} ];			(* edge (tEnum+2j-1) *)
			tES = If[ finger[[j,2]]>0,
				Append[tES, Join[{tDnum+j, tDnum+ModUp[j+1,L]}, Drop[tES[[Abs[finger[[j,2]]]]],2] ]],
				Append[tES, Join[{tDnum+ModUp[j+1,L], tDnum+j}, Drop[tES[[Abs[finger[[j,2]]]]],2] ]]
			];
		]; 
	];
	
	(*Print["alpha=", AS//MatrixForm]; Print["beta=",  BS//MatrixForm]; Print["gamma=", CS//MatrixForm]; *)	
	
	(** rewrite the local intersection number of edges with gamma curves, 
		for later use of get intersection numbers of a loop with gamma curves
	**)
	EC = Table[{}, {Length[ES]}];
	For[i=1, i<=HeeGenus, i++,
		For[j=1, j<= Length[CS[[i]]], j++,
			r = CS[[i,j]];
			EC[[ Abs[r] ]] = Append[ EC[[Abs[r]]], {i, Sign[r]}];
		];
	];
	(* Print["EC=", EC//MatrixForm]; *)

	(** find intersection points (of alpha and beta curves) **)
	VS=Table[{}, {HeeGenus}, {HeeGenus}];
	For[i=1, i<=HeeGenus, i++,
		For[j=1, j<=Length[AS[[i]]], j++,
			s = ES[[AS[[i,j]],1]];
			k = DS[[s, 3, ModUp[Position[DS[[s,3]], AS[[i,j]]][[1,1]]+1, Length[DS[[s,3]]] ] ]];
			If[k<0,	k=-k, 
				s = ES[[k, 4]]; k = BS[[ s, ModUp[Position[ BS[[s]], k][[1,1]]-1, Length[BS[[s]]]] ]];
			];
			VS[[ i, ES[[k, 4]]]] = Append[ VS[[ i, ES[[k, 4]]]], {AS[[i,j]], k}];
		];
	];
	(*Print[VS];Print[Map[Length, VS, {2}]//MatrixForm];*)

	(** get all generators {{alphaedge_1, betaedge_1}, {alphaedge_2, betaedge_2}, ...}**)
	If[HeeGenus==1, GS=Map[{#}&, Flatten[VS,2]],
		PBs = Permutations[Range[HeeGenus]];
		PBlen = Factorial[HeeGenus];
		For[GS={};i=1, i<=PBlen, i++,
			PB = PBs[[i]];
			SG = CartesianProduct[VS[[1, PB[[1]]]], VS[[2, PB[[2]]]]];
			For[j=3, (j<=HeeGenus) && (Length[SG]>0), j++, SG = CartesianProduct[SG, VS[[j, PB[[j]]]]] ];
			If[Length[SG]>0, GS=Join[GS, Map[Flatten, SG]] ]
		];
		GS = Map[Partition[#,2]&, GS];
	];
	(*Print["There are ", Length[GS], " generators"];
	Print[ Transpose[Prepend[Transpose[GS], Range[Length[GS]]]] //MatrixForm]; *)
	
	(** EToG Maps Edges to Generators (only alpha edges) **)
	EToG = Table[{}, {Length[ES]}];
	For[ i=1, i<=Length[GS], i++,
		For[j=1, j<=HeeGenus, j++,
			k = GS[[i,j,1]];
			EToG[[k]] = Append[ EToG[[k]], i];
		];
	];
	(*Print[EToG];*)
	
	(** 
		Compute the homology of the three-manifold.  The homology is generated by the gamma curves, 
		and each beta curve provide a relation (read as the intersection number of \beta_i with \alpha_j)
	**)
	BA = Table[0, {HeeGenus}, {HeeGenus}];		(** the intersection number of \beta & \alpha, beta = BA . gamma **)
	For[i=1, i<=HeeGenus, i++,
		For[j=1, j<=Length[AS[[i]]], j++,
			s = ES[[ AS[[i,j]], 1]];
			k = DS[[ s, 3, ModUp[ Position[DS[[s,3]], AS[[i,j]]][[1,1]]+1, Length[DS[[s,3]]] ] ]];	(* k is an beta edge *)
			BA[[ ES[[ Abs[k], 4]], i ]] += - Sign[k];
		];
	];
	BC = Table[0, {HeeGenus}, {HeeGenus}];		(** the intersection number of \beta & \gamma, beta = BC . alpha **)	
	For[i=1, i<=HeeGenus, i++,
		For[j=2, j<=Length[CS[[i]]], j++,
			BC[[ ES[[ Abs[CS[[i, j]]], 4]], i ]] += Sign[ CS[[i, j]] ]; 
		];
	];
	HomDiff = LinearAlgebra`MatrixManipulation`BlockMatrix[{
		{Table[ If[i==j, 1, 0], {i, 1, HeeGenus}, {j, 1, 2* HeeGenus}]},
		{BC, BA}
	}];
	{SFL, SFR, Homology} = GetSmithForm[HomDiff];			(** SFR=Change of Basis **)
	RHS = (Length[Select[Homology, (#==0)&]]==0);			(* where the three manifold is a rational homology sphere *)
	(*Print[SFL//MatrixForm, HomDiff//MatrixForm, SFR//MatrixForm, SFL . HomDiff . SFR //MatrixForm, Homology//MatrixForm];*)

	(** 
		Now we compute the grading, spinc, filtration of our chain complex, 
		we set the grading of the first generator to be zero.
		This is where we need the 3-manifold to be a rational homology sphere.
		use the loop to go along the direction of the alpha curve or beta curve.
	**)
	FGS = Table[{0,0,0,i}, {i,1,Length[GS]}];		(* spinc, grading, filtration, index of generator*)
	FGS[[1, 1]] = 1;								(* Set the first generator *)
	SP  = {Table[0, {2*HeeGenus}]};					(* SP = Set of spinc structures *)
	(** now assign spinc structures **)
	For[i=2, i<=Length[GS], i++,
		(** find the loop from GS[[1]] to GS[[i]], the part on alpha curves are already given **)
		LoopB = Table[{0,0}, {HeeGenus}];		(** LoopB **)
		For [j=1, j<=HeeGenus, j++,
			LoopB[[ ES[[GS[[1, j, 2]],4]], 1]] = GS[[1,j,2]];
			LoopB[[ ES[[GS[[i, j, 2]],4]], 2]] = GS[[i,j,2]];
		];
		LA = LC = Table[0, {HeeGenus}];			(** The intersection number of our loop with alpha curves and gamma curves **)
		
		(** The intersection number with alpha curves **)
		For[j=1, j<=HeeGenus, j++,				
			For[ k=Position[ BS[[j]], LoopB[[j, 1]]][[1,1]], True, k=ModUp[k+1, Length[BS[[j]]]],
				If[ LoopB[[j,1]]==LoopB[[j,2]], Break[]];
				s = ES[[ BS[[j, k]], 1 ]];
				r = DS[[ s, 3, ModUp[Position[DS[[s,3]], BS[[j,k]]][[1,1]] +1, Length[DS[[s,3]]] ] ]];		(* r is an alpha edge *)
				If[ MemberQ[ LoopB[[j]], BS[[j,k]] ],
					LA[[ ES[[Abs[r],4]] ]] += Sign[r],
					LA[[ ES[[Abs[r],4]] ]] += 2*Sign[r]
				];
				If[BS[[j, k]]==LoopB[[j,2]], Break[]];
			];
		];
		LA = LA/2;
		
		HL = Join[LC, LA];						(** the homology of our loop **)
		SHL = HL . SFR;							(** the homology in Smith base **)
		For[j=1, j<=2*HeeGenus, j++,
			If[Homology[[j]]==1, SHL[[j]]=0];									(** generator killed by relators **)
			If[Homology[[j]]>1,  SHL[[j]]=Mod[ SHL[[j]], Homology[[j]]]];		(** torsion subgroup **)
		];
		(*Print[GS[[1]], GS[[2]], "  LA=", LA, "  HL=", HL, "    SHL=", SHL];*)

		(** now assign spinc structures **)
		If[ MemberQ[SP, SHL], 	
			FGS[[i, 1]] = Position[SP, SHL][[1,1]],				(** if this spinc structures already exists, find it and assig **)
			SP = Append[SP, SHL];											(** if this spinc structures is not in the list, add it **)
			FGS[[i, 1]] = Length[SP];							(** and now assign it **)
		];
	];
	(*Print[Take[FGS]//MatrixForm, "    Spinc=", SP//MatrixForm];*)
	
	(** assign distance to each, to be used in assigning multiplicities, now we do not care whether alpha or beta. **)
	DiskMult = Table[{-1, n, 0, 0, 0}, {n, 1, Length[DS]}];
	(* distance, disk index, edge adjacent to lower distance disk, lower distance disk, multiplicity *)
	For[DiskMult[[1,1]]=0; HeeDist=0, Count[Transpose[DiskMult][[1]], -1]>0, HeeDist++,
		For[i=1, i<=Enum, i++,
			If[DiskMult[[ES[[i,1]],1]]==HeeDist && DiskMult[[ES[[i,2]],1]]==-1,
				DiskMult[[ES[[i,2]]]]={HeeDist+1, ES[[i,2]], -i, ES[[i,1]], 0};
			];
			If[DiskMult[[ES[[i,2]],1]]==HeeDist && DiskMult[[ES[[i,1]],1]]==-1,
				DiskMult[[ES[[i,1]]]]={HeeDist+1, ES[[i,1]], i, ES[[i,2]], 0};
			];
		];
	];
	DiskSorted = Sort[DiskMult];
	DiskSorted = Transpose[ Take[ Transpose[DiskSorted], 2]];
	MaxDist = Last[DiskSorted][[1]];
	(* Print[DiskMult //MatrixForm, DiskSorted //MatrixForm]; Abort[]; *)

	(** for rational homology spheres, we assign gradings and filtrations**)
	If[RHS, For[spinc=1, spinc<=Length[SP], spinc++,
		(* the index for generators for the Spin^c Structure (spinc) *)
		tmpG = Transpose[Select[FGS, (#[[1]]==spinc)&]][[4]];
		(* Print[tmpG//MatrixForm];*)

		(* assign the grading and filtration of the first generator, in fact we do not bother, but for this to be readable *)
		FGS[[ tmpG[[1]] ]] = {spinc, 0, 0, tmpG[[1]]};

		For[i=2, i<=Length[tmpG], i++,
			(*Print[GS[[ tmpG[[1]]]], GS[[tmpG[[i]]]]];*)
			tmpE = Table[0, {Length[ES]}];		(** the multiplicity of the loop on each edge **)
			tmpD = Table[{0,0}, {Length[DS]}];
			
			LoopB = Table[{0,0}, {HeeGenus}];		(** LoopB **)
			For [j=1, j<=HeeGenus, j++,
				LoopB[[ ES[[GS[[ tmpG[[1]], j, 2]],4]], 1]] = GS[[ tmpG[[1]],j,2]];
				LoopB[[ ES[[GS[[ tmpG[[i]], j, 2]],4]], 2]] = GS[[ tmpG[[i]],j,2]];
			];
			LA = LC = Table[0, {HeeGenus}];			(** The intersection number of our loop with alpha curves and gamma curves **)
			
			(* the intersection number with alpha curves we only need to care about the part on the beta curves*)
			For[j=1, j<=HeeGenus, j++,
				For[ k=Position[ BS[[j]], LoopB[[j, 1]]][[1,1]], True, k=ModUp[k+1, Length[BS[[j]]]],
					If[ LoopB[[j,1]]==LoopB[[j,2]], Break[]];
					s = ES[[ BS[[j, k]], 1 ]];
					r = DS[[ s, 3, ModUp[Position[DS[[s,3]], BS[[j,k]]][[1,1]] +1, Length[DS[[s,3]]] ] ]];		(* r is an alpha edge *)
					If[ MemberQ[ LoopB[[j]], BS[[j,k]] ],
						LA[[ ES[[Abs[r],4]] ]] += Sign[r],
						LA[[ ES[[Abs[r],4]] ]] += 2*Sign[r]
					];
					If[ BS[[j,k]]!=LoopB[[j,1]], tmpE[[ BS[[j,k]] ]] = tmpE[[ BS[[j,k]] ]]+1 ];			(* write the loop in tmpE *)
					If[ BS[[j,k]]==LoopB[[j,2]], Break[]];		
				];
			];
			LA = LA / 2;

			(* write the multiplicity of alpha edges in tmpE *)
			For[j=1, j<=HeeGenus, j++,
				r = GS[[ tmpG[[i]], j, 1]];
				s = GS[[ tmpG[[1]], j, 1]];
				For[ k=Position[ AS[[j]], r ][[1,1]], True, k = ModUp[k+1, Length[AS[[j]]]],
					If[ AS[[j,k]]!=r, tmpE[[ AS[[j,k]] ]] = 1];
					If[ AS[[j,k]]==s, Break[] ];
				];
			];
			(*Print[{Table[i, {i, 1, Length[tmpE]}],tmpE}//MatrixForm];*)

			(* the intersection number with gamma curves *)
			For[j=1, j<=Length[tmpE], j++,
				If[ tmpE[[j]]>0,
					For[k=1, k<=Length[EC[[j]]], k++,
						LC[[ EC[[j,k,1]] ]] = LC[[ EC[[j,k,1]] ]] + EC[[j,k,2]];
					];
				];
			];
			
			(* write the loop in a sum of alpha and beta curves *)
			HL = Join[LC, LA];
			SHL = HL . SFR;
			For[j=1, j<=2*HeeGenus, j++,
				If[ Homology[[j]]>0, SHL[[j]]=SHL[[j]]/Homology[[j]] ];
			];
			SHL = SHL . SFL;			(* HL = SHL . (alpha, beta) *)

			(*Print[{Table[i, {i, 1, Length[tmpE]}],tmpE}//MatrixForm];*)

			(* add alpha and beta curves to make it null homologous *)
			For[j=0, j<=HeeGenus, j++,
				If[SHL[[j]]==0, Continue[]];
				For[k=1, k<=Length[AS[[j]]], k++,
					tmpE[[ AS[[j,k]] ]] = tmpE[[ AS[[j,k]] ]] - SHL[[j]];
				];
			];
			For[j=0, j<=HeeGenus, j++,
				If[SHL[[HeeGenus+j]]==0, Continue[]];
				For[k=1, k<=Length[BS[[j]]], k++,
					tmpE[[ BS[[j,k]] ]] = tmpE[[ BS[[j,k]] ]] - SHL[[HeeGenus+j]];
				];
			];
			(*Print["EdgeMult=", {Table[i, {i, 1, Length[tmpE]}],tmpE}//MatrixForm];*)

			(* now assign a multiplicity to each disk and euler Number*)
			(*Print[DiskMult//MatrixForm, DiskSorted //MatrixForm];*)
			DiskMult = Transpose[ Append[Take[Transpose[DiskMult],4], Table[0, {Length[DS]}]]];
			euler = 0;
			For[j=2, j<=Length[DiskSorted], j++,
				k = DiskSorted[[j,2]];
				DiskMult[[k, 5]] = DiskMult[[ DiskMult[[k,4]], 5]] + Sign[DiskMult[[k,3]]] * tmpE[[ Abs[DiskMult[[k,3]]] ]];
				euler += (4 - 2 * DS[[k,2]]) * DiskMult[[k,5]];
			];
			(*Print["DiskMult=",DiskMult//MatrixForm];*)
			
			(* compute n *)
			For[j=1; n=0, j<=HeeGenus, j++,
				If[GS[[1,j]] == GS[[i,j]], Continue[]];
				r = GS[[1,j,1]];
				n += DiskMult[[ ES[[r,1]],5 ]] + DiskMult[[ ES[[r,2]],5 ]];
				r = AS[[ j,  ModUp[ Position[ AS[[j]], r ][[1,1]]+1, Length[AS[[j]]]] ]];
				n += DiskMult[[ ES[[r,1]],5 ]] + DiskMult[[ ES[[r,2]],5 ]];
				r = GS[[i,j,1]];
				n += DiskMult[[ ES[[r,1]],5 ]] + DiskMult[[ ES[[r,2]],5 ]];
				r = AS[[ j,  ModUp[ Position[ AS[[j]], r ][[1,1]]+1, Length[AS[[j]]]] ]];
				n += DiskMult[[ ES[[r,1]],5 ]] + DiskMult[[ ES[[r,2]],5 ]];
			];
			
			(** MaslovIndex, since the multiplicity of D1 is zero, this is also the grading difference **)
			MaslovIndex = (euler+n)/4;
			(*Print["euler=", euler, "    n=", n, "    MaslovIndex=", MaslovIndex];*)

			(* assign the grading and filtration*)
			FGS[[ tmpG[[i]] ]] = {spinc, MaslovIndex, DiskMult[[2,5]], tmpG[[i]]};
			
			(*Print[ "i=", i, "  ", FGS[[tmpG[[i]]]] ]; *)
		];
	];];
	(*Print["Generators=", Prepend[Transpose[Append[Transpose[FGS], GS]], 
		{"spinc", "grading", "filtration", "index", "vertices"}]//MatrixForm];*)

	(** construct some info about edges to improve the performance **)
	ELeft  = Table[{-1,-1}, {Length[ES]}];		(* first - positive left,  second - negative left  *)
	ERight = Table[{-1,-1}, {Length[ES]}];		(* first - positive right, second - negative right *)
	For[i=1, i<=Length[ES], i++,
		j = ES[[i, 1]];
		k = Position[ DS[[j, 3]], i][[1,1]];
		ELeft [[i,1]] = DS[[j, 3, ModUp[k+1, Length[DS[[j,3]]]] ]];
		ERight[[i,2]] = DS[[j, 3, ModUp[k-1, Length[DS[[j,3]]]] ]];
		j = ES[[i, 2]];
		k = Position[ DS[[j, 3]], -i][[1,1]];
		ELeft [[i,2]] = DS[[j, 3, ModUp[k+1, Length[DS[[j,3]]]] ]];
		ERight[[i,1]] = DS[[j, 3, ModUp[k-1, Length[DS[[j,3]]]] ]];
	];
	(*Print["ELeft=", Transpose[Prepend[Transpose[ELeft ], Range[Length[ES]]]]//MatrixForm, 
		"     ERight=", Transpose[Prepend[Transpose[ERight], Range[Length[ES]]]] //MatrixForm];*)
	EHead = Table[0, {Length[ES]}];
	ETail = Table[0, {Length[ES]}];
	For[i=1, i<=HeeGenus, i++,
		For[len=Length[AS[[i]]]; j=1, j<=len, j++,
			EHead[[AS[[i,j]]]] = AS[[i, ModUp[j+1, len]]];
			ETail[[AS[[i,j]]]] = AS[[i, ModUp[j-1, len]]];
		];
		For[len=Length[BS[[i]]]; j=1, j<=len, j++,
			EHead[[BS[[i,j]]]] = BS[[i, ModUp[j+1, len]]];
			ETail[[BS[[i,j]]]] = BS[[i, ModUp[j-1, len]]];
		];
	];
	(*Print["EHead=", Transpose[{Range[Length[ES]], EHead}]//MatrixForm, 
		"     ETail=", Transpose[{Range[Length[ES]], ETail}] //MatrixForm];*)

	(** first add all elementary disks **)
	IndexZeroDisk = Rest[Transpose[DS][[3]]];
	For[i=1, i<=Length[IndexZeroDisk], i++,
		If[Length[IndexZeroDisk[[i]]]==4,
			If[ Abs[ IndexZeroDisk[[i,1]] ] > Abs[ IndexZeroDisk[[i,3]] ], IndexZeroDisk[[i]] = RotateLeft[ IndexZeroDisk[[i]], 2] ];
		];
	];
	(*Print[IndexZeroDisk //MatrixForm];*)

	(** find all holo disks **)
	For[ len=Length[IndexZeroDisk]; i=1, True, i++,
		If[ i>len, Break[]];
		
		(*Print["IndexZeroDisk[[", i, "]]=", IndexZeroDisk[[i]]];*)

		(** Bigon **)
		If[ Length[IndexZeroDisk[[i]]]==2,
			For[j=1, j<=2, j++,
				edge1 = IndexZeroDisk[[i,j]];
				r = If[ edge1>0, ES[[edge1,2]], ES[[-edge1,1]] ]; (* the disk on the other side *)
				If[ DS[[r,2]]!=2, Continue[]];
				edge2start = edge2 = If[ DS[[r,3,j]]==-edge1, DS[[r,3,2+j]], DS[[r,3,j]] ];
				edgetarget = IndexZeroDisk[[i, ModUp[j+1,2]]];
				While[True,
					edge3 = ERight[[ Abs[edge1], If[edge1>0,1,2] ]];
					(* Print["here", {r, edge1, edge2, edge3, If[ edge1>0, ELeft[[edge1, 1]], ELeft[[-edge1,2]]]}]; *)
					(** if we reach the last edge, we are done, and add the disk **)
					If[ edgetarget == If[ edge1>0, ELeft[[edge1, 1]], ELeft[[-edge1,2]]],
						newdisk = {edge2start, edge3};
						If[Mod[j,2]==0, newdisk=RotateLeft[newdisk,1]];
						(*Print[newdisk];*)
						If[ !MemberQ[IndexZeroDisk, newdisk], IndexZeroDisk = Append[IndexZeroDisk, newdisk]; len++ ];
						Break[];
					];
					r =If[edge3>0, ES[[edge3, 2]], ES[[-edge3,1]]];
					If[ DS[[r,2]]!=2, Break[]];
					If[ edge3>0, 
						edge1 = -ERight[[edge3,  1]]; edge2 = ELeft[[ edge3, 2]],
						edge1 = -ERight[[-edge3, 2]]; edge2 = ELeft[[-edge3, 1]]
					];
				];	(* end While*)
			];
		];
		
		(** Square -- Change the first edge **)
		If[ Length[IndexZeroDisk[[i]]]==4,
			For[j=1, j<=4, j++,
				edge1 = IndexZeroDisk[[i,j]];
				r = If[ edge1>0, ES[[edge1,2]], ES[[-edge1,1]] ];	(* the disk on the other side *)
				If [ DS[[r,2]]!=2, Continue[]];
				edge2start = edge2 = DS[[r, 3, ModUp[ Position[ DS[[r,3]], -edge1][[1,1]]+2, 4] ]];
				edgetarget = IndexZeroDisk[[i, ModUp[j+1,4]]];
				While[True,
					edge3 = If[edge1>0, ERight[[edge1, 1]], ERight[[-edge1, 2]] ];
					(*Print["here", {r, edge1, edge2, edge3, edgetarget,If[ edge1>0, ELeft[[edge1, 1]], ELeft[[-edge1,2]]]}];*)
					If[ edgetarget == If[ edge1>0, ELeft[[edge1, 1]], ELeft[[-edge1, 2]] ],
						newdisk = {edge2start, edge3, IndexZeroDisk[[i, ModUp[j+2, 4]]], IndexZeroDisk[[i, ModUp[j+3,4]]] };
						If[ Mod[j,2]==0, newdisk=RotateLeft[newdisk, 1]];
						If[ Abs[newdisk[[1]]] > Abs[newdisk[[3]]], newdisk= RotateLeft[newdisk, 2]];
						(*Print[newdisk];*)
						If[ ! MemberQ[IndexZeroDisk, newdisk], IndexZeroDisk = Append[IndexZeroDisk, newdisk]; len++ ];
						Break[];
					];
					r =If[edge3>0, ES[[edge3, 2]], ES[[-edge3,1]]];
					If[ DS[[r,2]]!=2, Break[]];
					If[ edge3>0, 
						edge1 = -ERight[[edge3,  1]]; edge2 = ELeft[[ edge3, 2]],
						edge1 = -ERight[[-edge3, 2]]; edge2 = ELeft[[-edge3, 1]]
					];
				];
			];
		];

		(** en, we do not have other type of disks :) **)
	];
	(* Print["IndexZeroDisk=", Transpose[{Range[Length[IndexZeroDisk]], IndexZeroDisk}]//MatrixForm]; *)

	(** find all holomorphic disks, bigons give holomorphic disks, while squares do not necessarily **)
	HoloDisk = Table[{}];
	For[i=1, i<=Length[IndexZeroDisk], i++,
		If[ Length[IndexZeroDisk[[i]]]==4,
			(** if the two opposite alpha edges belong to the same alpha curve, not a differential **)
			If[ ES[[ Abs[IndexZeroDisk[[i,1]]], 4]] == ES[[ Abs[IndexZeroDisk[[i,3]]], 4]], Continue[]];
			(** if the two opposite alpha edges belong to the same alpha curve, not a differential **)
			If[ ES[[ Abs[IndexZeroDisk[[i,2]]], 4]] == ES[[ Abs[IndexZeroDisk[[i,4]]], 4]], Continue[]];
		];
		HoloDisk = Append[HoloDisk, IndexZeroDisk[[i]]];
	];
	(*Print["HoloDisk=", Transpose[{Range[Length[HoloDisk]], HoloDisk}]//MatrixForm];*)

	(** find all differentials from HoloDisk **)	
	For[Diff={}; i=1, i<=Length[HoloDisk], i++,
		If[Length[HoloDisk[[i]]]==2,			(* a bigon *)
			XAE = If[ HoloDisk[[i,1]]<0, -HoloDisk[[i,1]], ETail[[HoloDisk[[i,1]]]] ];
			YBE = If[ HoloDisk[[i,2]]<0, -HoloDisk[[i,2]], ETail[[HoloDisk[[i,2]]]] ];
			YAE = ELeft[[YBE, 1]];
			YAE = If[ YAE<0, -YAE, ETail[[YAE]]];
			AlphaIndex = ES[[XAE, 4]];
			
			XS = EToG[[ XAE ]];
			For[ j=1, j<=Length[XS], j++,
				GXIndex    = XS[[j]];
				GYVertices = GS[[GXIndex]];
				GYVertices[[ AlphaIndex ]] = {YAE, YBE};
				GYIndex = Position[ GS, GYVertices ][[1,1]];
				Diff = Append[ Diff, {GXIndex,GYIndex}];
			];
		];
		If[Length[HoloDisk[[i]]]==4, 
			XAE1 = If[ HoloDisk[[i,1]]<0, -HoloDisk[[i,1]], ETail[[HoloDisk[[i,1]]]] ];
			YBE1 = If[ HoloDisk[[i,2]]<0, -HoloDisk[[i,2]], ETail[[HoloDisk[[i,2]]]] ];
			YAE1 = ELeft[[YBE1, 1]];
			YAE1 = If[ YAE1<0, -YAE1, ETail[[YAE1]]];
			AlphaIndex1 = ES[[XAE1, 4]];

			XAE2 = If[ HoloDisk[[i,3]]<0, -HoloDisk[[i,3]], ETail[[HoloDisk[[i,3]]]] ];
			YBE2 = If[ HoloDisk[[i,4]]<0, -HoloDisk[[i,4]], ETail[[HoloDisk[[i,4]]]] ];
			YAE2 = ELeft[[YBE2, 1]];
			YAE2 = If[ YAE2<0, -YAE2, ETail[[YAE2]]];
			AlphaIndex2 = ES[[XAE2, 4]];
			
			XS = Intersection[ EToG[[XAE1]], EToG[[XAE2]] ];
			(*Print[XAE1, EToG[[XAE1]], "    ", XAE2, EToG[[XAE2]]];*)
			For[ j=1, j<=Length[XS], j++,
				GXIndex    = XS[[j]];
				GYVertices = GS[[GXIndex]];
				GYVertices[[ AlphaIndex1 ]] = {YAE1, YBE1};
				GYVertices[[ AlphaIndex2 ]] = {YAE2, YBE2};
				GYIndex = Position[ GS, GYVertices ][[1,1]];
				Diff = Append[ Diff, {GXIndex,GYIndex}];
			];
		];
	];
	(*Print[Length[Diff], Diff];*)
	
	(** Check whether this is a chain complex **)
	For[i=1, i<=Length[FGS], i++,
		FirstDiff = SecondDiff = Table[{n,0}, {n,1,Length[FGS]}];
		Boundary = Select[Diff, (#[[1]]==i)&];
		For[j=1, j<=Length[Boundary], j++,
			k = Boundary[[j,2]];
			FirstDiff[[k,2]] = Mod[ FirstDiff[[k,2]]+1,2]; 
		];
		FirstDiff = Select[FirstDiff, (#[[2]]==1)&];
		For[j=1, j<=Length[FirstDiff], j++,
			Boundary = Select[Diff, (#[[1]]==FirstDiff[[j,1]])&];
			For[l=1, l<=Length[Boundary], l++,
				k = Boundary[[l,2]];
				SecondDiff[[k,2]] = Mod[ SecondDiff[[k,2]]+1,2]; 
			];
		];
		SecondDiff = Select[SecondDiff, (#[[2]]==1)&];
		If[Length[SecondDiff]>0,
			Print[i, FirstDiff, SecondDiff];
			Break[];
		];
	];
	
	(** check whether there are two differentials from one to another **)
	For[i=1, i<=Length[Diff], i++,
		k = Length[Position[Diff, Diff[[i]]]];
		If[k>1, Print["There are " , k , " ", Diff[[i]]]; Abort[]];
	];
	
	Print["There are ", Length[SP], " Spin^c structures"];
	Print["Use HFhat[spinc_] to display the Heegaard homology, where spinc_ is 1 through ", Length[SP]];
	Print["Use HFKhat[spinc_] to display the Knot Floer homology"];

	CFhat={RHS, FGS, Diff, GS, HoloDisk, Homology, SP};
]

End[];

EndPackage[]