(* ::Package:: *)

BeginPackage["SqueezeRoutines`"]

(*

Eddy Ardonne and Nicolas Regnault

Squeeze Routines for polarized and spinful quantum Hall states,
bosonic as well as fermionic. Routines for the spinful case
are based on

The structure of spinful quantum Hall states: a squeezing perspective,
E. Ardonne, N. Regnault,
Phys. Rev. B 84, 205134 (2011);
arXiv:1107.2232

See this paper for references for the polarized case, and for reference
about particle and orbital entanglement spectra.

Routines making use of the additional symmetries in case of L=0 and/or
S=0 states are not implemented in this version. A routine to calculate
only the number of solutions to the consistency condition via an
eigenvalue problem is not included either. Please contact E. Ardonne
about these (and possible other) issues.

*)

bosonsqueeze::usage=
 "bonsonsqueeze[rp_] generetes the squeezed Hilbert space, from the
  root partition rp, assuming that the particles are spinless bosons."

fermionsqueeze::usage=
 "fermionsqueeze[rc_] generetes the squeezed Hilbert space, from the
  root configuration rc, assuming that the particles are spinless fermions.
  No checks are done to see if the root configuration qualifies!"

spinbosonsqueeze::usage=
 "spinbosonsqueeze[rc_,nup,ndown] generetes the squeezed Hilbert space, from the
  root configuration rc, assuming that the particles are spinful bosons, with
  nup spin-up particles, and ndown spin-down particles."

spinfermionsqueeze::usage=
 "spinfermionsqueeze[rc_,nup,ndown] generetes the squeezed Hilbert space, from the
  root configuration rc, assuming that the particles are spinful fermions, with
  nup spin-up particles, and ndown spin-down particles.
  No checks are done to see if the root configuration qualifies!"

bosonlplus::usage=
 "bosonlplus generates the L+ conditions for polarized boson states."

fermionlplus::usage=
 "fermionlplus generates the L+ conditions for polarized fermion states."

spinbosonlplus::usage=
 "spinbosonlplus generates the L+ conditions for spinful boson states."

spinfermionlplus::usage=
 "spinfermionlplus generates the L+ conditions for spinful fermion states."

spinbosonsplus::usage=
 "spinbosonsplus generates the S+ conditions for spinful boson states."

spinfermionsplus::usage=
 "spinfermionsplus generates the S+ conditions for spinful fermion states."

showstatenum::usage=
 "showstatenum[i] gives the state with label `i' in the number representation."

showstateorb::usage=
 "showstateorb[i] gives the state with label `i' in the orbital representation."

solveall::usage=
 "Solve the L+ (and S+ if applicable) conditions exactly. It is assumed that
  no additional conditions are present. The solutions are storedin b[i], and the
  number of solutions is returned."
 
solveallcond::usage=
 "Solve the L+ (and S+ if applicable) conditions, in addition to the extra conditions
  specified via `condlist' exactly. The solutions are stored in b[i], and the number
  of solutions is returned."

solveallnum::usage=
 "Solve the L+ (and S+ if applicable) conditions exactly. It is assumed that
  no additional conditions are present. The solutions are stored in b[i], and the
  number of solutions is returned."
 
solveallcondnum::usage=
 "Solve the L+ (and S+ if applicable) conditions, in addition to the extra conditions
  specified via `condlist' exactly. The solutions are stored in b[i], and the number
  of solutions is returned."

bosonl2::usage=
 "Generates the l^2 operator for the current hilbertspace, in the case of spinless bosons.
  Returned as a sparst matrix."

bosonpol::usage=
 "Generate the explicit symmetric polynomial, based on the stored basis, and the
  coefficients b[i], in terms of the variables z[i]."

fermionpol::usage=
 "Generate the explicit anti-symmetric polynomial, based on the stored basis, and the
  coefficients b[i], in terms of the variables z[i]."

spinbosonpol::usage=
 "Generate the explicit polynomial symmetric in the w[i] (spin up coordinates)
  and symmetric in the z[i] (spin down coordinates), based on the stored basis,
  and the coefficients b[i], in terms of the variables z[i]."

spinfermionpol::usage=
 "Generate the explicit polynomial anti-symmetric in the w[i] (spin up coordinates)
  and anti-symmetric in the z[i] (spin down coordinates), based on the stored basis,
  and the coefficients b[i], in terms of the variables z[i]."

partition::usage=
 "partition[orbocc] generates the partition corresponding to a set of
  orbital occupation numbers. Used by other packages."

orbocc::usage=
 "orbocc[partition,flux] generates the (flux+1) orbital occupation numbers
  corresponding to a partition. Used by other packages."

integerdigits::usage=
 "integerdigits[num,base,norb] gives the norb orbital occupation numbers, given num,
  which is the number representing the state, and the base used. Used by other packages."

spinintegerdigits::usage=
 "spinintegerdigits[nums,base,norb] gives the norb orbital occupation numbers, given nums,
  which are the numbers representing the spinful state, and the base used. Used by other packages."

fromdigits::usage=
 "fromdigits[rp,base] gives the number represeting the state specified by the orbital occupation
  numbers rp, with the specified base. Used by other packages."

spinfromdigits::usage=
 "spinfromdigits[rp,base] gives the numbers represeting the spinful state specified by the
  orbital occupation numbers rp, with the specified base. Used by other packages."

spindist::usage=
 "spindist[nu,nd] gives the possible `spin-distributions' with nu up particles, and
  nd down particles. Used by other packages."


orbtopart::usage=
 "orbtopart[part] converts orbital occupation numbers to the length of the rows of the
  corresponding partition. Used by other packages."

rootconfigs::usage=
 "rootconfigs[Ne,nflux,cellsize,maxincell] gives the root configurations, with
  non-negative lz values, which can be used to generate certain states with the
  squeezing rules."


Begin["`Private`"]

(* Some functions, needed internally *)

partition[orbocc_]:=
Catch[
 For[
  v={};n=Plus@@orbocc;j=Length[orbocc],
  j>0,
  j--,
  v=Flatten[Append[v,Table[j-1,{jj,1,orbocc[[j]]}]]];If[j==1,Throw[v]]
 ]
]

orbocc[part_,flux_]:=Table[Count[part,x_/;x==i],{i,0,flux}]

orbtopart[part_]:=
 Catch[
  For[
   i=1;res={},
   True,
   i++,
   Do[res=Append[res,i],{j,part[[i]]}];If[i==Length[part],Throw[res]]
  ]
 ]

integerdigits[num_,base_,norb_]:=
 PadLeft[IntegerDigits[num,base],norb]

spinintegerdigits[nums_,base_,norb_]:=
 {PadLeft[IntegerDigits[nums[[1]],base],norb],
  PadLeft[IntegerDigits[nums[[2]],base],norb]}

fromdigits[rp_,base_]:=
 FromDigits[rp,base]

spinfromdigits[drp_,base_]:=
 {FromDigits[drp[[1]],base],
  FromDigits[drp[[2]],base]}

spindist[nup_,ndown_]:=
 Permutations[Table[If[j<=nup,1,2],{j,1,nup+ndown}]]

showstatenum[i_]:=Global`state[i];
showstateorb[i_]:=
 With[{st=Global`state[i]},
  If[NumberQ[st],
   integerdigits[st,Global`base,Global`norbitals],
   spinintegerdigits[st,Global`base,Global`norbitals]]
]


(** Squeeze routines **)

(* spinless bosons *)

bosonsqueeze[rp_]:= With[{},
 norbitals=Length[rp];
 nparticles=Plus@@rp;
 base=nparticles+1;
 Clear[Global`state,Global`stateposition,statepresent,Global`b];
 hsdim=1;
 olddim=1;
 currentdim=1;
 Global`state[hsdim]=fromdigits[rp,base];
 Global`stateposition[Global`state[hsdim]]=hsdim;
 statepresent[Global`state[hsdim]]=True;
 new=True;
 While[new,
  new=False;
  Do[
   currentpart=integerdigits[Global`state[i],base,norbitals];
   left=Flatten[Position[currentpart,x_/;x!=0]];
   right=Flatten[Table[Cases[left,x_/;x>left[[i1]]+1,1,1],{i1,1,Length[left]-1}]];
   Do[
    newpart=If[right[[j]]==left[[j]]+2,
     ReplacePart[currentpart,{left[[j]]->currentpart[[left[[j]]]]-1,right[[j]]->currentpart[[right[[j]]]]-1,left[[j]]+1->currentpart[[left[[j]]+1]]+2}],
     ReplacePart[currentpart,{left[[j]]->currentpart[[left[[j]]]]-1,right[[j]]->currentpart[[right[[j]]]]-1,left[[j]]+1->currentpart[[left[[j]]+1]]+1,right[[j]]-1->currentpart[[right[[j]]-1]]+1}]
     ];
    newstate=FromDigits[newpart,base];
    If[statepresent[newstate],Null,Null,
    new=True;
    hsdim++;
    Global`state[hsdim]=newstate;
    Global`stateposition[newstate]=hsdim;
    statepresent[newstate]=True;
    ];
   ,{j,1,Length[right]}];
  ,{i,olddim,currentdim}];
 olddim=currentdim+1;currentdim=hsdim;
 ];
 Global`norbitals=norbitals;
 Global`nparticles=nparticles;
 Global`base=base;
 Global`hsdim=hsdim;
]

(* spinless fermions *)

fermionsqueeze[rp_]:= With[{},
 norbitals=Length[rp];
 nparticles=Plus@@rp;
 base=nparticles+1;
 Clear[Global`state,Global`stateposition,statepresent,Global`b];
 hsdim=1;
 olddim=1;
 currentdim=1;
 Global`state[hsdim]=fromdigits[rp,base];
 Global`stateposition[Global`state[hsdim]]=hsdim;
 statepresent[Global`state[hsdim]]=True;
 new=True;
 While[new,
  new=False;
  Do[
   currentpart=integerdigits[Global`state[i],base,norbitals];
   left=Cases[Flatten[Position[currentpart,1]],x_/;x<norbitals-2&&currentpart[[x+1]]!=1];
   preright=Cases[Flatten[Position[currentpart,1]],x_/;x>3&&currentpart[[x-1]]!=1];
   right=Flatten[Table[Cases[preright,x_/;x>left[[i]]+2,1,1],{i,1,Length[left]}]];
   Do[
    If[currentpart[[left[[j]]+1]]==currentpart[[right[[j]]-1]]==0,
     newpart=ReplacePart[currentpart,{left[[j]]->currentpart[[left[[j]]]]-1,right[[j]]->currentpart[[right[[j]]]]-1,left[[j]]+1->currentpart[[left[[j]]+1]]+1,right[[j]]-1->currentpart[[right[[j]]-1]]+1}];
     newstate=FromDigits[newpart,base];
     If[statepresent[newstate],Null,Null,
      new=True;
      hsdim++;
      Global`state[hsdim]=newstate;
      Global`stateposition[newstate]=hsdim;
      statepresent[newstate]=True;
    ]
   ]
   ,{j,1,Length[right]}];
  ,{i,olddim,currentdim}];
 olddim=currentdim+1;currentdim=hsdim;
 ];
 Global`norbitals=norbitals;
 Global`nparticles=nparticles;
 Global`base=base;
 Global`hsdim=hsdim;
]


(* spinful bosons *)

spinbosonsqueezeint[rp_]:=With[{},
 norbitals=Length[rp];
 nparticles=Plus@@rp;
 base=nparticles+1;
 Clear[intstate,intstateposition,intstatepresent];
 inthsdim=1;
 olddim=1;
 currentdim=1;
 intstate[inthsdim]=FromDigits[rp,base];
 intstateposition[intstate[inthsdim]]=inthsdim;
 intstatepresent[intstate[inthsdim]]=True;
 new=True;
 While[new,
  new=False;
  Do[
   currentpart=integerdigits[intstate[i],base,norbitals];
   left=Flatten[Position[currentpart,x_/;x!=0]];
   right=Flatten[Table[Cases[left,x_/;x>left[[i1]]+1,1,1],{i1,1,Length[left]-1}]];   
   Do[
    newpart=If[right[[j]]==left[[j]]+2,
     ReplacePart[currentpart,{left[[j]]->currentpart[[left[[j]]]]-1,right[[j]]->currentpart[[right[[j]]]]-1,left[[j]]+1->currentpart[[left[[j]]+1]]+2}],
     ReplacePart[currentpart,{left[[j]]->currentpart[[left[[j]]]]-1,right[[j]]->currentpart[[right[[j]]]]-1,left[[j]]+1->currentpart[[left[[j]]+1]]+1,right[[j]]-1->currentpart[[right[[j]]-1]]+1}]
     ];
    newstate=FromDigits[newpart,base];
  If[intstatepresent[newstate],Null,Null,
  new=True;inthsdim++;intstate[inthsdim]=newstate;
  intstateposition[newstate]=inthsdim;intstatepresent[newstate]=True;
 ]
 ,{j,1,Length[right]}];
,{i,olddim,currentdim}];
 olddim=currentdim+1;currentdim=inthsdim;
 ];
]

addspinboson[nu_,nd_]:=With[{},
 Clear[Global`state,Global`stateposition,statepresent,Global`b];
  nup=nu;
  ndown=nd;
  sd=spindist[nup,ndown];
  hsdim=0;
  Do[dist=integerdigits[intstate[i1],base,norbitals];orb=orbtopart[dist];Do[newpart=Catch[Module[{res=Table[0,{s,1,2},{i,1,Length[dist]}]},Do[res[[sd[[i,j]],orb[[j]]]]++,{j,1,nup+ndown}];Throw[res]]];
  newstate=spinfromdigits[newpart,base];
  If[statepresent[newstate],Null,Null,hsdim++;Global`state[hsdim]=newstate;Global`stateposition[newstate]=hsdim;statepresent[newstate]=True;]
  ,{i,1,Length[sd]}]
,{i1,1,inthsdim}]
]

spinbosonsqueeze[rp_,nup_,ndown_]:=
 With[{},
  spinbosonsqueezeint[rp];
  addspinboson[nup,ndown];
 Global`norbitals=norbitals;
 Global`nparticles=nparticles;
 Global`base=base;
 Global`hsdim=hsdim;
]



(* spinful fermions *)

spinfermionsqueezeint[rp_]:=With[{},
norbitals=Length[rp];
nparticles=Plus@@rp;
base=3;
Clear[intstate,intstateposition,intstatepresent];
inthsdim=1;
olddim=1;
currentdim=1;
intstate[inthsdim]=FromDigits[rp,base];
intstateposition[intstate[inthsdim]]=inthsdim;
intstatepresent[intstate[inthsdim]]=True;
new=True;
While[new,
new=False;
Do[
currentpart=integerdigits[intstate[i],base,norbitals];
left=Cases[Flatten[Position[currentpart,y_/;y>0]],x_/;x<norbitals-1&&currentpart[[x+1]]!=2];
preright=Cases[Flatten[Position[currentpart,y_/;y>0]],x_/;x>2&&currentpart[[x-1]]!=2];
right=Flatten[Table[Cases[preright,x_/;x>(left[[i]]+If[currentpart[[left[[i]]+1]]==1,2,1]),1,1],{i,1,Length[left]}]];
Do[
newpart=If[right[[j]]>left[[j]]+2,
newpart=ReplacePart[currentpart,{left[[j]]->currentpart[[left[[j]]]]-1,right[[j]]->currentpart[[right[[j]]]]-1,left[[j]]+1->currentpart[[left[[j]]+1]]+1,right[[j]]-1->currentpart[[right[[j]]-1]]+1}],
newpart=ReplacePart[currentpart,{left[[j]]->currentpart[[left[[j]]]]-1,right[[j]]->currentpart[[right[[j]]]]-1,left[[j]]+1->currentpart[[left[[j]]+1]]+2}]
];
newstate=FromDigits[newpart,base];
If[intstatepresent[newstate],Null,Null,
new=True;inthsdim++;intstate[inthsdim]=newstate;
intstateposition[newstate]=inthsdim;intstatepresent[newstate]=True;
];

,{j,1,Length[right]}];

,{i,olddim,currentdim}];
olddim=currentdim+1;currentdim=inthsdim;
];
]

addspinfermion[nu_,nd_]:=With[{},
Clear[Global`state,Global`stateposition,statepresent,Global`b];
base=2;
nup=nu;
ndown=nd;
sd=spindist[nup,ndown];
hsdim=0;
Do[dist=integerdigits[intstate[i1],3,norbitals];orb=orbtopart[dist];Do[newpart=Catch[Module[{res=Table[0,{s,1,2},{i,1,Length[dist]}]},Do[res[[sd[[i,j]],orb[[j]]]]++,{j,1,nup+ndown}];Throw[res]]];
If[FreeQ[newpart,2],
newstate=spinfromdigits[newpart,base];
If[statepresent[newstate],Null,Null,hsdim++;Global`state[hsdim]=newstate;Global`stateposition[newstate]=hsdim;statepresent[newstate]=True;]]
,{i,1,Length[sd]}]
,{i1,1,inthsdim}]
]

spinfermionsqueeze[rp_,nup_,ndown_]:=
 With[{},
 spinfermionsqueezeint[rp];
 addspinfermion[nup,ndown];
 Global`norbitals=norbitals;
 Global`nparticles=nparticles;
 Global`base=base;
 Global`hsdim=hsdim;
]

(** Generating the l^2 operators **)

(* spinless bosons *)

bosonl2:=
With[{},
 Clear[factor];
 Do[factor[i] = Times@@(Factorial/@IntegerDigits[Global`state[i], Global`base, Global`norbitals]),{i,1,Global`hsdim}];
 Do[
  Clear[used, currdata];
  usednumbers = {};
  currstate = IntegerDigits[Global`state[i], Global`base, Global`norbitals];
  midpos = DeleteCases[Flatten[Position[currstate, x_ /; x != 0]], 1];
  midlength = Length[midpos];
  midstates = Table[
   ReplacePart[currstate, {midpos[[i1]] -> currstate[[midpos[[i1]]]] - 1,midpos[[i1]] - 1 -> currstate[[midpos[[i1]] - 1]] + 1}], {i1, 1, midlength}];
  factors = Table[currstate[[midpos[[i1]]]] (midpos[[i1]] - 1), {i1, 1, midlength}];
  Do[
   finpos = DeleteCases[Flatten[Position[midstates[[j]], x_ /; x != 0]], Global`norbitals];
   finlength = Length[finpos];
   finstates = Table[
     ReplacePart[midstates[[j]], {finpos[[i1]] -> midstates[[j, finpos[[i1]]]] - 1, finpos[[i1]] + 1 -> midstates[[j, finpos[[i1]] + 1]] + 1}], {i1, 1, finlength}];
   finfactors = Table[
     midstates[[j, finpos[[i1]]]] (Global`norbitals - finpos[[i1]]), {i1, 1, finlength}];
   Do[
    finstate = FromDigits[finstates[[i1]], Global`base];
    If[statepresent[finstate],
     finnumber = Global`stateposition[finstate];
     If[used[finnumber],
       currdata[finnumber] += factors[[j]] finfactors[[i1]] factor[finnumber]/factor[i],
      Null,
      used[finnumber] = True;
      usednumbers = Append[usednumbers, finnumber];
      currdata[finnumber] = factors[[j]] finfactors[[i1]] factor[finnumber]/factor[i]
      ]
     ]
    , {i1, 1, finlength}];
   , {j, 1, midlength}];
  row[i] = Table[{i, i1} -> currdata[i1], {i1, usednumbers}];
  , {i, 1, Global`hsdim}];
 alldata = Flatten[Table[row[i1], {i1, 1, hsdim}]];
 Clear[used, currdata, usednumbers, row];
 SparseArray[alldata, {hsdim, hsdim}]
 ]

(** Generating L+ conditions **)

(* spinless bosons *)

bosonlplus:=With[{},
Clear[coef,terms,list];
listlength=0;
Do[With[{currentpart=integerdigits[Global`state[i],Global`base,Global`norbitals]},
Do[If[currentpart[[j]]>0,
With[{newpart=ReplacePart[currentpart,{j->currentpart[[j]]-1,j-1->currentpart[[j-1]]+1}],
       newstate=FromDigits[ReplacePart[currentpart,{j->currentpart[[j]]-1,j-1->currentpart[[j-1]]+1}],Global`base]},
 If[terms[newstate],
coef[newstate]+=(j-1)currentpart[[j]]Global`b[i]/Product[currentpart[[k]]!,{k,1,Global`norbitals}],Null,
terms[newstate]=True;
listlength++;
list[listlength]=newstate;
coef[newstate]=(j-1)currentpart[[j]]Global`b[i]/Product[currentpart[[k]]!,{k,1,Global`norbitals}]];
];
];
,{j,2,Global`norbitals}];
];
,{i,1,Global`hsdim}];
Global`listlength=listlength;
]

(* spinless fermions *)

fermionlplus:=With[{},
Clear[coef,terms,list];
listlength=0;
Do[With[{currentpart=integerdigits[Global`state[i],Global`base,Global`norbitals]},
Do[If[currentpart[[j]]>0&&currentpart[[j-1]]==0,
With[{newpart=ReplacePart[currentpart,{j->currentpart[[j]]-1,j-1->currentpart[[j-1]]+1}],
       newstate=FromDigits[ReplacePart[currentpart,{j->currentpart[[j]]-1,j-1->currentpart[[j-1]]+1}],Global`base]},
 If[terms[newstate],
coef[newstate]+=(j-1)currentpart[[j]]Global`b[i]/Product[currentpart[[k]]!,{k,1,Global`norbitals}],Null,
terms[newstate]=True;
listlength++;
list[listlength]=newstate;
coef[newstate]=(j-1)currentpart[[j]]Global`b[i]/Product[currentpart[[k]]!,{k,1,Global`norbitals}]];
];
];
,{j,2,Global`norbitals}];
];
,{i,1,Global`hsdim}];
Global`listlength=listlength;
]

(* spinful bosons *)

spinbosonlplus:=With[{},
Clear[coef,terms,list];
listlength=0;
Do[With[{currentpart=spinintegerdigits[Global`state[i],Global`base,Global`norbitals]},Do[If[currentpart[[type,j]]>0,With[{newpart=ReplacePart[currentpart,{{type,j}->currentpart[[type,j]]-1,{type,j-1}->currentpart[[type,j-1]]+1}],newstate=spinfromdigits[ReplacePart[currentpart,{{type,j}->currentpart[[type,j]]-1,{type,j-1}->currentpart[[type,j-1]]+1}],Global`base]},If[terms[newstate],coef[newstate]+=(j-1)currentpart[[type,j]]Global`b[i]/Product[currentpart[[t,k]]!,{k,1,Global`norbitals},{t,1,2}],Null,
terms[newstate]=True;
listlength++;
list[listlength]=newstate;
coef[newstate]=(j-1)currentpart[[type,j]]Global`b[i]/Product[currentpart[[t,k]]!,{k,1,Global`norbitals},{t,1,2}]]]],{j,2,Global`norbitals},{type,1,2}]],{i,1,Global`hsdim}];
Global`listlength=listlength;
]

(* spinful fermions *)

spinfermionlplus:=With[{},
Clear[coef,terms,list];
listlength=0;
Do[With[{currentpart=spinintegerdigits[Global`state[i],Global`base,Global`norbitals]},Do[If[currentpart[[type,j]]>0&&currentpart[[type,j-1]]==0,With[{newpart=ReplacePart[currentpart,{{type,j}->currentpart[[type,j]]-1,{type,j-1}->currentpart[[type,j-1]]+1}],newstate=spinfromdigits[ReplacePart[currentpart,{{type,j}->currentpart[[type,j]]-1,{type,j-1}->currentpart[[type,j-1]]+1}],2]},If[terms[newstate],coef[newstate]+=(j-1)currentpart[[type,j]]Global`b[i]/Product[currentpart[[t,k]]!,{k,1,Global`norbitals},{t,1,2}],Null,
terms[newstate]=True;
listlength++;
list[listlength]=newstate;
coef[newstate]=(j-1)currentpart[[type,j]]Global`b[i]/Product[currentpart[[t,k]]!,{k,1,Global`norbitals},{t,1,2}]]]],{j,2,Global`norbitals},{type,1,2}]],{i,1,Global`hsdim}];
Global`listlength=listlength;
]

(** Generating S+ conditions **)

(* spinful bosons *)

spinbosonsplus:=With[{},
Do[With[{currentpart=spinintegerdigits[Global`state[i],Global`base,Global`norbitals],uppos=Sort[Flatten[Position[spinintegerdigits[Global`state[i],Global`base,Global`norbitals][[1]],x_/;x!=0]],Greater],downpos=Sort[Flatten[Position[spinintegerdigits[Global`state[i],Global`base,Global`norbitals][[2]],x_/;x!=0]],Greater]},
Do[With[{newpart=ReplacePart[currentpart,{{1,downpos[[j]]}->currentpart[[1,downpos[[j]]]]+1,{2,downpos[[j]]}->currentpart[[2,downpos[[j]]]]-1}],newstate=spinfromdigits[ReplacePart[currentpart,{{1,downpos[[j]]}->currentpart[[1,downpos[[j]]]]+1,{2,downpos[[j]]}->currentpart[[2,downpos[[j]]]]-1}],Global`base]},
If[terms[newstate],
coef[newstate]+=(currentpart[[1,downpos[[j]]]]+1)Global`b[i],Null,
terms[newstate]=True;
listlength++;
list[listlength]=newstate;
coef[newstate]=(currentpart[[1,downpos[[j]]]]+1)Global`b[i];
];];
,{j,1,Length[downpos]}];
],{i,1,Global`hsdim}];
Global`listlength=listlength;
]

(* spinful fermions *)

spinfermionsplus:=With[{},
Do[With[{currentpart=spinintegerdigits[Global`state[i],Global`base,Global`norbitals],uppos=Sort[Flatten[Position[spinintegerdigits[Global`state[i],Global`base,Global`norbitals][[1]],x_/;x!=0]],Greater],downpos=Sort[Flatten[Position[spinintegerdigits[Global`state[i],Global`base,Global`norbitals][[2]],x_/;x!=0]],Greater]},
Do[If[Not[MemberQ[uppos,downpos[[j]]]],With[{newpart=ReplacePart[currentpart,{{1,downpos[[j]]}->1,{2,downpos[[j]]}->0}],newstate=spinfromdigits[ReplacePart[currentpart,{{1,downpos[[j]]}->1,{2,downpos[[j]]}->0}],Global`base],sign=Position[Sort[Append[uppos,downpos[[j]]],Greater],downpos[[j]]][[1,1]]},
If[terms[newstate],
coef[newstate]+=-(-1)^(j+sign)Global`b[i],Null,
terms[newstate]=True;
listlength++;
list[listlength]=newstate;
coef[newstate]=-(-1)^(j+sign)Global`b[i];
];];]
,{j,1,Length[downpos]}];
],{i,1,Global`hsdim}];
Global`listlength=listlength;
]


(** Solving the constraints **)

(* Solve the L+ (and S+ if they are present) conditions exactly *)

solveall:=Catch[With[{},sol=Quiet[Solve[Flatten[{Table[coef[list[i]]==0,{i,1,listlength}]}],Table[Global`b[i],{i,1,Global`hsdim}]],{Solve::svars}][[1]];bfound=Table[sol[[i,1,1]],{i,1,Length[sol]}];
Do[Global`b[bfound[[i]]]=sol[[i,2]],{i,1,Length[sol]}];
bs=Table[i,{i,1,hsdim}];
btbd=Complement[bs,bfound];
Throw[Length[btbd]];]]

(* Solve the L+ (and S+ if they are present) conditions, as well as the extra
conditions encoded via `condlist' exactly. *)

solveallcond:=Catch[With[{},sol=Quiet[Solve[Flatten[{Table[coef[list[i]]==0,{i,1,listlength}],Table[Global`b[Global`condlist[[i]]]==0,{i,1,Length[Global`condlist]}]}],Table[Global`b[i],{i,1,hsdim}]],{Solve::svars}][[1]];bfound=Table[sol[[i,1,1]],{i,1,Length[sol]}];
Do[Global`b[bfound[[i]]]=sol[[i,2]],{i,1,Length[sol]}];
bs=Table[i,{i,1,hsdim}];
btbd=Complement[bs,bfound];
Throw[Length[btbd]];]]

(* Solve the L+ (and S+ if they are present) conditions numerically *)

solveallnum:=Catch[With[{},sol=NSolve[Flatten[{Table[coef[list[i]]==0,{i,1,listlength}]}],Table[Global`b[i],{i,1,Global`hsdim}]][[1]];bfound=Table[sol[[i,1,1]],{i,1,Length[sol]}];
Do[Global`b[bfound[[i]]]=Chop[sol[[i,2]]],{i,1,Length[sol]}];
bs=Table[i,{i,1,Global`hsdim}];
btbd=Complement[bs,bfound];
Throw[Length[btbd]];]]

(* Solve the L+ (and S+ if they are present) conditions, as well as the extra
conditions encoded via `condlist' numerically. *)

solveallcondnum:=Catch[With[{},sol=NSolve[Flatten[{Table[coef[list[i]]==0,{i,1,listlength}],Table[Global`b[Global`condlist[[i]]]==0,{i,1,Length[Global`condlist]}]}],Table[Global`b[i],{i,1,Global`hsdim}]][[1]];bfound=Table[sol[[i,1,1]],{i,1,Length[sol]}];
Do[Global`b[bfound[[i]]]=sol[[i,2]],{i,1,Length[sol]}];
bs=Table[i,{i,1,hsdim}];
btbd=Complement[bs,bfound];
Throw[Length[btbd]];]]


(** Generation of explicit polynomials **)
(* Some initialization *)

pi[n_]:=Permutations[Table[i,{i,1,n}]];
p[v_]:=Permutations[v];

mw[rp_]:=With[{vec=Permutations[rp]},Sum[Product[Global`w[var]^(vec[[i,var]]),{var,1,Length[rp]}],{i,1,Length[vec]}]]
mz[rp_]:=With[{vec=Permutations[rp]},Sum[Product[Global`z[var]^(vec[[i,var]]),{var,1,Length[rp]}],{i,1,Length[vec]}]]
mwa[rp_]:=With[{vec=Permutations[rp],sigs=Signature/@Permutations[rp],n=Length[rp]},Sum[(-1)^(n (n-1)/2)sigs[[i]]Product[Global`w[var]^(vec[[i,var]]),{var,1,Length[rp]}],{i,1,Length[vec]}]]
mza[rp_]:=With[{vec=Permutations[rp],sigs=Signature/@Permutations[rp],n=Length[rp]},Sum[(-1)^(n (n-1)/2)sigs[[i]]Product[Global`z[var]^(vec[[i,var]]),{var,1,Length[rp]}],{i,1,Length[vec]}]]

(* The generation of the various polynomials *)

bosonpol:=
Sum[Global`b[i]mz[partition[integerdigits[Global`state[i],Global`base,Global`norbitals]]],{i,1,Global`hsdim}];

fermionpol:=
Sum[Global`b[i]mza[partition[integerdigits[Global`state[i],Global`base,Global`norbitals]]],{i,1,Global`hsdim}];

spinbosonpol:=
Sum[Global`b[i]mw[partition[spinintegerdigits[Global`state[i],Global`base,Global`norbitals][[1]]]]mz[partition[spinintegerdigits[Global`state[i],Global`base,Global`norbitals][[2]]]],{i,1,Global`hsdim}]

spinfermionpol:=
Sum[Global`b[i]mwa[partition[spinintegerdigits[Global`state[i],Global`base,Global`norbitals][[1]]]]mza[partition[spinintegerdigits[Global`state[i],Global`base,Global`norbitals][[2]]]],{i,1,Global`hsdim}]

(** Generation of root configurate, for `simple' Pauli rules, max n particles in
norb orbitals **)

(* The lz value of a certain orbital configuration *)
lz[vec_]:=With[{norb=Length[vec]},Sum[(-(norb-1)/2-1+i)vec[[i]],{i,1,norb}]]

(* The minimal number of flux needed for a particular Pauli rule *)
nfluxneeded[0,cellsize_,maxincell_]:=0;
nfluxneeded[Ne_/;Ne>0,cellsize_,maxincell_]:=cellsize( Ceiling[Ne/(maxincell)]-1)

(* The first possible configuration *)
firstconfig[Ne_,nflux_,cellsize_,maxincell_]:=If[nflux<nfluxneeded[Ne,cellsize,maxincell],{},Table[If[Mod[i,cellsize]==0,Max[Min[maxincell,(Ne-i/cellsize  maxincell)],0],0],{i,0,nflux}]]

(* Generation of the root configurations *)
rootconfigs[Ne_,nflux_,cellsize_,maxincell_]:=
With[{fc=firstconfig[Ne,nflux,cellsize,maxincell]},
If[fc=={},{},
maxp=partition[Reverse[fc]];
res={fc};
While[
lz[res[[-1]]]<-1/2,
tobeincreased=Position[partition[res[[-1]]]-maxp,x_/;x<0,1,1][[1,1]];
res=Append[res,orbocc[ReplacePart[partition[res[[-1]]],tobeincreased->partition[res[[-1]]][[tobeincreased]]+1],nflux]];
];
res
]]



End[ ]

EndPackage[ ]




