(* ::Package:: *)

(* ::Title:: *)
(*A Mathematica Package for Computing McClure' s Obstructions to H\[Infinity] Structure on BP*)


(* ::Text:: *)
(*by Niles Johnson and Justin Noel*)


(* ::Text:: *)
(*This document is a companion to the paper *)
(*"H_{\infty} Orientations on BP"*)
(*by Niles Johnson and Justin Noel*)
(*http://math.uchicago.edu/~justin/H-infty-BP.pdf*)


(* ::Section:: *)
(*Definitions*)
(*(Evaluate the entire section to clear/define all functions)*)


DistributeDefinitions[Context[],"Global`",Power,Times,Plus,arrayPBPvnH,arrayPBPvnJN,
  aseriesNegPowUNXP,formalMcClureSummands,formalMcClureSummandsUNXP,McClureSummandsUNXP,aClass,EulerClass];
SetSharedFunction[aseriesNegPowUNXP,formalMcClureSummands,formalMcClureSummandsUNXP,McClureSummandsUNXP,aClass];
SetSharedVariable[arrayPBPvnH,arrayPBPvnJN,idnty,one,EulerClass];

(* need a higher recursion limit for larger calculations; default is 256 *)
$RecursionLimit = 1000

(* ::Subsection::Closed:: *)
(*initial definitions*)


(*
Leave these for user input

Clear[pVal,order];
pVal=3;
order=15;
*)

Clear[numpowers,nVal];
(* short-hand *)
numpowers=order;
nVal=order+1;

Clear[one,idnty,shiftPower];
(* unit power series *)
one=Join[{1},Table[0,{i,0,order-1}]];
(* identity power series *)
idnty=Join[{0,1},Table[0,{i,0,order-2}]];
(* multiplying by xi^n has the effect of shifting coefficients by n *)
shiftPower[series1_,n_]:=shiftPower[series1,n]=Drop[PadLeft[series1,nVal+n],-n];


Clear[ParallelExpand];
(* apply ExpandAll to a list in parallel *)
(* this is extremely useful, since ExpandAll is a time-consuming process *)
(* it is important, because the rest of our operations are faster when *)
(* ExpandAll has been applied *)
ParallelExpand[Lst_List]:=ParallelTable[ExpandAll[Lst[[i]]],{i,1,Length[Lst]}];


(* ::Subsection::Closed:: *)
(*multiplication and composition of power series as lists*)


Clear[multj,mult,multParallel,multord]

(* Unless specified otherwise, the multiplication and composition functions
compute coefficients up through degree 'order' *)

(* multiply two power series *)
(* multj returns j-th component for mult *)
multj[series1_,series2_,j_]:=Sum[series1[[i+1]]series2[[j-i+1]],{i,0,j}];
(* mult compiles list of multj's *)
mult[series1_,series2_]:=Table[multj[series1,series2,j],{j,0, order}];
(* multParallel does the same as mult, but computes multj in parallel *)
multParallel[series1_,series2_]:=ParallelTable[multj[series1,series2,j],{j,0, order}];
(* multord computes product to smaller order, *) 
(* useful for twovarmult below *)
multord[series1_,series2_,ord_]:=PadRight[Table[multj[series1,series2,j],{j,0, ord}],nVal];

Clear[compi,comp,compseri,compser]
(* compose two power series; input is series and list of powers of another series *)
(* compi returns i-th component for composition *)
compi[series1_,powerlist_,i_]:=Sum[series1[[j+1]]powerlist[[j+1,i+1]],{j,0,order}];
comp[series1_,powerlist_]:=Table[compi[series1,powerlist,i],{i,0, order}];
(* compseri and compser take two series as input, and compute powers
of the second for use in the functions compi and comp *)
compseri[series1_,series2_,i_]:=compseri[series1,series2,i]=compi[series1,pow[series2],i];
compser[series1_,series2_]:=compser[series1,series2]=comp[series1,pow[series2]];

Clear[pow];
(* give list of powers of series *)
(*powers:=powers=FoldList[mult,one,Table[series,{i,1,numpowers}]];*)
pow[series1_]:=pow[series1,numpowers]=FoldList[mult,one,Table[series1,{i,1,numpowers}]];

Clear[tempseries,tc];
(* temp series with coefficients tc[i] for computing general multiplicative or compositional inverses *)
tempseries=Table[tc[i],{i,0,order}];

Clear[multcondi,multcond];
(* conditions for multiplicative inverse *)
(* multcondi returns i-th component for multcond *)
(*multconditionsi[i_]:=multconditionsi[i]=Solve[multj[tempseries,series,i]==one[[i+1]],tc[i]];*)
multcondi[i_,series1_]:=multcondi[i,series1]=Solve[multj[tempseries,series1,i]==one[[i+1]],tc[i]];
(*multconditions:=multconditions=Reverse[Flatten[ParallelTable[multconditionsi[i],{i,0,order}],2]];*)
multcond[series1_]:=multcond[series1]=Reverse[Flatten[Table[multcondi[i,series1],{i,0,order}],2]];

Clear[multinvi,multinv];
(* compute multiplicative inverse *)
(* multinvi returns i-th component for multinv *)
(*multinverseseries:=multinverseseries=Fold[ReplaceAll,tempseries,multconditions];*)
multinvi[i_,series1_]:=multinvi[i,series1]=Fold[ReplaceAll,tempseries[[i+1]],multcond[series1]];
multinv[series1_]:=multinv[series1]=Table[multinvi[i,series1],{i,0,order}];

Clear[compcond,compcond];
(* conditions for compositional inverse *)
(* compcondi returns the i-th component for compcond *)
(*compconditionsi[i_]:=compconditionsi[i]=Solve[compi[tempseries/.tc[0]->0,powers/.b[0]->0,i]==idnty[[i+1]],tc[i]]*)
compcondi[i_,series1_]:=compcondi[i,series1]=Solve[compi[tempseries/.tc[0]->0,pow[series1],i]==idnty[[i+1]],tc[i]]
(*compconditions:=compconditions=Reverse[Flatten[ParallelTable[compconditionsi[i],{i,1,order}],2]];*)
compcond[series1_]:=compcond[series1]=Reverse[Flatten[Table[compcondi[i,series1],{i,1,order}],2]];

Clear[compinvi,compinv];
(* compute compositional inverse *)
(* compinvi returns i-th component for compinv *)
(* compinverseseriesi[i_]:=compinverseseriesi[i]=Fold[ReplaceAll,tempseries[[i+1]],compconditions]/.tc[0]->0; *)
compinvi[i_,series1_]:=compinvi[i,series1]=Fold[ReplaceAll,tempseries[[i+1]],compcond[series1]]/.tc[0]->0;
(*compinverseseries:=compinverseseries=ParallelTable[compinverseseriesi[i],{i,0,order}];*)
compinv[series1_]:=compinv[series1]=Table[compinvi[i,series1],{i,0,order}];


Clear[arrayone];
(* multiplicative identity for two-variable series*)
arrayone=SparseArray[{1,1}->1,{nVal,nVal}];

Clear[twovarser,twovarmultj,twovarmult,twovarmultParallel];
(* two-variable powerseries 
made by adding two single-variable series*)
twovarser[series1_,series2_]:=twovarser[series1,series2]=SparseArray[Table[{i,1}->series1[[i]],{i,1,nVal}],{nVal,nVal}]+SparseArray[Table[{1,i}->series2[[i]],{i,1,nVal}],{nVal,nVal}]
(* j-th row of product *)
twovarmultj[array1_,array2_,j_]:=twovarmultj[array1,array2,j]=Sum[multord[array1[[i+1]],array2[[j-i+1]],order-j],{i,0,j}];
twovarmult[array1_,array2_]:=twovarmult[array1,array2]=Table[twovarmultj[array1,array2,j],{j,0,order}];
twovarmultParallel[array1_,array2_]:=twovarmultParallel[array1,array2]=ParallelTable[twovarmultj[array1,array2,j],{j,0,order}];

Clear[twovarcompi,twovarcomp,twovarcompParallel,twovarcompser,twovarcompserParallel];
Clear[twovarcompXPRi,twovarcompXPR];
(* compose two-variable power series with single-variable power series *)
(* f (g (y,z)) *)
twovarcompi[series1_,twovarpowerlist_,i_]:=twovarcompi[series1,twovarpowerlist,i]=Sum[series1[[j+1]]twovarpowerlist[[j+1,i+1]],{j,0,order}];
twovarcomp[series1_,twovarpowerlist_]:=twovarcomp[series1,twovarpowerlist]=Table[twovarcompi[series1,twovarpowerlist,i],{i,0, order}];
twovarcompParallel[series1_,twovarpowerlist_]:=twovarcompParallel[series1,twovarpowerlist]=ParallelTable[twovarcompi[series1,twovarpowerlist,i],{i,0, order}];
twovarcompser[series1_,array2_]:=twovarcompser[series1,array2]=twovarcomp[series1,twovarpow[array2]];
twovarcompXPRi[series1_,array2_,i_]:=twovarcompXPRi[series1,array2,i]=Sum[series1[[j+1]]twovarpowk[array2,j][[i+1]],{j,0,order}];
twovarcompXPR[series1_,array2_]:=twovarcompXPR[series1,array2]=Table[twovarcompXPRi[series1,array2,i],{i,0,order}];
twovarcompserParallel[series1_,array2_]:=twovarcompserParallel[series1,array2]=twovarcompParallel[series1,twovarpow[array2]];

Clear[twovarpow,twovarpowk,twovarpowParallel];
(* compute powers of two-variable series *)
twovarpow[array1_]:=twovarpow[array1]=FoldList[twovarmult,arrayone,Table[array1,{numpowers}]];
twovarpowk[array1_,0]:=arrayone;
twovarpowk[array1_,1]:=array1;
twovarpowk[array1_,k_]:=twovarpowk[array1,k]=twovarmult[twovarpowk[array1,k-1],array1];
twovarpowParallel[array1_]:=twovarpowParallel[array1]=FoldList[twovarmultParallel,arrayone,Table[array1,{numpowers}]];


(* ::Subsection::Closed:: *)
(*p-typification, Araki, Hazewinkel substitutions*)


Clear[m,v,w,pTypifyMUQ,vnSubstitutionsAraki,vnSubstitutionsHazewinkel,vnSubstitutionsJN,vnSubstitutionsJNName];
m[0]=1;
(*This gives the image of the idempotent map on [CP^n]=(n+1)m_n\inMU_**)
pTypifyMUQ[f_,p_]:=pTypifyMUQ[f,p]=f/.Table[m[i]->If[IntegerQ[Log[p,i+1]],L[Log[p,i+1]],0],{i,0,nVal+1}];
vnSubstitutionsAraki[p_Integer,n_Integer]:=vnSubstitutionsAraki[p,n]=Flatten[Solve[Join[{w[0]==p},Table[w[0] L[m]==Sum[L[i]w[m-i]^(p^i),{i,0,m}],{m,0,Floor[Log[p,n+1]]}]],Table[L[j],{j,1,Floor[Log[p,n+1]]}]]]
vnSubstitutionsHazewinkel[p_Integer,n_Integer]:=vnSubstitutionsHazewinkel[p,n]=Flatten[Solve[Join[{v[0]==p},Table[v[0] L[m]==Sum[L[i]v[m-i]^(p^i),{i,0,Max[m-1,0]}],{m,0,Floor[Log[p,n+1]]}]],Table[L[j],{j,1,Floor[Log[p,n+1]]}]]]
(* function for specifying alternate generators and *)
(* name of substitutions (for choosing filenames) *)
(* 
	when switching to other generators, don't forget to use the notation v[i]
	instead of w[i], since the rest of the formulas expect v's
*)
vnSubstitutionsJN:=vnSubstitutionsJN=vnSubstitutionsHazewinkel[pVal,nVal];
vnSubstitutionsJNName="Hazewinkel";
(*
vnSubstitutionsJN:=vnSubstitutionsJN=vnSubstitutionsAraki[pVal,nVal]/.w[m_]->v[m];
vnSubstitutionsJNName="Araki";
*)


(* drop v_i and w_i for n >= i >= m *)
truncate[f_,m_Integer,n_Integer]:=truncate[f,m,n]=f/.Join[Table[v[i]->If[i<m,v[i],0],{i,0,n}],Table[w[i]->If[i<m,w[i],0],{i,0,n}]]


(* ::Subsection::Closed:: *)
(*log, exp, and formal sum*)


Clear[L,arraylogBP];
L[0]=1;
arraylogBP:=PadRight[SparseArray[
Table[1+pVal^i,{i,0,Floor[Log[pVal,nVal]]}]->Table[L[i],{i,0,Floor[Log[pVal,nVal]]}]
],nVal];

Clear[arrayexpBP];
arrayexpBP:=arrayexpBP=SparseArray[compinv[arraylogBP]];

Clear[arrayiseriesBP,arrayiseriesBPIvar,arraypseriesBPQ,arraypseriesBPvnH,arraypseriesBPvnJN,arraypseriesRedBPvnH,arraypseriesRedBPvnJN,Ivariable,Unk];
arrayiseriesBPIvar:=arrayiseriesBPIvar=compser[arrayexpBP,Ivariable arraylogBP];
arrayiseriesBP[iVal_]:=arrayiseriesBPIvar/.Ivariable->iVal;
arraypseriesBPQ:=arraypseriesBPQ=arrayiseriesBP[pVal];
arraypseriesBPvnH:=arraypseriesBPvnH=ExpandAll[arraypseriesBPQ/.vnSubstitutionsHazewinkel[pVal,nVal]];
arraypseriesBPvnJN:=arraypseriesBPvnJN=ExpandAll[arraypseriesBPQ/.vnSubstitutionsJN];
(* reduced p-series *)
(* last coefficient is now unknown *)
arraypseriesRedBPvnH:=arraypseriesRedBPvnH=PadRight[Drop[arraypseriesBPvnH,1],nVal,Unk];
arraypseriesRedBPvnJN:=arraypseriesRedBPvnJN=PadRight[Drop[arraypseriesBPvnJN,1],nVal,Unk];


Clear[arrayPlusPolyBPi, arrayPlusPolyBPiUNXP, arrayPlusPolyBPiXPR];
(* formal sum of [i]x and y *)
(*arrayPlusPolyBPIvar:=arrayPlusPolyBPIvar=twovarcompser[arrayexpBP,twovarser[arraylogBP,Ivariable arraylogBP]];*)
(*arrayPlusPolyBPIsubi[i_]:=arrayPlusPolyBPIsubi[i]=arrayPlusPolyBPIvar/.Ivariable->i;*)
arrayPlusPolyBPi[0]:=SparseArray[{2,1}->1,{nVal,nVal}];
arrayPlusPolyBPi[i_]:=arrayPlusPolyBPi[i]=SparseArray[ParallelExpand[twovarcompser[arrayexpBP,twovarser[arraylogBP,i arraylogBP]]]];
arrayPlusPolyBPiUNXP[i_]:=arrayPlusPolyBPiUNXP[i]=twovarcompXPR[arrayexpBP,twovarser[arraylogBP,i arraylogBP]];
arrayPlusPolyBPiXPR[i_]:=arrayPlusPolyBPiXPR[i]=twovarcompXPR[arrayexpBP,twovarser[arraylogBP,i arraylogBP]];


(* ::Subsection::Closed:: *)
(*power operation, Euler class, a_i*)


Clear[arrayPBPQ,arrayPBPQUNXP,arrayPBPQformal,arrayPBPQsymmk,arrayPBPQsymm,arrayPBPvnH,arrayPBPvnHsymm,arrayPBPvnJN,arrayPBPvnJNUNXP,arrayPBPvnJNsymm];
(* Power operation P_p on BP *)
(* note: *)
(* McClure uses \prod_1^{p-1}; we use \prod_0^{p-1}*)
(* but: we implement this in such a way as to not loose accuracy in the last term *)
(* this is an array with nVal columns and nVal+1=order+2 rows *)
arrayPBPQ:=arrayPBPQ=Prepend[Fold[twovarmult,arrayone,Table[arrayPlusPolyBPi[i],{i,1,pVal-1}]],Table[0,{i,nVal}]];
arrayPBPQUNXP[n_]:=arrayPBPQUNXP[n]=Prepend[Fold[twovarmult,arrayone,Table[arrayPlusPolyBPiUNXP[i],{i,1,pVal-1}]],Table[0,{i,nVal}]][[n]];
arrayPBPQformal:=arrayPBPQformal=Prepend[Fold[twovarmult,arrayone,Table[arrayPlusPolyBPi[n[i]],{i,1,pVal-1}]],Table[0,{i,nVal}]];
(* here is our niave implementation; the term i = 0 causes a loss of accuracy *)
(*
arrayPBPQ:=arrayPBPQ=Fold[twovarmult,arrayone,Table[arrayPlusPolyBPi[i],{i,0,pVal-1}]];
arrayPBPQformal:=arrayPBPQformal=Fold[twovarmult,arrayone,Table[arrayPlusPolyBPi[n[i]],{i,0,pVal-1}]];
*)
(* compute symmetric reduction of k-th term and then use ParallelTable to compute all terms *)
(* warning: we use only component [[1]] of SymmetricReduction, *)
(*    if there were a remainder, we would be dropping it *)
arrayPBPQsymmk[k_]:=arrayPBPQsymmk[k]=Table[SymmetricReduction[arrayPBPQformal[[k,j]],Table[n[i],{i,0,pVal-1}],Table[S[i],{i,1,pVal}]][[1]],{j,1,nVal}];
arrayPBPQsymm:=arrayPBPQsymm=ParallelTable[arrayPBPQsymmk[k],{k,1,nVal}];
arrayPBPvnH:=arrayPBPvnH=ParallelExpand[arrayPBPQ/.vnSubstitutionsJN];
arrayPBPvnHsymm:=arrayPBPvnHsymm=ParallelExpand[arrayPBPQsymm/.vnSubstitutionsJN];
arrayPBPvnJN:=arrayPBPvnJN=ParallelExpand[arrayPBPQ/.vnSubstitutionsJN];
arrayPBPvnJNUNXP:=arrayPBPvnJNUNXP=arrayPBPQUNXP/.vnSubstitutionsJN;
arrayPBPvnJNsymm:=arrayPBPvnJNsymm=ParallelExpand[arrayPBPQsymm/.vnSubstitutionsJN];



Clear[arrayPBPvnJNRedi,arrayPBPvnJNRed];
arrayPBPvnJNRedi[i_]:=arrayPBPvnJNRedi[i]=arrayCoeffElimLast[arrayPBPvnJN[[i]]];
(*arrayPBPvnJNRedi[i_]:=arrayPBPvnJNRedi[i]=arrayCoeffElimAllLast[arrayPBPvnJN[[i]],order];*)
(* entries outside our range of accuracy are set to 'Unk' *)
arrayPBPvnJNRed:=arrayPBPvnJNRed=SparseArray[ReplacePart[ParallelTable[arrayPBPvnJNRedi[i],{i,order+2}],{a_,b_}/;b>nVal-a+2 -> 0]];



Clear[EulerClass,aClass,EulerClassRed,aClassRed,seriesList];
Clear[aClassXPRed];
(* Extract rows from arrayPBPvnJN, and replace head with 'seriesList' *)
(* Note that aClass[r] is undefined for r greater than 'order' *)
EulerClass:=EulerClass=ReplacePart[arrayPBPvnJN[[2]],0->seriesList];
aClass[n_]:=aClass[n]=ReplacePart[arrayPBPvnJN[[n+2]],0->seriesList];
aClassXPRed[n_]:=aClassXPRed[n]=seriesList@@ReplacePart[arrayCoeffElimLast[ExpandAll[arrayPBPQUNXP[n+2]/.vnSubstitutionsJN]],b_/;b>nVal-(n+2)+2 -> 0];

(* try UNXP versions...*)
EulerClassRed:=aClassXPRed[0];
aClassRed[n_]:=aClassXPRed[n];
(*
EulerClassRed:=EulerClassRed=ReplacePart[arrayPBPvnJNRed[[2]],0->seriesList];
aClassRed[n_]:=aClassRed[n]=ReplacePart[arrayPBPvnJNRed[[n+2]],0->seriesList];
*)



(* ::Subsection::Closed:: *)
(*dividing by p-series*)


Clear[FactorSummand,DivisionAlg];
Clear[arrayFirstUnreducedTerm,arrayCoeffElim,arrayCoeffElimFirst];
  (* functions for reducing one power series mod another *)
(* two improvements: 
   1. handle coefficients in Z_ (p) 
   2. handle terms including sums of v[i]-multinomials
*)

Clear[reductionFactor];
(* return factor by which <p> should be multiplied to reduce X modulo <p> *)

(* if input is list, determine list position of first non-zero factor *)
reductionFactor[X_List]:=reductionFactor[X]=
    Catch[Do[
	If[
	    reductionFactor[X[[i]]] =!= 0,
	    Throw[{i-1,reductionFactor[X[[i]]]}]
	  ],
	{i,1,Length[X]}
	    ];
	  Throw[False]];
(* if input is sum, determine first non-zero factor for summand *)
reductionFactor[X_Plus]:=reductionFactor[X]=
    If[
	reductionFactor[X[[1]]] =!= 0, 
	reductionFactor[X[[1]]],
	reductionFactor[Drop[X,1]]
      ];
(* if input is neither of the above, just determine factor *)
reductionFactor[X_]:=reductionFactor[X]=
    Quotient[Numerator[FactorTermsList[X][[1]]
		      ],pVal]*
    FactorTermsList[X][[2]]/
    Denominator[FactorTermsList[X][[1]]];



(* find first unreduced term *)


arrayFirstUnreducedTerm[L_List]:=Catch[Do[If[FactorTermsList[L[[i]]][[1]]<0||FactorTermsList[L[[i]]][[1]]>pVal-1,Throw[i-1]],{i,1,Length[L]}];Throw[False]];

(* factor summand number s *)
(*
FactorSummand[f_,s_]:=FactorSummand[f,s]=FactorTermsList[MonList[f,x][[s]]];
DivisionAlg[M_,N_]:=DivisionAlg[M,N]={Quotient[M,N],Mod[M,N]};
*)

Clear[arrayCoeffElimNew,arrayCoeffElimShowAll,arrayCoeffElimLast];
arrayCoeffElimNew[L_List]:=arrayCoeffElimNew[L]=
    If[reductionFactor[L] =!= False,
       ExpandAll[L-reductionFactor[L][[2]] shiftPower[arraypseriesRedBPvnJN,reductionFactor[L][[1]]]],
       L
      ];
		
						     

arrayCoeffElim[list1_,s_Integer]:=arrayCoeffElim[list1]=ExpandAll[list1-Quotient[FactorTermsList[list1[[s+1]]][[1]],pVal] FactorTermsList[list1[[s+1]]][[2]] shiftPower[arraypseriesRedBPvnJN,s]];
arrayCoeffElim[list1_,False]:=list1;
arrayCoeffElimFirst[list1_]:=arrayCoeffElimFirst[list1]=arrayCoeffElim[list1,arrayFirstUnreducedTerm[list1]];

Clear[vnMax];
(*Clear[arrayCoeffElimAll,arrayCoeffElimAllLast];*)
(*  now we use 'arrayCoeffElimLast' -- the new and improved reducer! *)

  (* returns list of successive reductions of f by pPolyRed *)
  (*   these are pairs {f_i,d_i} *)
  (*   where f_i = f_ {i+1} + {p}(x)*d_ {i+i} *)
(*
arrayCoeffElimAll[list1_,iMax_]:=Module[{h,i},h=list1; Print["Input"]; Print[list1]; i=1; Print["Reductions"]; While[arrayFirstUnreducedTerm[h]=!=False && i<iMax+1,h=arrayCoeffElimFirst[h];Print[h];i++]];
*)
(*arrayCoeffElimAllLast[list1_,iMax_]:=Module[{h,i},h=list1; i=1; While[arrayFirstUnreducedTerm[h]=!=False && i<iMax+1,h=arrayCoeffElimFirst[h];i++];h];*)

(* determine maximum power to which v_i can occur in deg. less than or equal to 2r *)
vnMax[r_]:=vnMax[r]=Table[Floor[r/(pVal^i - 1)],{i,Floor[Log[pVal,r+1]]}];

arrayCoeffElimShowAll[L_List,iMax_]:=
    Module[{h,i},h=L; Print["Input"]; Print[L//InputForm]; i=1; Print["Reductions"];
	   While[reductionFactor[h]=!=False && i < iMax+1,
		 h=arrayCoeffElimNew[h];Print[h//InputForm];i++
		]];

(* use rough upper bound for number of times we should have to reduce *)
(* just to be sure we do eventually stop *)
arrayCoeffElimLast[L_List]:=
    FixedPoint[arrayCoeffElimNew,L,Times@@Map[# + 1 &,vnMax[order]]];


(* ::Subsection::Closed:: *)
(*extend definition of Power and Times for expressions with head "seriesList"*)


(* Change head of output back to "seriesList" *)

(* Use head "Eu" for Euler class; Eu[n] is Eu^n *)
Unprotect[Power];
Power[x_seriesList,n_Integer]:=ReplacePart[pow[x][[n+1]],0->seriesList];
Power[x_Eu,n_Integer]:=Eu[x[[1]] n];
Protect[Power];

Unprotect[Times];
Times[x_seriesList,y_seriesList]:=ReplacePart[mult[x,y],0->seriesList];
Times[x_Integer,y_seriesList]:=ReplacePart[x ReplacePart[y,0->List],0->seriesList];
Times[x_Rational,y_seriesList]:=ReplacePart[x ReplacePart[y,0->List],0->seriesList];
Times[x_L,y_seriesList]:=ReplacePart[x ReplacePart[y,0->List],0->seriesList];
Times[x_m,y_seriesList]:=ReplacePart[x ReplacePart[y,0->List],0->seriesList];
Times[x_v,y_seriesList]:=ReplacePart[x ReplacePart[y,0->List],0->seriesList];
Times[x_v^n_Integer,y_seriesList]:=ReplacePart[x^n ReplacePart[y,0->List],0->seriesList];
Times[x_Eu,y_Eu]:=Eu[x[[1]]+y[[1]]];
Protect[Times];

Unprotect[Plus];
Plus[x_seriesList,y_seriesList]:=ReplacePart[ReplacePart[x,0->List]+ReplacePart[y,0->List],0->seriesList];
Protect[Plus];


(* ::Subsection::Closed:: *)
(*McClure[n]*)


Clear[seriesf,f,fZeroSubstitutions,genmultinv,Eu,aseriesNegPowUNXP,formalMcClureSummands,formalMcClureSummandsUNXP];
(* series with formal coefficients f[i] for computing multiplicative inverse *)
seriesf=Table[f[i],{i,0,order}];
(* detect which of the a_i are zero, and create a rule which sends the corresponding coefficients of seriesf to zero *)
fZeroSubstitutions:=fZeroSubstitutions=Inner[Rule,Table[f[i],{i,1,order}],Table[If[aClassRed[i]===seriesList@@Table[0,{nVal}],0,f[i]],{i,1,order}],List];

Eu[0]=1;
genmultinv:=genmultinv=multinv[seriesf/.f[0]->Eu[1]/.fZeroSubstitutions];
(* aseries^-(n) *)
aseriesNegPowUNXP[0]=1;
aseriesNegPowUNXP[1]:=genmultinv;
aseriesNegPowUNXP[num_]:=aseriesNegPowUNXP[num]=mult[aseriesNegPowUNXP[num-1],genmultinv];
(* no need to expand those parts which we'll multiply by 0 anyway *)

formalMcClureSummandsUNXP[n_]:=formalMcClureSummandsUNXP[n]=Eu[2n+1]Inner[Times,PadRight[SparseArray[pTypifyMUQ[Table[m[n-k]*(n-k+1),{k,0,n}],pVal]],nVal],aseriesNegPowUNXP[n+1],List]/.vnSubstitutionsJN;
formalMcClureSummands[n_]:=formalMcClureSummands[n]=ParallelExpand[formalMcClureSummandsUNXP[n]];

Clear[McClureSummandsXP,McClureSummandsUNXP,McClure];
McClureSummandsUNXP[n_]:=McClureSummandsUNXP[n]=formalMcClureSummands[n]/.f[i_]:>aClassRed[i]/.Eu[j_]:>EulerClassRed^j;
McClureSummandsXP[n_]:=McClureSummandsXP[n]=ParallelExpand[McClureSummandsUNXP[n]];
McClure[n_]:=McClure[n]=Apply[Plus,
	If[Head[#]===seriesList,Apply[List,#],PadRight[{#},nVal]]& /@ McClureSummandsXP[n]
];
RedMcClure[n_]:=arrayCoeffElimLast[McClure[n]];

(* simplified check for first obstruction, MC_ {2(p-1)]/((2p-1) * EulerClass^{2p-4}) *)
Clear[FirstObstruction,FirstObstructionRed];
FirstObstruction:=FirstObstruction=ExpandAll[List@@(
    pVal*(L[1]/.vnSubstitutionsJN)*(-1)*EulerClassRed*aClassRed[pVal-1]+
    (-1)*EulerClassRed*aClassRed[2*(pVal-1)]+
    pVal*aClassRed[pVal-1]^2)]
    
FirstObstructionRed:=arrayCoeffElimLast[FirstObstruction]


(* ::Subsection:: *)
(*output formatting*)


Clear[PrettyRedpSeries,PrettyMcClure,PrettyRedMcClure];

PrettyRedpSeries:=Dot[arraypseriesRedBPvnJN,Table[\[Xi]^i,{i,0,order}]]/.v[j_]->Subscript[v, j];
PrettyMcClure[n_]:=Dot[McClure[n],Table[\[Xi]^i,{i,0,order}]]/.v[j_]->Subscript[v, j];
PrettyRedMcClure[n_]:=Dot[arrayCoeffElimLast[McClure[n]],Table[\[Xi]^i,{i,0,order}]]/.v[j_]->Subscript[v, j];


Clear[niceOut];
niceOut[A_,x_]:=Dot[Table[x[i],{i,0,order}],A]/.v[j_]->Subscript[v, j]/.L[k_]->Subscript[\[ScriptL],k]/.x[i_]->x[t]^i;


NightRun[n_]:=Timing[{
DateString[],
{"n value",n},
{"McClure[n]",McClure[n]},
{"p-series",arraypseriesBPvnJN},
{"remainder",arrayCoeffElimLast[McClure[n]]}
}];
NightRunFO:=Timing[{
DateString[],
{"Computing first possible obstruction, MC_{2(p-1)}"},
{"p-series",arraypseriesBPvnJN},
{"EulerClassRed",EulerClassRed},
{"aClassRed[p-1]",aClassRed[pVal-1]},
{"aClassRed[2(p-1)]",aClassRed[2(pVal-1)]},
{"FirstObstructionRed",FirstObstructionRed}
}];
PBPOutFile:=outLoc[math]<>vnSubstitutionsJNName<>"vn.arrayPBPvnJN."<>IntegerString[pVal]<>"."<>IntegerString[order]<>".txt"
PBPOutFileOrochi:=outLoc[orochi]<>vnSubstitutionsJNName<>"vn.arrayPBPvnJN."<>IntegerString[pVal]<>"."<>IntegerString[order]<>".txt"
OutFile[n_]:=outLoc[math]<>vnSubstitutionsJNName<>"vn.output."<>IntegerString[pVal]<>"."<>IntegerString[n]<>"."<>IntegerString[order]<>".txt"
OutFileOrochi[n_]:=outLoc[orochi]<>vnSubstitutionsJNName<>"vn.output."<>IntegerString[pVal]<>"."<>IntegerString[n]<>"."<>IntegerString[order]<>".txt"