(* ::Package:: *)

Clear["Global`*"]
<<Combinatorica`

(* utility functions *)
Nub[list_]:=Fold[If[MemberQ[#1,#2],#1,Append[#1,#2]]&,{},list]
ReducePolyMod[poly_,n_]:=MapIndexed[If[IntegerQ[#],Mod[#,n],#]&,poly,2]
DigitalSum[n_,base_]:=Plus@@Most[DigitCount[n,base]*Range[1,base]]

(* generates the rational cocycle *)
GuaranteedCocycle[deg_,dim_]:=Fold[#1+(-1)^(dim-Length[#2])*(Plus@@#2)^deg&,0,Subsets[Table[Symbol["x"<>ToString[i]],{i,1,dim}]]]//Expand

(* utility functions for generating lots of monomials *)
GenPolyAtoms[deg_,tails_]:=With[{vars=Map[Symbol["x"<>ToString[#]]&,tails]},Map[Fold[Times,1,vars^#]&/@Permutations[Join[#,Table[0,{i,1,Length[tails]-Length[#]}]]]&,Select[Partitions[deg],Length[#]<=Length[tails]&]]]
GenSymPolys[deg_,dim_]:=Plus@@#&/@GenPolyAtoms[deg,Range[1,dim]]
GenMonomials[deg_,dim_]:=Join@@GenPolyAtoms[deg,Range[0,dim-1]]

(* the applicative part of the cocycle condition, encoded as a substitution *)
FirstCocycleRule[n_]:={}
SumCocycleRule[i_,n_]:=Join[Table[Symbol["x"<>ToString[j-1]]->Symbol["x"<>ToString[j-2]],{j,2,i}],{Symbol["x"<>ToString[i]]->Symbol["x"<>ToString[i-1]]+Symbol["x"<>ToString[i]]}]
LastCocycleRule[n_]:=Table[Symbol["x"<>ToString[i]]->Symbol["x"<>ToString[i-1]],{i,1,n}]
CocycleRules[n_]:=Join[{FirstCocycleRule[n]},Table[SumCocycleRule[i,n],{i,1,n}],{LastCocycleRule[n]}]

(* actual additive cocycle condition tools *)
BoundaryMapDescription[deg_,dim_,n_]:=Fold[Plus,0,(GenSymPolys[deg,dim]/.CocycleRules[n])*Table[(-1)^(i+1),{i,1,n+2}]]//Expand
(*BoundaryMapMatrix[deg_,dim_,n_]:=With[{bmd=BoundaryMapDescription[deg,dim,n]},Coefficient[bmd,#]&/@GenMonomials[deg,dim+1]]
FindNullspace[deg_,dim_,n_,char_]:=With[{mat=If[char==0,NullSpace[BoundaryMapMatrix[deg,dim,n]],NullSpace[BoundaryMapMatrix[deg,dim,n],Modulus->char]]},With[{gsp=GenSymPolys[deg,dim]},Plus@@(#*gsp)&/@mat]]*)

(* fast additive cocycle condition tools *)
SparseTerm[term_,gm_,row_]:=If[IntegerQ[First[term]],{{Flatten[Position[gm,term/First[term]]][[1]],row}->First[term]},{{Flatten[Position[gm,term]][[1]],row}->1}]
SparseRow[bmd_,gm_,row_]:=Map[SparseTerm[#,gm,row]&,bmd]
FormSparse[bmd_,gm_,row_]:=If[Length[bmd]==1,SparseRow[bmd[[1]],gm,row],Join[SparseRow[First[bmd],gm,row],FormSparse[Rest[bmd],gm,row+1]]]
BoundaryMapMatrix[deg_,dim_,n_]:=With[{bmd=BoundaryMapDescription[deg,dim,n]},With[{zeroes=Position[bmd,0]},With[{newbmd=Delete[bmd,Position[bmd,0]]},SparseArray[Flatten[Join[FormSparse[newbmd,GenMonomials[deg,dim+1],1]/.Plus->List,Table[{{1,zeroes[[i]][[1]]}->0},{i,1,Length[zeroes]}]]]]]]]
FindNullspace[deg_,dim_,n_,char_]:=With[{mat=If[char==0,NullSpace[BoundaryMapMatrix[deg,dim,n]],NullSpace[BoundaryMapMatrix[deg,dim,n],Modulus->char]]},With[{gsp=GenSymPolys[deg,dim]},Plus@@(#*gsp)&/@mat]]

(* utilities for generating and manipulating polynomials *)
NormalizePoly[poly_]:=If[IntegerQ[First[Factor[poly]]],Expand[Rest[Factor[poly]]],poly]
MakeSymmetric[poly_,dim_]:=Fold[#1+(poly/.#2)&,0,Table[Symbol["x"<>ToString[i]]->Symbol["x"<>ToString[#[[i]]]],{i,1,dim}]&/@Permutations[Range[1,dim]]]//NormalizePoly
MakePoly[exps_]:=MakeSymmetric[Times@@Table[Symbol["x"<>ToString[i]]^exps[[i]],{i,1,Length[exps]}],Length[exps]]

(* JohnMark's set of tools for pretty-printing the nullspace *)
Var[term_]:=If[Head[term]===Power,First[term],term]
Expon[term_]:=With[{copy=term},If[Head[copy]===Power,List[Rest[copy]],{1}]]
TermToExp[term_,vars_]:=If[Length[vars]==1,If[Flatten[List[term]]==={},{0},Expon[term]],If[Head[term]===Times,If[Var[First[term]]===First[vars],Join[Expon[First[term]],TermToExp[Rest[term],Rest[vars]]],Join[{0},TermToExp[term,Rest[vars]]]],If[Var[term]===First[vars],Join[Expon[term],TermToExp[{},Rest[vars]]],Join[{0},TermToExp[term,Rest[vars]]]]]]
TermCoefficient[term_,vars_]:=If[IntegerQ[First[term]],{{First[term]},TermToExp[Rest[term],vars]},{{1},TermToExp[term,vars]}]
MakeExponentListPoly[poly_,vars_]:=With[{list=Sort[Flatten[List[poly/.Plus->List]]],sortvars=Sort[vars]},TermCoefficient[#,sortvars]&/@list]
MakeExponentList[poly_,vars_]:=MapThread[MakeExponentListPoly[#1,#2]&,{poly,vars}]
MakeExponentListSorted[poly_]:=Map[Union,Map[Sort[#,Greater]&,MakeExponentList[poly,Map[Variables,poly]],3]]
FindNullspaceExponents[deg_,dim_,n_,characteristic_]:=With[{cocycles=FindNullspace[deg,dim,n,characteristic]},With[{vars=Map[Variables,cocycles]},Map[Union,Map[Sort[#,Greater]&,MakeExponentList[cocycles,vars],3]]]]

(* routines related to gathering maps *)
CollectMaps[n_]:=Join@@Table[Join[{Symbol["x"<>ToString[i]]->Symbol["x"<>ToString[j]]},Table[Symbol["x"<>ToString[k]]->Symbol["x"<>ToString[k-1]],{k,i+1,n}]],{i,2,n},{j,1,i-1}]
AsymmetricGatherings[poly_,dim_]:=(poly/.#)&/@CollectMaps[dim]
Gatherings[poly_,dim_]:=MakeSymmetric[Plus@@AsymmetricGatherings[poly,dim],dim-1]
GatherMod[poly_,dim_,n_]:=ReducePolyMod[#,n]&/@Gatherings[poly,dim]
ApplyCocycleMap[poly_,n_]:=Plus@@(((poly/.#)&/@CocycleRules[n])*Table[(-1)^i,{i,1,n+2}])

(* delete an instance of an element from a list (returns the original list if it doesn't contain the target) *)
DeleteElement[list_,ele_]:=If[list=={},{},If[First[list]==ele,Rest[list],Join[{First[list]},DeleteElement[Rest[list],ele]]]]
(* delete a list of instance elements from a parent list *)
DeleteFromList[biglist_,targetlist_]:=Fold[DeleteElement[#1,#2]&,biglist,targetlist]

(* number theoretic utilities *)
OverflowQ[a_,b_,p_]:=DigitalSum[a,p]+DigitalSum[b,p]>DigitalSum[a+b,p]
OverflowC[a_,b_,p_]:=CarryCount[{a,b},p]
CarryCount[index_,p_,c_:0]:=If[Divisible[Multinomial@@index,p^(c+1)],CarryCount[index,p,c+1],c]
MinimalCarry[weight_,length_,p_]:=Fold[If[CarryCount[#2,p]<#1,CarryCount[#2,p],#1]&,Infinity,Select[Partitions[weight],Length[#]==length&]]
Deg[poly_]:= Plus@@(Exponent[#,Variables[#]]&/@((poly//Expand)/.Plus->List))[[1]]
ModForm[poly_,char_] :=Fold[(Plus@@(Variables[poly]^#2))*#1&,1,PPartition[Deg[poly],char]]

(* find the constituent monomials in the image of a cocycle and whether or not we run into trouble *)
PPartition[n_,p_]:=If[n==0,{},With[{pull=Log[p,n]//Floor},Join[{p^pull},PPartition[n-p^pull,p]]]]
VisibleSplitsInt[n_,p_]:=With[{IntegerPart=PPartition[n,p]},{Plus@@#,Plus@@DeleteFromList[IntegerPart,#]}&/@Subsets[IntegerPart]//Rest//Most//Nub]
VisibleSplitsIndex[multiindex_,p_,n_]:=Join@@Table[Join[Take[multiindex,i-1],#,Take[multiindex,i-Length[multiindex]]]&/@VisibleSplitsInt[multiindex[[i]],p],{i,1,n}]
VisibleSplitsPoly[multiindex_,p_,n_]:=Join@@(VisibleSplitsIndex[#,p,n]&/@Permutations[multiindex])
CannotHaveDualSplit[index_,p_,n_]:=(Plus@@Table[If[OverflowQ[index[[i]],index[[i+1]],p],0,1],{i,1,n}])<=1
ImageContainsAbsentMonomial[index_,p_,n_]:=Or@@(CannotHaveDualSplit[#,p,n]&/@VisibleSplitsPoly[index,p,n])
DitchCocycles[weight_,length_,p_]:=With[{mincarry=MinimalCarry[weight,length,p]},Select[Select[Partitions[weight],Length[#]==length&],CarryCount[#,p]!=mincarry&]]
DitchSingularities[weight_,length_,p_,n_]:=Select[DitchCocycles[weight,length,p],!ImageContainsAbsentMonomial[#,p,n]&]

(* this will return true if we didn't manage to cull the index with \ref{NoSingularMonomials} *)
TestIndex[index_,p_]:=!Or@@(Or[And[OverflowC[#1,#2,p]==1,DigitalSum[#1+#2,p]>1],OverflowC[#1,#2,p]>1]&@@#&/@Select[Subsets[index],Length[#]==2&])
ReallyDitchSingularities[weight_,length_,p_]:=Select[DitchCocycles[weight,length,p],TestIndex[#,p]&]

(* find the support of a polynomial *)
Regather[image_]:=Nub[Sort[Join[DeleteFromList[image,#],{Plus@@#}],Greater]&/@{{image[[1]],image[[2]]},{image[[2]],image[[3]]}}]
Step[index_,char_]:=Nub[Join@@(Regather/@VisibleSplitsPoly[index,char,2])]
Support[index_,char_]:=FixedPoint[(SortBy[Nub[Join@@(Step[#,char]&/@#)],First]//Reverse)&,{index}]
