FusionAtlasPaths={"~/projects/fusionatlas/package/","~/Documents/research/fusionatlas/package/","C:\\Documents and Settings\\James Tener\\My Documents\\Research\\fusionatlas\\package\\"};
$Path=$Path~Join~FusionAtlasPaths;
<<FusionAtlas`

EliminatingVinesCode.nb

Calculating N

Utilities

LowestPower[p_,q_]:=-Exponent[p,q,-#&]

HalfIfEven[n_Integer]:=If[Mod[n,2]==0,n/2,n];

FindFj[g_GradedBigraph,q_,j_]:=Module[{x},Expand[(-1)^(GraphRank[g]+j)CharacteristicPolynomial[GraphAdjacencyMatrix[Translate[g,j]],x]/.x->(q+q^(-1))]]

FindModifiedFj[g_GradedBigraph,q_,j_]:=Expand[(q-q^(-1))FindFj[g,q,j]]

Working with A, B,C

CheckACandidate[g_GradedBigraph,candidate_,q_]:=Module[{rank=GraphRank[g]},
And[Expand[q^rank candidate - q^(-rank)(candidate /. q->q^(-1))]==FindModifiedFj[g,q,0],Expand[q^(rank+1) candidate - q^(-rank-1)(candidate /. q->q^(-1))]==FindModifiedFj[g,q,1]]
];

FindAands[g_GradedBigraph,q_, maxsteps_Integer:50] :=
Module[{rank=GraphRank[g],n,oldcandidate=0,candidate},
For[n=rank,n≤rank+maxsteps,n++,
candidate = Expand[q^(-rank-n)( FindModifiedFj[g,q,n]/.q^(k_/;k<0)->0)];
If[And[candidate==oldcandidate,CheckACandidate[g,candidate,q]],Return[{candidate,Max[-LowestPower[candidate,q],0]+1}];];
oldcandidate=candidate;
];
(* Find A, then s is max{-(smallest exponent of A),0} +1 *)
Message["Not found.  Try more steps?"];
Return[$Failed];
];
FindA[g_GradedBigraph,q_, maxsteps_Integer:50] :=FindAands[g,q, maxsteps][[1]]
Finds[g_GradedBigraph, maxsteps_Integer:50] :=Module[{q},FindAands[g,q, maxsteps][[2]]]

(* Returns True if A has only one sign change, and False otherwise.  If True, this means that A is in the "Salem Case" *)
ASalemFormQ[g_GradedBigraph]:=Module[{q},ASalemFormQ[FindA[g,q],q]];
ASalemFormQ[A_,q_]:=SignChanges[A,q]==1;

(* Returns {l,B,C} where A(t)=t^(-l)B(t)C(t), B and C are polynomials, and B is a maximal self-reciprocal factor *)
FactorReciprocal[F_,q_]:=Module[{factorlist=FactorList[F],result,B,laurentC},
{B,laurentC}=(Times @@ Power @@@ # &) /@ GatherBy[factorlist,And[#[[2]]>0,ReciprocalPolynomialQ[#[[1]],q]]&];
{-LowestPower[laurentC,q],B,Expand[q^(-LowestPower[laurentC,q])laurentC]}
];

ReciprocalPolynomialQ[F_,q_]:=Module[{coefficientlist=CoefficientList[F,q],powerlist=Exponent[#,q]&/@ MonomialList[F]},
And[Or[coefficientlist==Reverse[coefficientlist],coefficientlist==-Reverse[coefficientlist]],powerlist==Exponent[F,q]-Reverse[powerlist]]
];

SignChanges[A_,q_]:=Length[Split[Select[CoefficientList[A*q^(Exponent[A,q,-#&]),q],#≠0&],#1*#2>0&]]-1

Calculating R

Calculating R1

BoundR1[A_,q_,maxsteps_Integer:50]:=Module[{B,C,l,m,length,s,kthrootbound,dividends,k,zetak,r,goodk={},period,n,Breciprocalsign,R1BoundList,R1Bound},
DebugPrint["Entering BoundR1 with A = ",A];
{l,B,C}=FactorReciprocal[A,q];
Breciprocalsign=If[CoefficientList[B,q][[1]]/CoefficientList[B,q][[-1]]>0,1,-1];
m=Exponent[C,q];
length=Length[MonomialList[Expand[C]]];
s=Fold[Times,1,Table[Prime[k],{k,1,PrimePi[2*length]}]];
kthrootbound=2*s*m;
dividends=Table[j*s,{j,1,4*m}];
DebugPrint["Upper bound on k for primitive kth root of unity to be a root is ",kthrootbound];

For[k=3,k≤kthrootbound,k++,
If[Mod[k,10000]==0,DebugPrint["At k=",k];];
If[Or@@(Mod[#,k]==0&/@dividends),
zetak=Exp[2Pi I /k];
r=(C/.q->Conjugate[zetak])/(C/.q->zetak);
If[RootReduce[Together[r^(2k)-1==0]],
For[n=1,n≤HalfIfEven[k],n+=1,
If[RootReduce[Together[r-Breciprocalsign * zetak^(2(n-l)+Exponent[B,q])==0]],
DebugPrint["zeta ",k," is a root of F",n];
goodk=Union[goodk,{{k,n}}]];
];];];];

period =HalfIfEven[LCM @@ goodk[[All,1]]];

R1BoundList=1/2*(Total[(Exponent[Cyclotomic[#[[1]],q],q] &)/@ #]&) /@ Table[Select[goodk,Mod[2*#[[2]]-2*p,#[[1]]]==0&],{p,1,period}]+GetRForBRootsOfUnity[B,q];
R1Bound=Max@@R1BoundList;
DebugPrint["|R1| achieves its maximum ",R1Bound," at P",Select[Range[1,Length[R1BoundList]],R1BoundList[[#]]==R1Bound&,1][[1]],"."];
R1Bound
];

Calculating R2 (modulo finding the derivative bound for roots on the circle)

GetRForB[B_,q_]:=Module[{roots=q/.Solve[B==0,q],roottally,rootsofunity,nonrootsofunity,Rforrootsofunity,Rfornonrootsofunity},
If[Not[ListQ[roots]],Return[{0,0}]]; (* B is constant *)
roottally=Tally[roots];
nonrootsofunity=Select[roottally,Not[RootOfUnityQ[#[[1]]]]&];
rootsofunity=Select[roottally,RootOfUnityQ[#[[1]]]&];
(* Count roots of unity besides ±1 as half, to account for complex conjugate *)
Rforrootsofunity=Fold[#1+If[Or[#2[[1]]==1,#2[[1]]==-1],#2[[2]],#2[[2]]/2]&,0,rootsofunity];
Rfornonrootsofunity=Fold[#1+If[#2[[2]]>1,#2[[2]],0]&,0,nonrootsofunity];
{Rforrootsofunity,Rfornonrootsofunity}
];
GetRForBRootsOfUnity[B_,q_]:=GetRForB[B,q][[1]];
GetRForBRepeatedRoots[B_,q_]:=GetRForB[B,q][[2]];

Putting it all together

BoundR[g_GradedBigraph,maxsteps_Integer:50]:=
Module[{A,l,B,C,q,R1Bound,R2BoundCircle,R2BoundLine,RBound,N},
DebugPrint["Entering BoundR."];
A=FindA[g,q,maxsteps];
{l,B,C}=FactorReciprocal[A,q];
R1Bound=BoundR1[A,q,maxsteps];
If[SignChanges[A,q]==1,Return[R1Bound];];
R2BoundCircle=GetRForBRepeatedRoots[B,q];
R2BoundLine=2SignChanges[A,q]+1 ;
RBound=R1Bound+R2BoundCircle+R2BoundLine;
DebugPrint["|R| bound is ",RBound];
RBound
];

FindKGamma[g_GradedBigraph,maxsteps_Integer:50]:=Module[{polyList,q},polyList=#[[1]]^#[[2]]&/@FactorList[FindFj[g,q,Max[0,Finds[g,maxsteps]-GraphRank[g]]]];
Round[Total[DeleteCases[Flatten[(q/.NSolve[#==0,q]&)/@polyList],q]^4]]]

FindN[g_GradedBigraph,maxsteps_Integer:50]:=Module[{q,KGamma,RBound,N,s},
DebugPrint["Entering FindN with g = ",GraphToString[g]];
KGamma=FindKGamma[g,maxsteps];
s=Finds[g,maxsteps];
RBound=BoundR[g,maxsteps];
N = 4*KGamma + 9*RBound;
DebugPrint["K(Γ): ",KGamma,". |R|≤",RBound,". N=",N,". s=",s,"."];
Max[4*KGamma+9*RBound,s]
];

The cyclotomic test

RemoveFactors[polynomial_,factors_List,x_]:=Fold[RemoveFactor[#1,#2,x]&,polynomial,factors];
RemoveFactorOnce[polynomial_,factor_,x_]:=If[Exponent[factor,x]==0,polynomial,
PolynomialQuotientRemainder[polynomial,factor,x]/.{q_,0}:>q/.{_,_}:>polynomial];
RemoveFactor[polynomial_,factor_,x_]:=FixedPoint[RemoveFactorOnce[#,factor,x]&,polynomial];

GetMinimalPolynomialList[g_GradedBigraph,maxVertices_Integer,x_]:=
Module[{charPoly,factors,trialDivisors,guess},
charPoly[k_]:=charPoly[k]=Expand[-x charPoly[k-1]-charPoly[k-2]];
charPoly[t:0|1]:=CharacteristicPolynomial[GraphAdjacencyMatrix[Translate[g,t]],x];
factors[k_]:=FactorList[charPoly[k]/(MinimalPolynomial[DimensionOfGenerator[Translate[g,k]]][x])];
DebugPrint["Entering GetMinimalPolynomialList with g = ",GraphToString[g]," and maxVertices = ",maxVertices,"."];
trialDivisors=DeleteCases[Union[Flatten[Table[First/@factors[k],{k,0,10}]]],_Integer];
Table[
If[Mod[t,20]==0,DebugPrint["Finding minimal polynomial for graph with ",GraphRank[g]+t," vertices."];];
guess=RemoveFactors[charPoly[t],trialDivisors,x];
If[IrreduciblePolynomialQ[guess],
If[And@@(EvenQ/@(Exponent[#,x]&/@MonomialList[guess])),
guess/.x^r_:>x^(r/2),
DebugPrint["Minimal polynomial of d was not even with n = ",GraphRank[g]+t,"."];
MinimalPolynomial[GraphIndex[Translate[g,t]]][x]
],
DebugPrint["Failed to guess the minimal polynomial with n=",GraphRank[g]+t,"."];
MinimalPolynomial[GraphIndex[Translate[g,t]]][x]
],
{t,0,maxVertices-GraphRank[g]}]
]

LeadingTerm[polynomial_,x_]:=
Coefficient[polynomial,x^(Exponent[polynomial,x])]

If[$VersionNumber<6.0,
Discriminant[polynomial_,x_]:=
Together[Resultant[polynomial, D[polynomial,x],x]/LeadingTerm[polynomial,x]]
]

(* Selects from the list of the first n primes the ones that do not divide the discriminant of polynomial *)
GoodPrimes[polynomial_,x_,n_]:=
Module[{Discriminant=Abs[Discriminant[polynomial,x]]}, Select[Table[Prime[k],{k,1,n}],Not[Mod[Discriminant,#]==0]&]];

(* Tests whether or not the given polynomial is the minimal polynomial of a cyclotomic integer.  Can return false positives. *)
PolyCyclotomicTest[polynomial_,primeBound_Integer:200, x_] :=
Module[{j,goodprimelist=GoodPrimes[polynomial,x,primeBound],AllTheSame,FactorPowersModp},
AllTheSame[a_]:=Length[Union[a]]==1;
FactorPowersModp[poly_,p_Integer]:=
Exponent[#[[1]],x]&/@Rest[FactorList[poly,Modulus->p]];

If[Exponent[polynomial,x]==1,Return[True];,
For[j=1,j≤Length[goodprimelist],j++,
If[Not[AllTheSame[FactorPowersModp[polynomial,goodprimelist[[j]]]]],
Return[False];
];];];

(* If the first n primes can't rule it out, we'll return the (potentially false positive) True *)
True
];

(* Runs PolyCyclotomicTest on the list of minimal polynomials coming from the indices of graphs in a vine *)
VineCyclotomicTest[g_GradedBigraph,maxVertices_Integer,primeBound_Integer:200]:=
Module[{graphList,minimalPolyList,goodPolyList,zipList,x},
graphList=Table[Translate[g,t],{t,0,maxVertices-GraphRank[g]}];
minimalPolyList=GetMinimalPolynomialList[g,maxVertices,x];
zipList=Table[{graphList[[t]],minimalPolyList[[t]]},{t,1,Length[graphList]}];
First/@Select[zipList,PolyCyclotomicTest[#[[2]],primeBound,x]&]
];

(* Like PolynomialCyclotomicTest, but returns many primes all of which verify False results *)
MultiPrimePolynomialCyclotomicTest[polynomial_,primeBound_Integer:200, resultPrimeMax_Integer:10,x_] :=
Module[{j,goodprimelist=GoodPrimes[polynomial,x,primeBound],AllTheSame,FactorPowersModp,resultprimes={}},
AllTheSame[a_]:=Length[Union[a]]==1;
FactorPowersModp[poly_,p_Integer]:=
Exponent[#[[1]],x]&/@Rest[FactorList[poly,Modulus->p]];

If[Exponent[polynomial,x]==1,Return[{}];,
For[j=1,j≤Length[goodprimelist],j++,
If[Not[AllTheSame[FactorPowersModp[polynomial,goodprimelist[[j]]]]],
resultprimes=resultprimes~Join~{goodprimelist[[j]]};
If[Length[resultprimes]==resultPrimeMax,Return[resultprimes];];
];];];

(* If the first n primes can't rule it out, we'll return the (potentially false positive) True *)
resultprimes
];

(* Like VineCyclotomicTest, but attempts to find 2 primes that eliminate all of the graphs that fail the test *)
MultiPrimeVineCyclotomicTest[g_GradedBigraph,maxVertices_Integer,primeBound_Integer:200,resultPrimeMax_Integer:10]:=
Module[{graphList,minimalPolyList,goodPolyList,zipList,resultlist,x,FrequentlyOccuring,graphlist,primelist,i,j,totalprimelist},
FrequentlyOccuring[list_]:=Union[Select[Flatten[list],Count[Flatten[list],#]≥Length[list]/2&]];
graphList=Table[Translate[g,t],{t,0,maxVertices-GraphRank[g]}];
minimalPolyList=GetMinimalPolynomialList[g,maxVertices,x];
zipList=Table[{graphList[[t]],minimalPolyList[[t]]},{t,1,Length[graphList]}];
resultlist={#[[1]],MultiPrimePolynomialCyclotomicTest[#[[2]],primeBound,resultPrimeMax,x]}&/@zipList;
{graphlist,totalprimelist,primelist}={First/@Select[resultlist,#[[2]]=={}&],#[[2]]&/@Select[resultlist,#[[2]]≠{}&],FrequentlyOccuring[#[[2]]&/@Select[resultlist,#[[2]]≠{}&]]};

If[Intersection@@totalprimelist≠{},Return[{graphlist,{(Intersection@@totalprimelist)[[1]]}}]];

For[i=1,i≤Length[primelist],i++,
For[j=i+1,j≤Length[primelist],j++,
If[Select[totalprimelist,Intersection[primelist[[i;;i+1]],#]=={}&]=={},Return[{graphlist,primelist[[{i,j}]]}];];
];
];
Return[{graphlist,{}}];
];