text:
#
# ---------- Group Arithmetics ----------
#

InstallMethod(Subword,
              "Subword of a unipotent element by total degree of coefficients",
              [IsUnipotent,IsPosInt],
              function(u,degree)
              local sys,
              coeffs,poly_degree;

              sys:=chevalleyAdj(u);
              if not IsPolynomialRing(ring(sys)) then
                  Error("Subword: The ring is not a polynomial ring to evaluate!\n");
              fi;

              coeffs:=Filtered(coefficients(u),i->TotalDegree(i[2])=degree);

              return Unipotent(sys,coeffs,Ordering(u));
end);

#
# ---------- Algebraic Operations ----------
#

InstallMethod(Value,
              "Evaluation of the parameter polynomials",
              [IsUnipotent,IsList,IsList],
              function(u,vars,vals)
              local sys,
              lista,i,len;

              sys:=chevalleyAdj(u);
              if not IsPolynomialRing(ring(sys)) then
                  Error("Value: The ring is not a polynomial ring to evaluate.");
              fi;

              lista:=List(coefficients(u),i->[i[1],One(ring(sys))*Value(i[2],vars,vals)]);
              i:=1;
              len:=Length(lista);
              while i <= len do
                  if lista[i][2]=Zero(ring(sys)) then
                      Remove(lista,i);
                      len:=len-1;
                  else i:=i+1; fi;
              od;

              return Unipotent(sys,lista,Ordering(u));
              #return Unipotent(chevalleyAdj(sys),lista,Ordering(u)); # which one is better? I multiply by One(ring(sys))
end);

InstallMethod(Descend,
              "Changes the parameters from the polynomial ring to the coefficients ring",
              [IsUnipotent],
              function(u)
              local sys,
              check;

              sys:=chevalleyAdj(u);
              if not IsPolynomialRing(ring(sys)) then
                  Error("The ring is not a polynomial ring to evaluate.");
              fi;

              if IsAlgebraicU(sys) then sys:=chevalleyAdj(u);
              else sys:=ChevalleyAdj(sys,CoefficientsRing(ring(sys))); fi;

              check:=List(coefficients(u),i->LeadingMonomial(i[2]));
              check:=List(check,i->Sum(i{2*[1..Length(i)/2]}));
              check:=Filtered(check,i->i<>0);
              if check=[] then
                  return Unipotent(sys,
                                   List(coefficients(u),i->[i[1],LeadingCoefficient(i[2])]),
                                   Ordering(u),
                                   Order(u));
              fi;

              return fail;
end);

InstallMethod(Descend,
              "If possible changes the parameters from the polynomial ring to the coefficients ring",
              [IsUnipotent,IsChevalleyAdj],
              function(u,sys)
              local sysu,
              check;

              sysu:=chevalleyAdj(u);

              if not IsPolynomialRing(ring(sysu)) then
                  Error("The ring is not a polynomial ring to evaluate.");
              elif CoefficientsRing(ring(sysu))<>ring(sys) then
                  Error("Cannot descend unipotent to given simple adjoint system");
              fi;

              check:=List(coefficients(u),i->LeadingMonomial(i[2]));
              check:=List(check,i->Sum(i{2*[1..Length(i)/2]}));
              check:=Filtered(check,i->i<>0);
              if check=[] then
                  return Unipotent(sys,
                                   List(coefficients(u),i->[i[1],LeadingCoefficient(i[2])]),
                                   Ordering(u),
                                   Order(u));
              fi;

              return u;
end);

InstallMethod(Ascend,
              "If possible embed Unipotent into AlgebraicU",
              [IsUnipotent,IsAlgebraicU],
              function(u,sys)
              local sysu,
              check;

              if chevalleyAdj(u)<>chevalleyAdj(sys) then
                  Error("Cannot ascend Unipotent to AlgebraicU: different chevalleyAdj!\n");
              fi;

              return Unipotent(sys,
                               List(coefficients(u),i->[i[1],One(ring(sys))*i[2]]),
                               Ordering(u),
                               Order(u));
end);

InstallMethod(Ascend,
              "If possible embed the parameters in a polynomial ring with coefficients ring the base ring",
              [IsUnipotent,IsChevalleyAdj],
              function(u,sys)
              local sysu,
              check;

              sysu:=chevalleyAdj(u);

              if not IsPolynomialRing(ring(sys)) then
                  Error("The new ring is not a polynomial ring to embed in.");
              elif CoefficientsRing(ring(sys))<>ring(sysu) then
                  Error("CoefficientsRing and ring(sysu) do not coincide.");
              fi;

              return Unipotent(sys,
                               List(coefficients(u),i->[i[1],One(ring(sys))*i[2]]),
                               Ordering(u),
                               Order(u));
end);

InstallMethod(Ascend,
              "If possible embed the parameters in a polynomial ring with coefficients ring the base ring",
              [IsUnipotent,IsPolynomialRing],
              function(u,inel)
              local sys,
              check;

              sys:=chevalleyAdj(u);

              if not IsPolynomialRing(inel) then
                  Error("The new ring is not a polynomial ring to embed in.");
              elif CoefficientsRing(inel)<>ring(sys) then
                  Error("CoefficientsRing and ring(sys) do not coincide.");
              fi;

              return Unipotent(ChevalleyAdj(sys,inel),
                               List(coefficients(u),i->[i[1],One(inel)*i[2]]),
                               Ordering(u),
                               Order(u));
end);

InstallMethod(UnipotentCentralizer,
              "Centralizer of unipotent element in algebraicU",
              [IsUnipotent],
              function(u)
              local result,
              sys,pr,pr_len,uu,
              APR,avars,avars_free,generic,
              relations,sol,torsion,
              i,new_indices;

              sys:=chevalleyAdj(u);
              if not IsAlgebraicU(sys) then
                  sys:=AlgebraicU(sys);
                  u:=Ascend(u,sys);
              fi;
              pr:=positiveRoots(sys);
              pr_len:=Length(pr);
              APR:=ring(sys);
              avars:=IndeterminatesOfPolynomialRing(APR){[1..pr_len]};

              generic:=Unipotent(sys,List([1..pr_len],i->[i,avars[pr_len-i+1]]),[1..pr_len],0);

              result:=List([1..pr_len],i->avars[i]);

              torsion:=[];
              for i in [2..pr_len] do
                  #relations:=coefficients(CommMod(generic,u,i));
              relations:=coefficients(CommMod(u,generic,i));   #!!!ATENTIE CommMod(a,b,i) foloseste Ordering lui b care tre sa fie [1..pr_len]
                  relations:=List(relations,i->i[2]);
              #Print("\nsolving ",i,": ",relations,"\n");
              #if i=8 then Error("!"); fi;
                  sol:=SolveRelations(relations,avars,result,APR);
                  if sol[2]<>[] then Print(" >>>",i," UnipotentCentralizer: could not solve all relations!",relations,"\n"); fi;
                  Add(torsion,sol[3]);
                  result:=sol[1];
                  result:=List(result,i->i*One(APR));
              #Print("result :",result,"\n");
                  generic:=Unipotent(sys,List([1..pr_len],i->[i,result[pr_len-i+1]]),[1..pr_len],0);
              od;

              avars_free:=result{Filtered([1..pr_len],i->result[i] in avars)};
              avars_free:=avars_free{Filtered([1..Length(avars_free)],i->not avars_free[i] in avars_free{[1..i-1]})};

              new_indices:=Reversed([1..Length(avars_free)]);
              result:=List(result,i->Value(i,avars_free,avars{new_indices})*One(APR));
              generic:=Unipotent(sys,List([1..pr_len],i->[i,result[pr_len-i+1]]),[1..pr_len],0);

              return [generic,u,Concatenation(torsion),relations];
end);

#
#
#  THIS DOES NOT WORK
#
#
InstallMethod(UnipotentCentralizer2,
              "Centralizer of unipotent element in algebraicU",
              [IsUnipotent],
              function(u)
              local result,
              sys,pr,pr_len,uu,
              APR,avars,avars_free,generic,
              relations,sol,torsion,
              root_support,other_roots,
              i,new_indices;

              sys:=chevalleyAdj(u);
              if not IsAlgebraicU(sys) then
                  sys:=AlgebraicU(sys);
                  u:=Ascend(u,sys);
              fi;
              pr:=positiveRoots(sys);
              pr_len:=Length(pr);
              APR:=ring(sys);
              avars:=IndeterminatesOfPolynomialRing(APR){[1..pr_len]};

              root_support:=[];
              for i in [1..pr_len] do
                  generic:=Unipotent(sys,[[i,avars[1]]]);
                  generic:=Conj(generic,u);
                  if Length(coefficients(generic))=1 and coefficients(generic)[1][2]=avars[1] then Add(root_support,i); fi;
              od;
              other_roots:=Difference([1..pr_len],root_support);

              #generic:=Unipotent(sys,List([1..pr_len],i->[i,avars[pr_len-i+1]]),[1..pr_len],0);
              #generic:=Subword(generic,other_roots);
              generic:=Unipotent(sys,List(other_roots,i->[i,avars[pr_len-i+1]]),[1..pr_len],0);

              result:=List([1..pr_len],i->avars[i]);
              for i in root_support do result[pr_len-i+1]:=Zero(APR); od;

              torsion:=[];
              for i in [2..pr_len] do
                  #relations:=coefficients(CommMod(generic,u,i));
              relations:=coefficients(CommMod(u,generic,i));   #!!!ATENTIE CommMod(a,b,i) foloseste Ordering lui b care tre sa fie [1..pr_len]
                  relations:=List(relations,i->i[2]);
              #Print("\nsolving ",i,": ",relations,"\n");
              #if i=8 then Error("!"); fi;
                  sol:=SolveRelations(relations,avars,result,APR);
                  if sol[2]<>[] then Print(" >>>",i," UnipotentCentralizer: could not solve all relations!",relations,"\n"); fi;
                  Add(torsion,sol[3]);
                  result:=sol[1];
                  result:=List(result,i->i*One(APR));
              #Print("result :",result,"\n");
                  generic:=Unipotent(sys,List([1..pr_len],i->[i,result[pr_len-i+1]]),[1..pr_len],0);
              od;
              Error("!");
              for i in root_support do result[pr_len-i+1]:=avars[pr_len-i+1]; od;

              avars_free:=result{Filtered([1..pr_len],i->result[i] in avars)};
              avars_free:=avars_free{Filtered([1..Length(avars_free)],i->not avars_free[i] in avars_free{[1..i-1]})};

              new_indices:=Reversed([1..Length(avars_free)]);
              result:=List(result,i->Value(i,avars_free,avars{new_indices})*One(APR));
              generic:=Unipotent(sys,List([1..pr_len],i->[i,result[pr_len-i+1]]),[1..pr_len],0);

              return [generic,u,Concatenation(torsion),relations];
end);

InstallMethod(TangentSpace,
              "Returns a basis of Nilpotents representing for the tangent space",
              [IsUnipotent],
              function(u)
              local result,
              sysAPR,APR,avars,
              dim,BPR,xvars,
              order,exp,vals,tmp,tmp_len,i,j,monoms;

              sysAPR:=chevalleyAdj(u);
              APR:=ring(sysAPR);
              avars:=IndeterminatesOfPolynomialRing(APR){[1..Length(positiveRoots(sysAPR))]};
              dim:=Length(Set(List(Filtered(coefficients(u),i->i[2] in avars),j->j[2])));
              avars:=avars{[1..dim]};

              vals:=List([1..dim],i->Zero(APR));
              result:=[];
              for i in [1..dim] do
                  vals[i]:=avars[1];
                  tmp:=Value(u,avars,vals);
                  #SetOrder(tmp,Order(u));
                  tmp:=ShallowCopy(coefficients(tmp));
                  tmp_len:=Length(tmp);
                  j:=1;
                  while j <= tmp_len do
                       monoms:=ExtRepPolynomialRatFun(tmp[j][2]);
                       monoms:=List(2*[1..Length(monoms)/2]-1,i->[monoms[i],monoms[i+1]]);
                       monoms:=Filtered(monoms,i->i[1][2]=1);
                       if monoms<>[] then
                          tmp[j]:=[tmp[j][1],monoms[1][2]*One(APR)];
                          j:=j+1;
                      else
                          Remove(tmp,j);
                          tmp_len:=tmp_len-1;
                      fi;
                  od;
                  Add(result,tmp);
                  vals[i]:=Zero(APR);
              od;

              result:=List(result,i->NilpotentChv(sysAPR,i));
              return result;
end);

InstallMethod(Zentrum,
              "Centre of connected unipotent algebraic group",
              [IsUnipotent],
              function(u)
              local result,
              APR,avars,
              relations,
              tang,adu_res,unu;

              APR:=ring(chevalleyAdj(u));
              avars:=IndeterminatesOfPolynomialRing(APR);
              avars:=Filtered(coefficients(u),i->i[2] in avars);
              avars:=Set(List(avars,i->i[2]));

              tang:=TangentSpace(u);
              #Error("!");
              adu_res:=RestrictedAdu(u,tang);

              unu:=DiagonalMat(List([1..Length(avars)],i->One(APR)));

              relations:=Concatenation(adu_res-unu);
              relations:=Filtered(relations,i->i<>Zero(APR));

#return relations;

              #Print(avars," with ",relations,"\n");

              result:=SolveRelations(relations,avars,avars,APR);

              if result[2]<>[] then Print("Zentrum: could not solve all relations! ",result[2],"\n"); fi;

              return Value(u,avars,result[1]);
end);

InstallMethod(ParameterPieces,
              "Subvarieties for each parameter",
              [IsUnipotent],
              function(u)
              local result,
              sysAPR,APR,avars,avar,
              dim,BPR,xvars,
              order,exp,vals,tmp,tmp_len,i,j,monoms;

              sysAPR:=chevalleyAdj(u);
              APR:=ring(sysAPR);
              avars:=IndeterminatesOfPolynomialRing(APR){[1..Length(positiveRoots(sysAPR))]};
              avar:=avars[1];
              avars:=Set(List(Filtered(coefficients(u),i->i[2] in avars),j->j[2]));
              Sort(avars,function(a,b) return ExtRepPolynomialRatFun(a)[1][1]<ExtRepPolynomialRatFun(b)[1][1]; end);
              dim:=Length(avars);
              #avars:=avars{[1..dim]};

              vals:=List([1..dim],i->Zero(APR));
              result:=[];
              for i in [1..dim] do
                  vals[i]:=avar;#avars[1];
                  tmp:=Value(u,avars,vals);
                  #Print(coefficients(tmp));
                  #SetOrder(tmp,Order(u));
                  tmp:=ShallowCopy(coefficients(tmp));
                  tmp_len:=Length(tmp);
                  Add(result,tmp);
                  vals[i]:=Zero(APR);
              od;

              result:=List(result,i->Unipotent(sysAPR,i,Ordering(u),-2));
              return result;
end);

InstallMethod(Dimension,
              "Dimension of subvariety in unipotent group",
              [IsUnipotent],
              function(u)
              local avars_free;

              avars_free:=List(coefficients(u),i->i[2]);
              avars_free:=List(avars_free,i->ExtRepPolynomialRatFun(i));
              avars_free:=Concatenation(avars_free);
              avars_free:=avars_free{2*[1..Length(avars_free)/2]-1};
              avars_free:=Concatenation(avars_free);
              avars_free:=avars_free{2*[1..Length(avars_free)/2]-1};

              return Length(Set(avars_free));
end);

#
# ---------- Witt Group Operations ----------
#

InstallMethod(GenericWittGroup,
              "Witt group with generic polynomials along the unipotent group gr",
              [IsUnipotent,IsWitt,IsList],
              function(gr,w,values)
              local result,
              coeffs,
              APR,avars,avars_free,avars_free_len,avars_ext,
              BPR,bvars,
              CPR,cvarnames,cvars,
              poly,monom,i,j,k,pos;

              coeffs:=coefficients(gr);

              # APR's
              APR:=ring(chevalleyAdj(gr));
              avars:=IndeterminatesOfPolynomialRing(APR);
              avars_free:=Filtered(coeffs,i->i[2] in avars);
              avars_free:=Set(List(avars_free,i->i[2]));
              Sort(avars_free,function(a,b) return ExtRepPolynomialRatFun(a)[1][1]<ExtRepPolynomialRatFun(b)[1][1]; end);
              avars_free_len:=Length(avars_free);
              avars_ext:=List(avars_free,i->ExtRepPolynomialRatFun(i));
              avars_ext:=List(avars_ext,i->i[1][1]);

              # BPR's
              BPR:=ring(chevalleyAdjBPR(w));

              # CPR's
              cvarnames:=List([1..avars_free_len],i->Concatenation("c_",String(i)));
              CPR:=PolynomialRing(BPR,cvarnames);
              cvars:=IndeterminatesOfPolynomialRing(CPR);

              coeffs:=List(coeffs,i->[i[1],List(ExtRepPolynomialRatFun(i[2]))]);
              result:=[];
              for i in [1..Length(coeffs)] do
                  poly:=Zero(CPR);
                  for j in [1..Length(coeffs[i][2])/2] do
                      monom:=One(CPR);
                      for k in [1..Length(coeffs[i][2][2*j-1])/2] do
                          pos:=Position(avars_ext,coeffs[i][2][2*j-1][2*k-1]);
                          monom:=monom*cvars[pos]^coeffs[i][2][2*j-1][2*k];
                      od;
                      monom:=monom*((coeffs[i][2][2*j]*One(APR))*One(BPR));
                      poly:=poly+monom;
                  od;
                  Add(result,[coeffs[i][1],poly]);
              od;

              result:=List(result,i->[i[1],Value(i[2],cvars,values)]);

              return Unipotent(chevalleyAdjBPR(w),result,Ordering(gr),-2);
end);

InstallMethod(SolveWittRoot,
              "Solution for witt group construction at certain root",
              [IsUnipotent,IsWitt,IsPosInt,IsList,IsInt],
              function(wx,w,root,values,step)
              local result,
              APR,avars,avars_free,avars_free_len,
              BPR,dim,xvar,yvar,wy,
              gws,gen_poly,gen_poly_len,xy_poly,total_degree,relations,i;

              gws:=GenericSum(w);

              # BPR's
              BPR:=ring(chevalleyAdjBPR(w));
              dim:=Dimension(w);
              xvar:=IndeterminatesOfPolynomialRing(BPR){[1..dim]};
              yvar:=IndeterminatesOfPolynomialRing(BPR){dim+[1..dim]};

              wy:=Value(wx,xvar,yvar);
              xy_poly:=coefficients(Subword(wx*wy,[root]));
              if xy_poly=[] then
                  return [Zero(BPR),values,step];
              else
                  xy_poly:=xy_poly[1][2];
              fi;

              # APR's
              APR:=ring(algebraicU(w));
              avars:=IndeterminatesOfPolynomialRing(APR);

              result:=ShallowCopy(values);

              # construct generic polynomial
              total_degree:=TotalDegree(xy_poly);
              gen_poly:=Concatenation(GenericLevelPolynomials(w){[1..total_degree]});
              gen_poly_len:=Length(gen_poly);
              gen_poly:=Sum(List([1..gen_poly_len],i->result[step+i]*gen_poly[i]));

              # contruct relations
              relations:=Value(gen_poly,xvar,gws)-gen_poly-Value(gen_poly,xvar,yvar)-xy_poly;
              relations:=ExtRepPolynomialRatFun(relations);
              relations:=relations{2*[1..Length(relations)/2]};

              # solve relations
              result:=SolveRelations(relations,avars,result,APR);
              if result[2]<>[] then Error("SolveWittRoot: coulde not solve all relations!"); return fail; fi;
              result:=result[1];

              # rearrange free variables
              avars_free:=Filtered([1..gen_poly_len+step],i->result[i] in avars);
              avars_free:=result{avars_free};
              avars_free_len:=Length(avars_free);
#    result:=List(result,i->i*One(APR));
#    result:=List(result,i->Value(i,avars_free,avars{[1..avars_free_len]})*One(APR));

              gen_poly:=Concatenation(GenericLevelPolynomials(w){[1..total_degree]});
              gen_poly:=Sum(List([1..Length(gen_poly)],i->result[i+step]*gen_poly[i]));

              return [gen_poly,result,step+gen_poly_len];
end);

InstallMethod(WittGroup,
              "Witt group of unipotent inside group",
              [IsUnipotent,IsUnipotent,IsWitt],
              function(u,centralizer,w)
              local result,u_ordered,
              APR,avars,avars_free,avars_ext,
              BPR,dim,xvars,xvars_ext,
              coeffs,
              i,j,k,pos,sol,step,values,update_poly,monom,
              wg;

              # APR's
              #APR:=ring(algebraicU(w)); asta de ce nu merge?
              APR:=ring(chevalleyAdj(u));
              avars:=IndeterminatesOfPolynomialRing(APR);

              # BPR's
              BPR:=ring(chevalleyAdjBPR(w));
              dim:=Dimension(w);
              xvars:=IndeterminatesOfPolynomialRing(BPR){[1..dim]};

              # centralizer
              #centralizer:=UnipotentCentralizer(algebraicAdj(w),u,distP);
              avars_free:=Set(List(Filtered(coefficients(centralizer),i->i[2] in avars),j->j[2]));
              Sort(avars_free,function(a,b) return ExtRepPolynomialRatFun(a)[1][1]<ExtRepPolynomialRatFun(b)[1][1]; end);
              avars_ext:=List(avars_free,i->ExtRepPolynomialRatFun(i));
              avars_ext:=List(avars_ext,i->i[1][1]);

              result:=List(avars_free,i->Zero(BPR));
              coeffs:=coefficients(centralizer);
              u_ordered:=Unipotent(chevalleyAdj(u),#algebraicU(w), asta de ce nu merge?
                                   coefficients(u),#coefficients(Ascend(u,simpleAdjAPR(algebraicAdj(w)))),
                                   [1..Length(positiveRoots(chevalleyAdjBPR(w)))],
                                   -2);
              for i in coefficients(u_ordered) do
                  k:=Filtered(coeffs,j->j[1]=i[1]);
                  if k = [] then return fail; fi;
                  k:=k[1];
                  if k[2] in avars_free then
                      pos:=Position(avars_free,k[2]);
                      result[pos]:=xvars[1];
                  fi;
              od;
              wg:=GenericWittGroup(centralizer,w,result);

              # function for handeling free variables udpate
              xvars_ext:=List(xvars,i->ExtRepPolynomialRatFun(i));
              xvars_ext:=List(xvars_ext,i->i[1][1]);
              update_poly:=function(poly,values)
                  local result,
                  monom,ext,j,k,pos;
                  result:=Zero(BPR);
                  ext:=ExtRepPolynomialRatFun(poly);
                  for j in [1..Length(ext)/2] do
                      monom:=One(BPR);
                      for k in [1..Length(ext[2*j-1])/2] do
                          pos:=Position(xvars_ext,ext[2*j-1][2*k-1]);
                          monom:=monom*xvars[pos]^ext[2*j-1][2*k];
                      od;
                      monom:=monom*(Value(ext[2*j],avars,values)*One(APR));
                      result:=result+monom;
                  od;
                  return result;
              end;

              values:=ShallowCopy(avars);
              step:=0;
              for i in coeffs do
                  if i[2] in avars_free then
                      pos:=Position(avars_free,i[2]);
                      if result[pos]=Zero(BPR) then
                          # value of free variable
                          sol:=SolveWittRoot(wg,w,i[1],values,step);
                          #if sol=fail then return fail; fi;
                          values:=sol[2];
                          step:=sol[3];
                          # update the construction made so far
                          result:=List(result,i->update_poly(i,values));
                          # set free variable
                          result[pos]:=sol[1];
                          # rebuild group
                          wg:=GenericWittGroup(centralizer,w,result);
                      fi;
                  fi;
              od;

              return wg;
end);

InstallMethod(TestWittGroup,
              "Test if the given unipotent element has the structure of a given witt group",
              [IsUnipotent,IsWitt],
              function(ux,w)
              local dim,p,BPR,xvar,yvar,
              uy,uxy,xcoeff,xycoeff,
              gws,i;

              dim:=Dimension(w);
              BPR:=ring(chevalleyAdjBPR(w));
              xvar:=IndeterminatesOfPolynomialRing(BPR){[1..dim]};
              yvar:=IndeterminatesOfPolynomialRing(BPR){dim+[1..dim]};

              uy:=Value(ux,xvar,yvar);
              gws:=GenericSum(w);

              uxy:=ux*uy;
              for i in [1..Length(positiveRoots(chevalleyAdj(ux)))] do
                  xcoeff:=coefficients(Subword(ux,[i]));
                  xycoeff:=coefficients(Subword(uxy,[i]));
                  if xcoeff=[] and xycoeff<>[] then return i;
                  elif xcoeff<>[] and xycoeff=[] then return i;
                  elif xcoeff<>[] and xycoeff<>[] and Value(xcoeff[1][2],xvar,gws)<>xycoeff[1][2] then return i;
                  fi;
              od;

              return true;
end);

InstallMethod(WittDecomposition,
              "Decomposition of abelian algebraic group into witt group factors",
              [IsUnipotent,IsWitt],
              function(abelian,w)
              local result,
              APR,avars,p,vals,i,tmp,dim,exp,
              handled,support;

              APR:=ring(chevalleyAdj(abelian));
              avars:=IndeterminatesOfPolynomialRing(APR);
              avars:=Filtered(coefficients(abelian),i->i[2] in avars);
              avars:=Set(List(avars,i->i[2]));
              Sort(avars,function(a,b) return ExtRepPolynomialRatFun(a)[1][1]<ExtRepPolynomialRatFun(b)[1][1]; end);
              dim:=Length(avars);
              #    Print("dim=",dim,"\n");

              p:=Characteristic(APR);
              exp:=Maximum(List(ParameterPieces(abelian),i->OrderOp(i)));
              exp:=LogInt(exp,p);

              vals:=List([1..dim],i->Zero(APR));
              result:=[];
              handled:=[];
              for i in [1..dim] do
                  vals[i]:=One(APR);#avars[1];
                  tmp:=Value(abelian,avars,vals);
                  support:=List(coefficients(tmp),i->i[1]);
                  if not IsSubset(handled,support) then
                      #tmp:=WittGr(tmp,abelian,w[LogInt(OrderOp(tmp),p)]);
                      tmp:=WittGroup(tmp,abelian,w);
                      Add(result,tmp);
                      support:=List(coefficients(tmp),i->i[1]);
                      handled:=Concatenation(handled,support);
                  fi;
                  vals[i]:=Zero(APR);
              od;

#    result:=List(result,i->UnipotentLieElement(sysAPR,i,0,ordering(u),1));
              return result;
end);


text:
#
# cls 17: A_2
#
handle17char2:=function()
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=AllClasses(orbs)[17];
    info:=infos[17];

    rep:=Representative(o);
    #pr[15]=[0,0,1,1,1,0]
    #pr[13]=[0,1,1,1,0,0]
    #pr[26]=[1,1,1,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[15,13,26]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1]]);
    tmp1:=Conj(tmp1,uu);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[15,13,26]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1]]);
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;

#
# cls 10: D_4(a_1)
#
handle10char2:=function()
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=AllClasses(orbs)[10];
    info:=infos[10];

    rep:=Representative(o);
    #pr[4] =[0,0,0,1,0,0]
    #pr[32]=[1,1,2,2,1,1]
    #pr[33]=[1,1,1,2,2,1]
    tmp1:=ApplyRootsReflections(rep,[4,32,33]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,-1],[5,1]]);
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[4,1,2,6],[-1,-1,-1,-1]);

    #pr[20]=[0,1,0,1,1,1]
    #pr[21]=[0,0,1,1,1,1]
    tmp3:=ApplyRootsReflections(rep,[20,21]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1],[4,1],[5,-1],[9,1],[10,-1]]);
    tmp3:=Conj(tmp3,uu);
    tmp3:=ConjugateByTori(tmp3,[2,4],[-1,-1,-1,-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[4,32,33]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,-1],[5,1]]);
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[4,1,2,6],[-1,-1,-1,-1]);

    tmp4:=ApplyRootsReflections(cent,[20,21]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1],[4,1],[5,-1],[9,1],[10,-1]]);
    tmp4:=Conj(tmp4,uu);
    tmp4:=ConjugateByTori(tmp4,[2,4],[-1,-1,-1,-1]);

    return [tmp1,tmp3,tmp2,tmp4];
end;

#
# cls 4: E_6(a_3) for [1,8,9,11,14,19]
#
handle4char2:=function()
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=AllClasses(orbs)[4];
    info:=infos[4];

    rep:=Representative(o);
    tmp1:=ConjNegUa(rep,3,1);
    tmp1:=ConjNegUa(tmp1,2,1);
    uu:=Unipotent(chevalleyAdj(rep),[[5,1],[6,1],[7,1],[8,1],[13,1],[14,1],[19,1],[20,1],[25,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjNegUa(cent,3,1);
    tmp2:=ConjNegUa(tmp2,2,1);
    uu:=Unipotent(chevalleyAdj(cent),[[5,1],[6,1],[7,1],[8,1],[13,1],[14,1],[19,1],[20,1],[25,1]],Ordering(cent));
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;
#
# THIS IS FOR THE E_6(a_3) class of E_6 in char 2
#
#vals:=ShallowCopy(avars);
#vals[3]:=0;
#vals[2]:=0;
#vals[9]:=vals[1];
#vals[5]:=One(APR);
#vals[4]:=0;
#vals[7]:=vals[1]+One(APR);
#vals[6]:=vals[1]+One(APR);
#vals[8]:=One(APR);
#vals[13]:=vals[7];
#vals[14]:=vals[1]+One(APR);
#vals[15]:=0;
#vals[11]:=vals[1];
#vals[18]:=vals[1]^2+vals[1]+vals[16];
#vals[19]:=One(APR);
#vals[20]:=vals[1]^2+vals[12]+One(APR);
#vals[21]:=vals[12];
#vals[22]:=vals[12];
#vals[24]:=vals[1]*vals[10]+vals[1]+vals[10]+vals[12]+vals[18];
#vals[25]:=vals[1]^2+vals[17]+One(APR);
#vals[27]:=vals[17];
#vals[28]:=vals[1]^2+vals[1]+vals[12]+vals[23];
#vals[29]:=vals[23];
#vals[31]:=vals[1]^2*vals[10]+vals[1]^2+vals[1]*vals[12]+vals[10]+vals[16]+vals[23]+vals[26];
#vals[33]:=vals[1]^4+vals[1]^3+vals[1]^2*vals[16]+vals[1]^2+vals[1]*vals[12]+vals[12]^2+vals[1]+vals[16]+vals[23]+vals[26];
#vals[34]:=vals[1]^4+vals[1]^3+vals[1]^2*vals[16]+vals[1]^2+vals[12]^2+vals[12]*vals[17]+vals[1]+vals[16]+vals[30]+vals[32];
#
#vv:=Value(generic,avars,vals);
#vvv:=Value(vv,[avars[1],avars[10],avars[12],avars[16],avars[17],avars[23],avars[26],avars[30],avars[32],avars[35],avars[36]],[0,0,0,0,0,0,0,0,0,0,0]);

ext: gi


text:
#
# cls 14: \tilde A_1
#
handle14char3:=function()
    local tmp1,tmp2,o,info,rep,cent;
    o:=AllClasses(orbs)[14];
    info:=infos[14];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[12]);#pr[12]=[1,0,2,1]-->[0,1,2,1]
    tmp1:=ConjugateByTorus(tmp1,3,-1);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[12]);#pr[12]=[1,0,2,1]-->[0,1,2,1]
    tmp2:=ConjugateByTorus(tmp2,3,-1);

    return [tmp1,tmp2];
end;

#
# cls 12: A_2
#
handle12char3:=function()
    local tmp1,tmp2,o,info,rep,cent,u;
    o:=AllClasses(orbs)[12];
    info:=infos[12];

    rep:=Representative(o);
    #pr[7]=[0,0,1,1]-->[0,1,1,0]
    #pr[13]=[0,1,2,1]-->[1,1,2,0]
    tmp1:=ApplyRootsReflections(rep,[7,13]);
    u:=Unipotent(algebraicU(orbs),[[4,-1]]);
    tmp1:=Conj(tmp1,u);
    tmp1:=ConjugateByTorus(tmp1,2,-1);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[7,13]);
    tmp2:=Conj(tmp2,u);
    tmp2:=ConjugateByTorus(tmp2,2,-1);

    return [tmp1,tmp2];
end;

#
# cls 9: B_2
#
handle9char3:=function()
    local tmp1,tmp2,o,info,rep,cent;
    o:=AllClasses(orbs)[9];
    info:=infos[9];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[9]);#pr[9]=[0,1,1,1]-->[1,1,1,0]
    tmp1:=ConjugateByTorus(tmp1,4,-1);# merge si cu 1

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[9]);#pr[9]=[0,1,1,1]-->[1,1,1,0]
    tmp2:=ConjugateByTorus(tmp2,4,-1);# merge si cu 1

    return [tmp1,tmp2];
end;

#
# cls 7: C_3(a_1)
#
handle7char3:=function()
    local tmp1,tmp2,o,info,rep,cent;
    o:=AllClasses(orbs)[7];
    info:=infos[7];

    rep:=Representative(o);
    tmp1:=ConjugateByTorus(rep,1,-1);# h_1 --> h_4
    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTorus(cent,1,-1);# h_1 --> h_4
    return [tmp1,tmp2];
end;

#
# cls 6: F_4(a_3) [4,6,10,18] (finished?)
#
handle6char3:=function()
    local tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,o,info,rep,cent,uu;
    o:=AllClasses(orbs)[6];
    info:=infos[6];

    rep:=Representative(o);
    tmp1:=ApplyRootsReflections(rep,[3]);#pr[3]=[0,0,1,0]-->[0,0,1,0]
    uu:=Unipotent(chevalleyAdj(rep),[[2,-1],[6,1]],Ordering(rep));
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[1,2,3,4],[1,-1,aaa,aaa^2]);

    tmp3:=ApplyRootsReflections(rep,[5,1]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1]],Ordering(rep));#[1..24]);
    tmp3:=Conj(tmp3,uu);
    uu:=Unipotent(chevalleyAdj(rep),[[16,-1]],Ordering(rep));
    tmp3:=Conj(tmp3,uu);
    tmp3:=ConjNegUa(tmp3,2,-1);
    uu:=Unipotent(chevalleyAdj(rep),[[13,1],[2,-1],[4,1]],Ordering(rep));
    tmp3:=Conj(tmp3,uu);
    tmp3:=ConjugateByTori(tmp3,[3],[aaa]);

    tmp5:=ConjugateByTori(rep,[3],[-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[3]);#pr[3]=[0,0,1,0]-->[0,0,1,0]
    uu:=Unipotent(chevalleyAdj(cent),[[2,-1],[6,1]]);
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[1,2,3,4],[1,-1,aaa,aaa^2]);

    tmp4:=ApplyRootsReflections(cent,[5,1]);
    uu:=Unipotent(chevalleyAdj(rep),[[2,1]],Ordering(cent));#[1..24]);
    tmp4:=Conj(tmp4,uu);
    uu:=Unipotent(chevalleyAdj(rep),[[16,-1]],Ordering(cent));
    tmp4:=Conj(tmp4,uu);
    tmp4:=ConjNegUa(tmp4,2,-1);
    uu:=Unipotent(chevalleyAdj(rep),[[13,1],[2,-1],[4,1]],Ordering(cent));
    tmp4:=Conj(tmp4,uu);
    tmp4:=ConjugateByTori(tmp4,[3],[aaa]);

    tmp6:=ConjugateByTori(cent,[3],[-1]);

    return [tmp1,tmp3,tmp5,tmp2,tmp4,tmp6];
end;

#
# cls 3: F_4(a_2) (not finished)
#
handle3char3:=function()
    local tmp1,tmp2,o,info,rep,cent;
    o:=AllClasses(orbs)[3];
    info:=infos[3];

    rep:=Representative(o);
    tmp1:=ConjugateByTori(rep,[4],[-1]);# h_1 --> h_4

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTori(cent,[4],[-1]);# h_1 --> h_4
    return [tmp1,tmp2];
end;

#
# cls 2: F_4(a_1)
#
handle2char3:=function()
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=AllClasses(orbs)[2];
    info:=infos[2];

    rep:=Representative(o);
    tmp1:=ConjugateByTori(rep,[3,1],[-1,-1]);# h_1 --> h_4, h_3 --> h_3
    uu:=Unipotent(chevalleyAdj(rep),[[3,-1]],Ordering(rep));#pr[3]=[0,0,1,0]-->[0,0,1,0]
    tmp1:=Conj(tmp1,uu);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTori(cent,[3,1],[-1,-1]);# h_1 --> h_4, h_3 --> h_3
    uu:=Unipotent(chevalleyAdj(cent),[[3,-1]],Ordering(cent));#pr[3]=[0,0,1,0]-->[0,0,1,0]
    tmp2:=Conj(tmp2,uu);

    return [tmp1,tmp2];
end;


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

IsASPAlgebra := x -> IsRec(x) and IsBound(x.isASPAlgebra) and x.isASPAlgebra;

Class(AlgebraOps, PrintOps, rec(
    \= := (b1, b2) -> Cond(
        not IsASPAlgebra(b1) or not IsASPAlgebra(b2), false,
        ObjId(b1)=ObjId(b2) and b1.rChildren() = b2.rChildren())
));

Class(ASPAlgebra, rec(
    isASPAlgebra:=true,
    codeletShape := self >> Concatenation([ObjId(self)], List(self.rChildren(), CodeletShape)),
    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), rch),
    unparse := "RewritableObjectExp",
));

Declare(XN_skew, XN_plus_1, XN_min_1, XN_min_w);

#F XN_min_1(n) -- represents x^n - 1 polynomial
#F   Spectral decomposition: 1d and 2d components
#F   <n> even or odd
#F
Class(XN_min_1, ASPAlgebra, rec(
    __call__ := (self, n, rot) >> let(nn:=_unwrap(n), Checked(IsPosIntSym(nn), IsIntSym(rot), WithBases(self,
        rec(n := nn, rot := rot, operations := PrintOps)))),

    hashAs := self >> ObjId(self)(self.n, 1),
    rChildren := self >> [self.n, self.rot],
    rSetChild := rSetChildFields("n", "rot"),
    print := self >> Print(self.__name__, "(", self.n, ", ", self.rot,")"),

    conj := self >> CopyFields(self, rec(rot := -self.rot)),

    prettyPrint := self >> PrintEval("x$1 - 1", When(self.n=1, "", Concat("^", StringInt(self.n)))),

    rspectrum := self >> let(rot := EvalScalar(self.rot), n := EvalScalar(self.n), Cond(
        n <= 1, [self], # if n<=1 then then can't decompose further
        Concatenation(
            [ XN_min_1(1, rot) ],
            Cond(IsOddInt(n), [], [XN_plus_1(1, rot)]),
            List([1..Int((n-1)/2)], i -> XN_skew(2, i/n, rot)))
    )),

    cspectrum := self >> let(rot := EvalScalar(self.rot), n := EvalScalar(self.n), Cond(
        n = 1, [self],
        n = 2, [XN_min_1(1,self.rot), XN_plus_1(1, self.rot)],
        n > 2, ConcatList(self.rspectrum(), x -> x.cspectrum())
    )),

    # rfactor(1) == rspectrum()
    # rfactor(n) == self
    rfactor := (self, m) >> Checked((self.n mod m) = 0,
        List(XN_min_1(self.n/m, self.rot).rspectrum(), x -> CopyFields(x, rec(n := x.n * m)))),

    shift := self >> Z(self.n, -1)
));

#F XN_min_1U(n) -- represents x^n - 1 polynomial
#F   Spectral decomposition: 2d components only
#F   <n> must be even
#F
Class(XN_min_1U, XN_min_1, rec(
    __call__ := (self, n, rot) >> let(nn:=_unwrap(n), Checked(IsPosIntSym(nn), IsSymbolic(nn) or IsEvenInt(nn),
        IsIntSym(rot), WithBases(self,
            rec(n := nn, rot := rot, operations := PrintOps)))),

    rspectrum := self >> let(rot := EvalScalar(self.rot), n := EvalScalar(self.n), Cond(
        self.n <= 1, [self], # if n==1 then then can't decompose further
        Concatenation(
            Cond(IsOddInt(n), [XN_min_1(1, rot)], [XN_min_1(2, rot)]),
            List([1..Int((n-1)/2)], i -> XN_skew(2, i/n, rot)))
    ))
));

Class(XN_skew_base, ASPAlgebra, rec(
    prettyPrint := self >> Cond(IsOddInt(self.n) or self.a=1/4,
        PrintEval("x$1 + 1", When(self.n=1, "", Concat("^", StringInt(self.n)))),
        PrintEval("x^$1 - 2 x$2 cos($3*pi) + 1", self.n, When(self.n=2, "", Concat("^", StringInt(self.n/2))), 2*self.a)
    ),

    conj := self >> CopyFields(self, rec(rot := -self.rot)),

    shift := self >> let(n:=self.n, cos:=CosPi(self.rot*2*self.a), i := Ind(n),
        HStack(VStack(O(1,n-1),I(n-1)),
            Mat([List([0..n-1], i -> Cond(i=0, -1, i=n/2, 2*cos, 0))]).transpose())),

    cspectrum := self >> let(rot := EvalScalar(self.rot), n := EvalScalar(self.n), Cond(
        n = 1, [self],
        n = 2, [XN_min_w(1,self.a, rot), XN_min_w(1,1-self.a, rot)],
        n > 2, ConcatList(self.rspectrum(), x -> x.cspectrum())
    ))
));

#F XN_skew(n, a, rot) -- represents x^n - 2x^{n/2} cospi(2*a) + 1 polynomial
#F   Spectral decomposition: 2d components only
#F   <n> must be even
#F
Class(XN_skew, XN_skew_base, rec(
    __call__ := (self, n, a, rot) >> let(nn := _unwrap(n),
	Checked(IsPosIntSym(nn), IsSymbolic(nn) or IsEvenInt(nn), IsIntSym(rot),
	        IsRatSym(a), IsSymbolic(a) or (0 < a and a < 1),
		WithBases(self,
		    rec(n := nn, a := a, rot := rot, operations := PrintOps)))),

    hashAs := self >> XN_skew(self.n, 1/16, 1),
    rChildren := self >> [self.n, self.a, self.rot],
    rSetChild := rSetChildFieldsF(_unwrap, "n", "a", "rot"),
    print := self >> Print(self.__name__, "(", self.n, ", ", self.a, ", ", self.rot, ")"),

    rspectrum := self >> let(rot := EvalScalar(self.rot), n := EvalScalar(self.n), Cond(
        n <= 2, [self], # if n<=2 then then can't decompose further
        List([0..Int(n/2)-1], i -> XN_skew(2, (i + self.a)/(n/2), rot)))),

    # rfactor(1) == rspectrum(), rfactor(n) == self
    rfactor := (self, m) >> Checked(IsEvenInt(m), m > 1, (self.n mod m) = 0,
        List(XN_skew(2*self.n/m, self.a, self.rot).rspectrum(), x -> CopyFields(x, rec(n := x.n * m/2)))),
));

#F XN_plus_1(n, rot) -- represents x^n + 1 polynomial
#F   Spectral decomposition: 1d and 2d components
#F   <n> even or odd
#F
Class(XN_plus_1, XN_skew_base, rec(
    __call__ := (self, n, rot) >> let(nn:=_unwrap(n), Checked(IsPosIntSym(nn), IsIntSym(rot),
        WithBases(self, rec(n := nn, a := 1/4, rot := rot, operations := PrintOps)))),

    rspectrum := self >> let(rot := EvalScalar(self.rot), n := EvalScalar(self.n), Cond(
        n <= 2, [self], # if n<=2 then then can't decompose further
        Concatenation(
            List([0..Int(n/2)-1], i -> XN_skew(2, (i + 1/2)/(n), rot)),
            Cond(IsEvenInt(n), [], [XN_plus_1(1, rot)]))
    )),

    rfactor := (self, mm) >> Checked(mm >= 1, (self.n mod mm) = 0, Cond(
        IsOddInt(mm) and not IsOddInt(self.n), Error("<mm> can only be odd if polynomial degree (self.n) is odd"),
        let(m := When(IsOddInt(self.n), mm, mm/2),
            List(XN_plus_1(self.n/m, self.rot).rspectrum(), x -> CopyFields(x, rec(n := x.n * m)))))),

    hashAs := self >> XN_plus_1(self.n, 1),
    rChildren := self >> [self.n, self.rot],
    rSetChild := rSetChildFields("n", "rot"),
    print := self >> Print(self.__name__, "(", self.n, ", ", self.rot, ")"),
));

#F XN_min_w(n, a, rot) -- represents x^n - w_a polynomial
#F   Spectral decomposition: 1d components only over complex numbers
#F
Class(XN_min_w, XN_skew_base, rec(
    __call__ := (self, n, a, rot) >> let(nn:=_unwrap(n), Checked(
	IsPosIntSym(nn), IsRatSym(a), IsSymbolic(a) or (0 < a and a < 1), IsIntSym(rot),
        WithBases(self, rec(n := nn, a := a, rot := rot, operations := PrintOps)))),

    hashAs := self >> XN_min_w(self.n, 1/16, 1),
    rChildren := self >> [self.n, self.a, self.rot],
    rSetChild := rSetChildFields("n", "a", "rot"),
    print := self >> Print(self.__name__, "(", self.n, ", ", self.a, ", ", self.rot, ")"),

    rspectrum := self >> Error("Polynomial algebra is complex, no real spectrum"),

    cspectrum := self >> let(rot := EvalScalar(self.rot), n := EvalScalar(self.n), Cond(
        n = 1, [self], # if n<=2 then then can't decompose further
        List([0..n-1], i -> XN_min_1(1, (i + self.a)/n, rot)))),

    prettyPrint := self >> PrintEval("x^$1 - w_$2", self.n, self.a),

    shift := self >> let(n:=self.n, cos:=CosPi(self.rot*2*self.a), i := Ind(n),
        DirectSum(_omega(EvalScalar(self.a))*I(1), When(n=1, [], I(n-1)))*
        Prm(Z(n,-1)))
));


text:
NrQueens := function(n)
    local a, up, down, m, sub;
    a := [1 .. n];
    up := ListWithIdenticalEntries(2*n - 1, true);
    down := ListWithIdenticalEntries(2*n - 1, true);
    m := 0;
    sub := function(i)
        local j, k, p, q;
        for k in [i .. n] do
            j := a[k];
            p := i + j - 1;
            q := i - j + n;
            if up[p] and down[q] then
                if i = n then
                    m := m + 1;
                else
                    up[p] := false;
                    down[q] := false;
                    a[k] := a[i];
                    a[i] := j;
                    sub(i + 1);
                    up[p] := true;
                    down[q] := true;
                    a[i] := a[k];
                    a[k] := j;
                fi;
            fi;
        od;
    end;
    sub(1);
    return m;
end;

Queens := function(n)
    local a, up, down, v, sub;
    a := [1 .. n];
    up := ListWithIdenticalEntries(2*n - 1, true);
    down := ListWithIdenticalEntries(2*n - 1, true);
    v := [];
    sub := function(i)
        local j, k, p, q;
        for k in [i .. n] do
            j := a[k];
            p := i + j - 1;
            q := i - j + n;
            if up[p] and down[q] then
                if i = n then
                    Add(v, ShallowCopy(a));
                else
                    up[p] := false;
                    down[q] := false;
                    a[k] := a[i];
                    a[i] := j;
                    sub(i + 1);
                    up[p] := true;
                    down[q] := true;
                    a[i] := a[k];
                    a[k] := j;
                fi;
            fi;
        od;
    end;
    sub(1);
    return v;
end;

NrQueens(8);
a := Queens(8);;
PrintArray(PermutationMat(PermList(a[1]), 8));

[ [  1,  0,  0,  0,  0,  0,  0,  0 ],
  [  0,  0,  0,  0,  1,  0,  0,  0 ],
  [  0,  0,  0,  0,  0,  0,  0,  1 ],
  [  0,  0,  0,  0,  0,  1,  0,  0 ],
  [  0,  0,  1,  0,  0,  0,  0,  0 ],
  [  0,  0,  0,  0,  0,  0,  1,  0 ],
  [  0,  1,  0,  0,  0,  0,  0,  0 ],
  [  0,  0,  0,  1,  0,  0,  0,  0 ] ]

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
# return a list of two lists :
# first is the list of months with five weekends between years y1 and y2 (included)
# second is the list of years without such months, in the same interval
FiveWeekends := function(y1, y2)
  local L, yL, badL, d, m, y;
  L := [ ];
  badL := [ ];
  for y in [y1 .. y2] do
    yL := [ ];
    for m in [1, 3, 5, 7, 8, 10, 12] do
      if WeekDay([1, m, y]) = "Fri" then
        d := StringDate([1, m, y]);
        Add(yL, d{[4 .. 11]});
      fi;
    od;
    if Length(yL) = 0 then
      Add(badL, y);
    else
      Append(L, yL);
    fi;
  od;
  return [ L, badL ];
end;

r := FiveWeekends(1900, 2100);;
n := Length(r[1]);
# 201
Length(r[2]);
# 29
r[1]{[1 .. 5]};
# [ "Mar-1901", "Aug-1902", "May-1903", "Jan-1904", "Jul-1904" ]
r[1]{[n-4 .. n]};
# [ "Mar-2097", "Aug-2098", "May-2099", "Jan-2100", "Oct-2100" ]

max_stars_repo_name: mullikine/RosettaCodeData


text:
################################################################################
##
#W verlinde.gi              GroupTheoretical Package
##
#W Paul Bruillard, Cesar Galindo, Siu-Hung Ng, Julia Plavnik, Eric Rowell,
#W Zhenghan Wang
##
## Installation file for verlinde functions of the GroupTheoretical Package
##
#Y Copyright (C) 2016, Battelle Memorial Institute
##
################################################################################

################################################################################
##
#F compute_fusion_rules(<list>,<number>). . . . . . . .use verlinde to compute N
##
InstallGlobalFunction(compute_fusion_rules, function(S,FPdimC)
  local N,r,i,j,k,m;
  r:=Length(S[1]);
  N:=List([1..r],i->NullMat(r,r));
  for i in [1..r] do
    for j in [1..r] do
      for k in [1..r] do
        for m in [1..r] do
          N[i][j][k]:=N[i][j][k]+S[i][m]*S[j][m]*ComplexConjugate(S[k][m])/S[1][m];
        od;
        N[i][j][k]:=N[i][j][k]/FPdimC;
      od;
    od;
  od;
  return N;
end);

################################################################################
##
#F display_fusion_rules(<list>). . . . . . . . . display the fusion rules
##
InstallGlobalFunction(display_fusion_rules, function(N)
  local x;
  for x in [1..Size(N)] do
    Print("\nN_",x,"\n");
    Display(N[x]);
  od;
  Print("\n");
end);

#E verlinde.gi . . . . . . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: pnnl/GroupTheoretical

max_issues_repo_name: pnnl/GroupTheoretical

max_forks_repo_name: pnnl/GroupTheoretical


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

Class(ScratchpadGlobals, rec(
    getOpts := meth(arg)
        local lssize, opts, swp, brules, nrules, br, nrsgmts, vlen, globalUnrolling,size, ttype;

	lssize := When (Length(arg) >= 2, arg[2], 2);
	nrsgmts := When (Length(arg) >= 3, arg[3], 1);
	vlen := When (Length(arg) >= 4, arg[4], 1);
	size := When (Length(arg) >= 5, arg[5], 2);
    ttype := When (Length(arg) >= 6, arg[6], 'R');
	swp := When (Length(arg) >= 7, arg[7], false);
	globalUnrolling := When(Length(arg) >=8, arg[8], 1);

    brules := When(IsRec(SpiralDefaults.breakdownRules),
        UserRecFields(SpiralDefaults.breakdownRules),
        Filtered(Dir(SpiralDefaults.breakdownRules), i->not i in SystemRecFields));
    nrules := rec();
    for br in brules do
        nrules.(br) := List(SpiralDefaults.breakdownRules.(br), i->CopyFields(i));
    od;
    opts := CopyFields(SpiralDefaults);
    opts.breakdownRules := nrules;

    opts.breakdownRules.TCompose := [ TCompose_tag];
    opts.breakdownRules.DFT := [DFT_Base, DFT_CT, DFT_tSPL_CT, DFT_PD, DFT_Rader];
    opts.breakdownRules.TTwiddle := [ TTwiddle_Tw1];

	opts.breakdownRules.WHT := [WHT_tSPL_BinSplit, WHT_Base, WHT_BinSplit];

	opts.breakdownRules.TTensor := [AxI_IxB,IxB_AxI];
	opts.breakdownRules.TTensorI := Concat([IxA_scratch_push, IxA_base, AxI_base, IxA_L_base, L_IxA_base], [IxA_scratch, AxI_scratch, IxAL_scratch]);

	opts.tags := [Cond( ttype = 'R', ALStore(lssize,nrsgmts,vlen), ALStoreCx(lssize,nrsgmts,vlen)) ];

	opts.formulaStrategies.sigmaSpl := [ MergedRuleSet(RulesSumsScratch, RulesFuncSimpScratch, RulesDiag, RulesDiagStandalone, RulesStrengthReduce, RulesRCScratch,RulesII,OLRules) ];
    opts.formulaStrategies.preRC := [ MergedRuleSet(RulesSumsScratch, RulesFuncSimpScratch, RulesDiag, RulesDiagStandalone, RulesStrengthReduce, RulesRCScratch, RulesII, OLRules), (s,o) -> ScratchModel.updateInfo(s) ];
	opts.formulaStrategies.rc := [ MergedRuleSet(RulesSums, RulesFuncSimp, RulesDiag, RulesDiagStandalone, RulesStrengthReduce, RulesRCScratch, RulesII, OLRules) ];
    #opts.formulaStrategies.postProcess := [(s, opts) -> compiler.BlockSums(opts.globalUnrolling, s)];
    opts.size := size;
	opts.swp := swp;
    opts.globalUnrolling := globalUnrolling;

    opts.sumsgen := ScratchSumsGen;
	opts.codegen := ScratchCodegen;
    opts.unparser := CScratchUnparserProg;

    opts.memModifier := "__memory";
    opts.scratchModifier := "__scratch";
    opts.arrayDataModifier := "__rom";
    opts.romModifier := "__rom";

	opts.includes := [];
    Add(opts.includes, "\"scratch.h\"");

    opts.dmaSignal := (self, opts) >> "DMA_signal";
    opts.dmaWait := (self, opts) >> "DMA_wait";
    opts.cpuSignal := (self, opts) >> "CPU_signal";
    opts.cpuWait := (self, opts) >> "CPU_wait";
    opts.dmaFence := (self, opts) >> "DMA_fence";
    opts.dmaLoad := (self, opts) >> "DMA_load";
    opts.dmaStore := (self, opts) >> "DMA_store";
    opts.model := ScratchModel;

    return opts;
    end
));


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

# -----------------------------
# Special extension definitions
# ------------------------------

# periodic
Class(per, rec(
    __call__ := (self, k) >> WithBases(self, rec(n:=k, operations:=PrintOps)),
    print := self >> Print(self.name, "(", self.n, ")"),
    left  := (self, N) >> fAdd(N, self.n, N - self.n),
    right := (self, N) >> fAdd(N, self.n, 0),
));

# half-point symmetric
Class(hs, rec(
    __call__ := (self, k) >> WithBases(self, rec(n:=k, operations:=PrintOps)),
    print := self >> Print(self.name, "(", self.n, ")"),

    left  := (self, N) >> fCompose(fAdd(N, self.n, 0), J(self.n)),
    right := (self, N) >> fCompose(fAdd(N, self.n, N-self.n),J(self.n))
));

# whole-point symmetric
Class(ws, rec(
    __call__ := (self, k) >> WithBases(self, rec(n:=k, operations:=PrintOps)),
    print := self >> Print(self.name, "(", self.n, ")"),

    left  := (self, N) >> fCompose(fAdd(N, self.n, 1), J(self.n)),
    right := (self, N) >> fCompose(fAdd(N, self.n, N-self.n-1),J(self.n))
));

# zero padding - doing nothing, just returning the extension length
Class(zer, rec(
    __call__ := (self, k) >> WithBases(self, rec(n:=k, operations:=PrintOps)),
    print := self >> Print(self.name, "(", self.n, ")"),
    left  := (self, N) >> self.n,
    right := (self, N) >> self.n,
));

# ******************************
#  Operations on polynomials
# ******************************

ListPoly := (p) -> [p.coefficients, p.valuation];

ListSum := function(l1,l2)
  local i, d1, d2, dm, l;

  l := [];
  d1 := Length(l1);
  d2 := Length(l2);
  dm := Minimum(d1,d2);
  for i in [1..dm] do
     l[i] := l1[i]+l2[i];od;
  if d2 > dm then Append(l,l2{[dm+1..d2]});
  else Append(l,l1{[dm+1..d1]});
  fi;
  return l;
end;

ListPolyMat := function(M)
   local S,L,V;
   S:=List(M, j-> ListPoly(j));
   L:=List(S, i->i[1]);
   V:=List(S, i->i[2]);

  return [L,V];
end;

Shift := function(h, sh)
local g;

if IsPolynomial(h) then
g := ShallowCopy(h);
g.valuation := h.valuation + sh;
return g;
else
return [h[1],h[2]+sh];
fi;
end;

Upsample := function(h)
local i, L, g;
L:=[];
 g:=ShallowCopy(h);
 if Length(h.coefficients)=1 then
  Add(L,h.coefficients[1]);
 elif Length(h.coefficients)>1 then
  for i in [1..Length(h.coefficients)-1] do
    Add(L, h.coefficients[i]);
    Add(L, h.baseRing.zero);
  od;
  Add(L, h.coefficients[i+1]);
 fi;
  g.coefficients := L;
  g.valuation := h.valuation * 2;
  return g;
end;

# Downsample creates apolynomial of even coefficients
# h_e(z^2) =1/2*( h(z)+h(-z) )

Downsample := function(h)
local i, L, g, even, odd, v, l, k, coef;

 if IsPolynomial(h) then L:= ListPoly(h); coef:=L[1]; l:=L[2];
 else coef:=h[1]; l:=h[2]; fi;
 L:=[];
 k := l mod 2;
    even := Filtered([1..Length(coef)], i-> i mod 2 = 0);
    odd := Filtered([1..Length(coef)], i-> i mod 2 = 1);

  if (k=0) then
    L:= Sublist(coef,odd);
    v := l/2;
  else
    L := Sublist(coef,even);
    v := Int((l+1)/2);
  fi;

  #g.coefficients := L;
  #g.valuation :=v;
  return [L,v];
end;

DownsampleTwo := function(h)
local he, ho;

if IsPolynomial(h) then h:=ListPoly(h); fi;
he:=Downsample(h);
ho:=Downsample(Shift(h,1));
if h[2] mod 2 = 0 then
ho:=Shift(ho,-1);
fi;

return[he,ho];
end;

DownsampleTwoNeg := function(h)
local he, ho;

he:=Downsample(h);
ho:=Downsample(Shift(h,-1));
if h[2] mod 2 = 0 then
ho:=Shift(ho,1);
fi;

return[he,ho];
end;

PosInt := (int) -> (int + AbsInt(int))/2;

#F Positive(n)
#F returns positive value of n thresholded at 0
#F
Positive := (n) -> When(n>0, n, 0);

#F Poly(L,v)
#F Creates a polynomial L[1]*z^(v+k-1) + ... + L[k]*z^v
#F

Poly := function(arg)
local z,L,v;
z:=Indeterminate(Cyclotomics);
z.name:="z";
L := When(IsList(arg[1]),arg[1],[arg[1]]);
v := When(Length(arg)<2,0,arg[2]);
return Polynomial(Cyclotomics, L, v);
end;

#F PolyMat(L,V)
#F Creates a polynomial matrix out of coefficient and valuation matrices
#F

PolyMat := function(L,V)
local M,dim;
dim:=Dimensions(V);
M := List([1..dim[1]],i->List([1..dim[2]],j->Poly(L[i][j],V[i][j])));
return M;
end;

#F PolyExtension(p(z))
#F Computes extension lengths for poynomial p(z)
#F
PolyExtension := p -> let(
    l := Positive(-p.valuation),
    r := Positive(LaurentDegree(p)+p.valuation),
    [l, r]
);

#F FillZeros := function(L,v)
#F fills zeros in the list of polynomial coefficients up to z^0
#F
FillZeros := function(p)
local w,Ls,L,v;
if IsList(p) then
L:=p[1];
v:=p[2];
else
L:= p.coefficients;
v:= p.valuation;
fi;

w:= v+Length(L)-1;
Ls:= L;
if w<=0 then
  Ls := Concat(L,List([1..-w], i->0));
  #k := -w;
elif v>=0 then
  Ls := Concat(List([1..v], i->0), L);
  #k := v;
fi;
return [Ls,-Positive(-v)];
end;

#F FillZerosMat := function(L,v)
#F fills zeros in the list of polynomial coefficients up to z^0
#F

FillZerosMat := function(p)
local S,L,V,dim;

dim := Dimensions(p);
L:=List(p, i->List(i, j->  j.coefficients));
V:=List(p, i->List(i, j->  j.valuation));

S:=List([1..dim[1]],i->List([1..dim[2]],j-> FillZeros([L[i][j],V[i][j]])));
L:=List(S, i->List(i, j-> j[1]));
V:=List(S, i->List(i, j-> j[2]));

return [L,V];
end;

TimeReverse := (L,v) -> let(
               deg := Length(L),
               [Reversed(L),-(v+deg)+1]);

TimeReverseMat := function(arg)
   local S,L,V,Ls,Vs,deg;

if Length(arg)=1 then L:=arg[1][1]; V:=arg[1][2];
else L:=arg[1]; V:=arg[2]; fi;

deg := Length(V);
S:=List([1..deg],i->TimeReverse(L[i],V[i]));
Ls:=List(S,i->i[1]);
Vs:=List(S,i->i[2]);
return [Ls,Vs];
end;

Alternate := (L) -> [List([1..Length(L[1])], i->(-1)^(L[2]+i-1)*L[1][i]),L[2]];

Synth2AnalFilts := function(L,V)

 local ht,gt,h,g;

 h:=TimeReverse(L[1],V[1]);
 g:=TimeReverse(L[2],V[2]);

 h:=Alternate(h);
 g:=Alternate(g);

 ht:=[];
 gt:=[];

 ht[1] := g[1];
 ht[2] := g[2]+1;
 gt[1] := -h[1];
 gt[2] := h[2]+1;

 return [[ht[1],gt[1]],[ht[2],gt[2]]];

end;

Anal2SynthFilts := function(L,V)

 local ht,gt,g,h;
 h:=TimeReverse(L[1],V[1]);
 g:=TimeReverse(L[2],V[2]);

 h:=Alternate(h);
 g:=Alternate(g);

 ht:=[];
 gt:=[];

 ht[1] := g[1];
 ht[2] := g[2]+1;
 gt[1] := -h[1];
 gt[2] := h[2]+1;

 return [[ht[1],gt[1]],[ht[2],gt[2]]];

end;

# wraps a polynomial given as Poly or as List+valuation periodically
# for circulants
CircularWrap := function(arg)
local L,v,l,r,n;
if (Length(arg) = 2) then
if IsPolynomial(arg[1]) then
l:=ListPoly(arg[1]);
else l:=arg[1];fi;
L:=l[1];
v:=l[2];
n:=arg[2];
else
L:=arg[1];
v:=arg[2];
n:=arg[3];
fi;

r:= Length(L)+v-1;
if v>0 then L:=Concat(Replicate(v,0),L);fi;

if r <n then
L:=Concat(L, Replicate(n-r-1,0));
L:=L*MatSPL(Ext(n,per(PosInt(-v)),per(0)));
else
L:=L*MatSPL(Ext(n,per(PosInt(-v)),per(r-n+1)));
fi;

# Circulant needs first column as param instead of first row (need to be changed)

l:= Reversed(L{[2..Length(L)]});
L:=Concat([L[1]],l);

return L;
end;

#
# Helper functions for non-terminal arguments
#
toSize := o -> Cond(
    IsPosInt(o), o,
    IsList(o), Length(o),
    IsPolynomial(o), Length(o.coefficients),
    IsFunction(o), o.domain(),
    Error("List, polynomial, or generating function expected"));

toFunc := o -> Cond(
    IsPosInt(o), FUnk(o),
    IsList(o), Checked(o<>[], FList(UnifyTypesV(o), o)),
    IsPolynomial(o), FList(UnifyTypesV(o.coefficients), o.coefficients),
    IsFunction(o),   o,
    Error("List, polynomial, or generating function expected"));

toValuation := o -> Cond(
    IsPosInt(o), 0,
    IsList(o), 0,
    IsPolynomial(o), o.valuation,
    IsFunction(o), 0,
    Error("List, polynomial, or generating function expected"));

QuantizeQuarters := rat -> Cond(0  <= rat and rat <= 3/8, 1/4,
                                3/8 < rat and rat <= 5/8, 1/2,
                                5/8 < rat and rat <= 7/8, 3/4,
				1);

toFiltFunc := o -> Cond(
    IsPosInt(o),
        [ fUnk(TReal, o), -o+1],
    IsList(o), Checked(o<>[],
	[ toFunc(o), -Length(o)+1 ]),
    IsPolynomial(o), let(zp := FillZeros(o),
        [ toFunc(zp[1]), zp[2] ]),
    IsFunction(o) or IsFuncExp(o),
        [ o, -o.domain()+1 ],
    Error("List, polynomial, or generating function expected"));

toFiltFuncV := (o,v) -> Cond(
    IsPosInt(o), let(w := v+o-1,
	When(w < 0 or v > 0,
	    Error("Valuation does not match domain of <o>, can't zero-pad a gener. function"),
        [ fUnk(TReal, o), v])),
    IsList(o), Checked(o<>[], let(zp := FillZeros([o,v]),
	[ toFunc(zp[1]), zp[2] ])),
    IsFunction(o) or IsFuncExp(o), let(dom := o.domain(), w := v+dom-1,
	When(w < 0 or v > 0,
	    Error("Valuation does not match domain of <o>, can't zero-pad a gener. function"),
	    [ o, -o.domain()+1 ])),
    Error("Only list or generating function can be specified with separate valuation"));


text:
# All arrangements of 4 elements in 1 .. 4
Arrangements([1 .. 4], 4);
# All permutations of 1 .. 4
PermutationsList([1 .. 4]);


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

# VGath_sv/VScat_sv may have bigger dimensions than their function domain so we need padding
# to have Diag with correct dimensions when commuting diag through VGath_sv and VScat_sv
_fComposeVGathScatSv_Diag := (m, dim, d) -> let(
    g := fTensor(m.func, fId(m.sv)),
    mdim := dim(m),
    rem := Cond(m.rem _is Unk, g.domain() mod mdim, m.rem=0, 0, m.v - m.rem * m.sv),
    f := fCompose(d.element, g),
    When(rem = 0, f, diagDirsum(f, fConst(TReal, mdim-rem, 0.0)))
);

Class(RulesVDiag, RuleSet);
RulesVDiag.__transparent__ := [ [Compose, [vRC]] ];  # creates extra rules, A*B->C is made into vRC(A)*vRC(B) -> vRC(C)

RewriteRules(RulesVDiag, rec(
    DiagConst_AnyScat := ARule(Compose, [[Diag, [fConst, @(1), @(2), @(3)]], @(4).cond(IsVScat)],
	e -> [ @(4).val, VDiag(fConst(@(1).val, Cols(@(4).val), @(3).val), getV(@(4).val))]),

    AnyGath_DiagConst := ARule(Compose, [@(4).cond(IsVGath), [Diag, [fConst, @(1), @(2), @(3)]]],
	e -> [ VDiag(fConst(@(1).val, Rows(@(4).val), @(3).val), getV(@(4).val)), @(4).val]),

    DiagConst_NoDiagPullin := ARule(Compose, [[Diag, [fConst, @(1), @(2), @(3)]], @(4, NoDiagPullin)],
	e -> [ @(4).val, Diag(fConst(@(1).val, Cols(@(4).val), @(3).val))]),

    # VGath * IDirSum
     CommuteVGathIDirSum := ARule( Compose,
	 [ [@(1, VGath), @(3)],
	   [@(0, [PushL, PushLR]), @(4, IDirSum, e -> (@(1).val.v mod Rows(e.child(1))) = 0)] ],
       e -> let(
             ratio := @(1).val.v / Rows(@(4).val.child(1)),
             f := When(ratio=1, @(3).val, fTensor(@(3).val, fId(ratio))),
             j := @(4).val.var,
             jj := Ind(f.domain()),

             [ ObjId(@(0).val)(IDirSum(jj,
                   Data(j, f.at(jj), SubstTopDown(@(4).val.child(1), @(5,fBase,e->e.params[2]=j),
                                       e -> fCompose(f, fBase(jj)))))),
               @(1).val ])
    ),

    # IDirSum * VScat
    CommuteIDirSumVScat := ARule( Compose,
       [ [@(0,[PushR, PushLR]), @(4, IDirSum)],
         [@(1, VScat, e -> (e.v mod Cols(@(4).val.child(1)))=0), @(3)] ],

       e -> let(
             ratio := @(1).val.v / Cols(@(4).val.child(1)),
             f := When(ratio=1, @(3).val, fTensor(@(3).val, fId(ratio))),
             j := @(4).val.var,
             jj := Ind(f.domain()),

             [ @(1).val,
               ObjId(@(0).val)(IDirSum(jj,
                   Data(j, f.at(jj), SubstTopDown(@(4).val.child(1), @(5,fBase,e->e.params[2]=j),
                                       e -> fCompose(f, fBase(jj)))))) ])
    ),

    # Gath * RCDiag
    CommuteVGathVRCDiag_fTensor := ARule( Compose,
       [ [@(1, VGath), [@(0,fTensor), ..., [fId, @(2).cond(_isEven)]]], @(4, VRCDiag) ],
       e -> [ VRCDiag(fCompose(@(4).val.element, @(0).val), @(4).val.v), @(1).val ]),

    # Gath * RCDiag
    CommuteVGathVRCDiag_H := ARule( Compose,
       [ [@(1, VGath), [@(2, H), @, @.cond(_isEven), @.cond(_isEven), _1]], @(4, VRCDiag) ],
       e -> [ VRCDiag(fCompose(@(4).val.element, @(2).val), @(4).val.v), @(1).val ]),

    # RCVDiag * VScat
    CommuteRCVDiagVScat := ARule( Compose,
       [ @(1, RCVDiag), @(2, VScat) ], # <-1 <-2 o
        e -> [ @2.val, RCVDiag(fCompose(@1.val.element, @2.val.func), @2.val.v).attrs(@(1).val) ]),

    # RCVDiagSplit * VScat
    CommuteRCVDiagSplitVScat := ARule( Compose,
       [ @(1, RCVDiagSplit), @(2, VScat) ], # <-1 <-2 o
        e -> [ @2.val, RCVDiagSplit(fCompose(@1.val.element, fTensor(@2.val.func, fId(2))), @2.val.v).attrs(@(1).val) ]),

    # VRCDiag * VScat
    CommuteVRCDiagVScat_fTensor := ARule( Compose,
        [ @(4, VRCDiag), [@(1, VScat), [@(0,fTensor), ..., [fId, @(2).cond(_isEven)]]] ],
        e -> [ @(1).val, VRCDiag(fCompose(@(4).val.element, @(0).val), @(4).val.v) ]),

    # VRCDiag * VScat
    CommuteVRCDiagVScat_H := ARule( Compose,
        [ @(4, VRCDiag), [@(1, VScat), [@(2, H), @, @.cond(_isEven), @.cond(_isEven), _1]] ],
        e -> [ @(1).val, VRCDiag(fCompose(@(4).val.element, @(2).val), @(4).val.v) ]),

    # VGath * Diag
    CommuteVGathDiag := ARule( Compose,
       [ @(1, VGath), @(2, Diag) ], # o 1-> 2->
        e -> [ VDiag(fCompose(@2.val.element, fTensor(@1.val.func, fId(@1.val.v))), @1.val.v).attrs(@(2).val), @1.val ]),

    # VGath * VDiag
    CommuteVGathVDiag := ARule( Compose,
       [ @(1, VGath), @(2, VDiag) ], # o 1-> 2->
        e -> [ VDiag(fCompose(@2.val.element, fTensor(@1.val.func, fId(@1.val.v))), @1.val.v).attrs(@(2).val), @1.val ]),

    # Diag * VScat
    CommuteDiagVScat := ARule( Compose,
       [ @(1, Diag), @(2, VScat) ], # <-1 <-2 o
        e -> [ @2.val, VDiag(fCompose(@1.val.element, fTensor(@2.val.func, fId(@2.val.v))), @2.val.v).attrs(@(1).val) ]),

    # VDiag * VScat
    CommuteVDiagVScat := ARule( Compose,
       [ @(1, VDiag), @(2, VScat) ], # <-1 <-2 o
        e -> [ @2.val, VDiag(fCompose(@1.val.element, fTensor(@2.val.func, fId(@2.val.v))), @2.val.v).attrs(@(1).val) ]),

    # RCVDiag * VScat
    CommuteRCVDiagVScat := ARule( Compose,
       [ @(1, RCVDiag), @(2, VScat) ], # <-1 <-2 o
        e -> [ @2.val, RCVDiag(fCompose(@1.val.element, @2.val.func), @2.val.v).attrs(@(1).val) ]),

    # RCVDiagSplit * VScat
    CommuteRCVDiagSplitVScat := ARule( Compose,
       [ @(1, RCVDiagSplit), @(2, VScat) ], # <-1 <-2 o
        e -> [ @2.val, RCVDiagSplit(fCompose(@1.val.element, fTensor(@2.val.func, fId(2))), @2.val.v).attrs(@(1).val) ]),

    # VGath * VDiag_x_I
    CommuteVGathVDiag_x_I := ARule( Compose,
        [ @(1, VGath), @(2, VDiag_x_I) ], # o 1-> 2->
    e -> [ VDiag_x_I(fCompose(@2.val.element, @1.val.func), @2.val.v).attrs(@(2).val), @1.val ]),

    # VDiag_x_I * VScat
    CommuteVDiag_x_IVScat := ARule( Compose,
        [ @(1, VDiag_x_I), @(2, VScat) ], # <-1 <-2 o
    e -> [ @2.val, VDiag_x_I(fCompose(@1.val.element, @2.val.func), @1.val.v).attrs(@(1).val) ]),

    #------------------------
    # VGath_sv * Diag
    CommuteVGathsvDiag := ARule( Compose,
       [ @(1, VGath_sv), @(2, Diag) ], # o 1-> 2->
       e -> [ VDiag(_fComposeVGathScatSv_Diag(@1.val, Rows, @2.val), @1.val.v).attrs(@2.val), @1.val ]
    ),

    # VGath_sv * VDiag
    CommuteVGathsvVDiag := ARule( Compose,
       [ @(1, VGath_sv), @(2, VDiag) ], # o 1-> 2->
       e -> [ VDiag(_fComposeVGathScatSv_Diag(@1.val, Rows, @2.val), @1.val.v).attrs(@2.val), @1.val ]
    ),

    # Diag * VScat_sv
    CommuteDiagVScatsv := ARule( Compose,
       [ @(1, Diag), @(2, VScat_sv) ], # <-1 <-2 o
       e -> [ @2.val, VDiag(_fComposeVGathScatSv_Diag(@2.val, Cols, @1.val), @2.val.v).attrs(@1.val) ]
    ),

    # VDiag * VScat_sv
    CommuteVDiagVScatsv := ARule( Compose,
       [ @(1, VDiag), @(2, VScat_sv) ], # <-1 <-2 o
       e -> [ @2.val, VDiag(_fComposeVGathScatSv_Diag(@2.val, Cols, @1.val), @2.val.v).attrs(@1.val) ]
    ),

    #-----------------------------------------------------------
    # VStretchGath * Diag
    CommuteVStretchGathDiag := ARule( Compose,
       [ @(1, VStretchGath), @(2, Diag) ], # o 1-> 2->
       e -> let( vg := @1.val, den := vg.func.domain()/vg.part, num := _roundup(den, vg.v),
            [ VDiag(fStretch(fCompose(@2.val.element, @1.val.func), num, den), @1.val.v).attrs(@2.val), @1.val ])
    ),

    # Diag * VStretchScat
    CommuteVStretchScatDiag := ARule( Compose,
       [ @(1, Diag), @(2, VStretchGath)],
       e -> let( vg := @2.val, den := vg.func.domain()/vg.part, num := _roundup(den, vg.v),
            [ @2.val, VDiag(fStretch(fCompose(@1.val.element, @2.val.func), num, den), @2.val.v).attrs(@1.val) ])
    ),

    #-----------------------------------------------------------
    # Diag * IxVScat_pc
    CommuteDiagIxVScat_pc := ARule( Compose,
	[ @(1, [Diag, VDiag]), @(2, IxVScat_pc) ], # <-1 <-2 o
        e -> let(s:=@(2).val, [ s, VDiag(fStretch(@(1).val.element, _roundup(s.n, s.v), s.n), s.v) ])),

    # IxVGath_pc * Diag
    CommuteIxVGath_pcDiag := ARule( Compose,
	[ @(1, IxVGath_pc), @(2, [Diag, VDiag]) ], # o 1-> 2->
        e -> let(g:=@(1).val, [ VDiag(fStretch(@(2).val.element, _roundup(g.n, g.v), g.n), g.v), g ])),

    # Diag * VScat_pc
    CommuteDiagVScat_pc := ARule(Compose,
       [ @(1, [Diag, VDiag]), @(2, VScat_pc) ], # <-1 <-2 o
       e -> let(s := @(2).val, f := @(1).val.element,
	   [ s, VDiag(fStretch(fCompose(f, fAdd(f.domain(), s.n, s.ofs)),
			       _roundup(s.n, s.v), s.n), s.v) ])),

    # VGath_pc * Diag
    CommuteVGath_pcDiag := ARule( Compose, # o 1-> 2->
       [ @(1, VGath_pc), @(2, [Diag, VDiag]) ],
       e -> let(g := @(1).val, f := @(2).val.element,
	   [ VDiag(fStretch(fCompose(f, fAdd(f.domain(), g.n, g.ofs)),
		            _roundup(g.n, g.v), g.n), g.v), g ])),

    #-----------------------------------------------------------
    # YSV: Below 2 rules have limitations which can be removed,
    #      But it seems that fStretch has built-in assumptions which makes it difficult
    #      to remove these limitations..

    # NOTE: Remove limitations, the rule as is, is invalid for g.ofs <> 0, or if g.N > g.v
    #        Hence the precondition
    ## Diag * IxRCVScat_pc
    CommuteDiagIxRCVScat_pc := ARule(Compose,
       [ @(1, [Diag, VDiag]), @(2, IxRCVScat_pc, e->e.ofs = 0 and not IsSymbolic(e.n) and e.n <= e.v) ], # <-1 <-2 o
       e -> let(s := @(2).val,
	   [ s, VDiag(fStretch(@(1).val.element,
			       _roundup(s.n, 2*s.v), 2*s.n), s.v) ])),

    # NOTE: Remove limitations, the rule as is, is invalid for g.ofs <> 0, or if g.N > g.v
    #        Hence the precondition
    # IxRCVGath_pc * Diag
    CommuteIxRCVGath_pcDiag := ARule( Compose, # o 1-> 2->
       [ @(1, IxRCVGath_pc, e->e.ofs = 0 and not IsSymbolic(e.n) and e.n <= e.v), @(2, [Diag, VDiag]) ],
       e -> let(g := @(1).val,
	   [ VDiag(fStretch(@(2).val.element,
		            _roundup(g.n, 2*g.v), 2*g.n), g.v), g ])),
    #-----------------------------------------------------------

    #-----------------------------------------------------------
# NOTE: unclear how to do. Rewriting rule missing: f^{n->N} o \pi^{n} -> \rho^{N} o h^{n->N}
#    # BlockVPerm * VDiag
#    CommuteBlockVPermVDiag := ARule( Compose,
#       [ @(1, BlockVPerm), @(2, VDiag) ], # o 1-> 2->
#         e -> let(vdiag := Copy(@2.val), bperm := Copy(@1.val),
#         [ VDiag(fCompose(vdiag.element, fTensor(fId(bperm.n), bperm.perm)), vdiag.v).attrs(vdiag), BlockVPerm2(bperm.n, bperm._children[1], bperm.perm) ])),
#
#    # VDiag * BlockVPerm
#    CommuteVDiagBlockVPerm := ARule( Compose,
#       [ @(1, VDiag), @(2, BlockVPerm) ], # <-1 <-2 o
#         e -> let(vdiag := Copy(@1.val), bperm := Copy(@2.val),
#        [ BlockVPerm2(bperm.n, bperm._children[1], bperm.perm), VDiag(fCompose(vdiag.element, fTensor(fId(bperm.n), bperm.perm)), vdiag.v).attrs(vdiag) ])),

# this is a hack that *DOES NOT* produce correct code
#    # BlockVPerm * VDiag: HACK that *DOES NOT* produce correct code
#    CommuteBlockVPermVDiag := ARule( Compose,
#       [ @(1, BlockVPerm), @(2, VDiag) ], # o 1-> 2->
#         e -> let(vdiag := Copy(@2.val), bperm := Copy(@1.val),
#         [ vdiag, BlockVPerm2(bperm.n, bperm._children[1], bperm.perm) ])),
#
#    # VDiag * BlockVPerm: HACK that *DOES NOT* produce correct code
#    CommuteVDiagBlockVPerm := ARule( Compose,
#       [ @(1, VDiag), @(2, BlockVPerm) ], # <-1 <-2 o
#         e -> let(vdiag := Copy(@1.val), bperm := Copy(@2.val),
#        [ BlockVPerm2(bperm.n, bperm._children[1], bperm.perm), vdiag ])),

    #--------------------------------------------------
    #   vectorize constructs
    ParDiagVPerm := ARule(Compose, [@(1, Diag), @(2, VPerm)],
        e -> [VDiag(@1.val.element, @(2).val.vlen), @(2).val]
    ),
    ParVPermDiag := ARule(Compose, [@(1, VPerm), @(2, Diag)],
        e -> [@(1).val, VDiag(@2.val.element, @(1).val.vlen)]
    ),
    ParScatVDiag := ARule(Compose,
    [[Scat, [@(0,fTensor), ..., @(2,fId)]], @(1,VDiag,e->@(2).val.size=Rows(e))],
    e -> let(ch := @(0).val.children(), X := @(1).val, v := X.v,
               [VTensor(Scat(fTensor(DropLast(ch,1), fId(Rows(X)/v))), v), X])
    ),
    ParVTensorDiag := ARule(Compose, [@(1, VTensor), @(2, Diag)],
        e -> [ @(1).val, VDiag(@2.val.element, @(1).val.vlen)]
    ),
    ParDiagVTensor := ARule(Compose, [@(1, Diag), @(2, VTensor)],
        e -> [ VDiag(@1.val.element, @(2).val.vlen), @(2).val ]
    ),

    ParVSDiag := ARule(Compose, [@(1, VS), @(2, Diag)],
        e -> [ @(1).val, VDiag(@2.val.element, @(1).val.v)]
    ),

    ParDiagVS := ARule(Compose, [@(1, Diag), @(2, VS)],
        e -> [ VDiag(@1.val.element, @(2).val.v), @(2).val ]
    )
));

RewriteRules(RulesFuncSimp, rec(

    RCData_VDupOnline := Rule( [RCData, VDupOnline],
        e -> VDupOnline(RCData(e.func.func), e.func.v)),

    RCVData_VDupOnline := Rule( [RCVData, VDupOnline],
        e -> VDupOnline(RCVData(VDup(e.func.func, 1)), e.func.v)),

    RCVDataSplitVec_VDupOnline := Rule( [RCVDataSplitVec, VDupOnline],
        e -> VDupOnline(RCVDataSplitVec(VDup(e.func.func, 1)), e.func.v)),

    fCompose_VDupOnline_FList := ARule(fCompose, [@(1, VDupOnline), @(2, FList)],
        e -> [VDupOnline(fCompose(@(1).val.func, @(2).val), @(1).val.v)]),

    fPrecompute_VDupOnline := Rule( [fPrecompute, VDupOnline],
        e -> VDupOnline(fPrecompute(e.child(1).func), e.child(1).v)),

    RCVData_VDup_CRData := Rule([RCVData, [VDup, [CRData, @(1)], @(2)]],
        e -> VDup( VData(@(1).val, 2), @(2).val) ),
));


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

# Class(DistMixin1, rec(
# #   includes := ["<include/threads.h>", "<include/smp2.h>"],
#     dist_loop := meth(self, o, i, is)
#         local v, lo, hi;
#         v := o.var;
#         lo := o.range[1];
#         hi := Last(o.range);
#         Print(Blanks(i), "for(int ", v, " = tid + ", lo, "; ", v, " <= ", hi, "; ", v, "+= ", o.p, ") {\n");
#         Print(Blanks(i+is), "int ", o.tidvar, " = ", self(o.tidexp, i+is+is, is), ";\n");
#         self(o.cmd,i+is,is);
#         Print(Blanks(i), "}\n");
#         Print(Blanks(i), When(IsBound(self.opts.smp),
#                               self.opts.smp.barrier,
#                               "//barrier(num_threads, tid, &GLOBAL_BARRIER);\n"));
#     end,
#
#     threadId := (self, o, i, is) >> Print("tid")
#
# ));

Class(DistMixin, rec(
#  includes := ["<spumacros.h>"],

   dist_loop := meth(self, o, i, is)
       Print(Blanks(i),    "{\n");
       #Print(Blanks(i+is), "unsigned int ", o.var, " = spe_info.spuid;\n");
       self(o.cmd,i+is,is);
       Print(Blanks(i),    "}\n");
       #Print(Blanks(i),    "//ALL_TO_ALL_BARRIER;\n");
       #Print(Blanks(i),    "BLOCK_ON_READ();\n");
    end,

    dist_barrier := (self,o,i,is) >> Print(Blanks(i), "BLOCK_ON_CPUDMA; ALL_TO_ALL_BARRIER;\n"),

    call := (self, o, i, is) >> Print(Blanks(i), o.args[1].id, self.pinfix(Drop(o.args, 1), ", "), ";\n"),

    fcall := (self, o, i, is) >> Print(self(o.args[1],0,0), "(", self.infix(Drop(o.args, 1), ", "), ")"),

));

Class(DistUnparser, DistMixin, CUnparserProg);
#Class(DistCellUnparser, DistMixin, CellUnparser);


text:
ImportAll(simt);

Class(mpi_rcperm, call);

Class(MPICUDACodeGenMixin, rec(
    MPIRCPrm := (self,o,y,x,opts) >> simt_block(ApplyFunc(mpi_rcperm, [y,x]::o.func.params{[1..2]}::[Length(o.func.params[3])]::o.func.params[3]))
));

Class(MPICUDAUnparserMixin, rec(
    mpi_rcperm := (self, o, i, is) >> Print(Blanks(i), "fftx_mpi_rcperm", self.pinfix(o.args, ", "), ";\n")
));

FixUpMPICUDAPerm := function(c)
    local mpicalls, mpinames, cucalls, mc, cc, ic;

    mpicalls := Collect(c, @(1, specifiers_func, e->IsBound(e.cmd.cmds) and ObjId(e.cmd.cmds[1]) in [call, mpi_rcperm]));
    mpinames := List(mpicalls, i->i.id);
    cucalls := Collect(c, cu_call);

    for mc in mpicalls do
        SubstVars(mc, FoldR(Zip2(mc.params, Filtered(cucalls, cc->cc.func = mc.id)[1].args), (a,b) -> CopyFields(a, rec((b[1].id) := V(b[2]))), rec()));
    od;

    c:= SubstTopDown(c, @(1, cu_call, e->e.func in mpinames),
        e->Filtered(mpicalls, i->i.id = @(1).val.func)[1].cmd.cmds[1]
    );

    c := SubstTopDown(c, @(1, specifiers_func, e->IsBound(e.cmd.cmds) and ObjId(e.cmd.cmds[1]) in [call, mpi_rcperm]),
        e->skip());

    c := SubstBottomUp(c, @(1, chain), e-> chain(Filtered(@(1).val.cmds, c->ObjId(c) <> skip)));

    cc := Collect(c, @(1, call, e->IsBound(e.args[1].codegen)));
    ic := List(cc, _c -> _c.args[1].codegen.init());
    c := SubstBottomUp(c, @@(1, chain, (e,cx)->IsBound(cx.func) and Length(cx.func) = 1 and cx.func[1].id = "init"),
           e->chain(@@(1).val.cmds::ic));

    c := FoldL(List(cc, d->d.args[1].codegen), (a,b)->b.data(a), chain(c.cmds));
    c := program(decl(c.free(), c));
    c := SubstBottomUp(c, @(1, decl),
        e -> decl(Filtered(@(1).val.vars, v -> v in @(1).val.cmd.free()), @(1).val.cmd));

    return c;
end;

max_stars_repo_name: spiral-software/spiral-package-mpi

max_issues_repo_name: spiral-software/spiral-package-mpi

max_forks_repo_name: spiral-software/spiral-package-mpi


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

# -------------------------------------------------------------------------------
# NOTE: Make at()/lambda() methods consistent
#
#  Either implement at via lambda consistently or the other way around
#  Currently some classes implement at() via lambda(), and some lambda() via at()
#  This is error prone!
# -------------------------------------------------------------------------------

#################################################################################
#F VData( <diagfunc>, <vlen> ) - new diagonal gen. function returning vectors
#F
#F Given a diagonal generating function F with domain N, and range (=element
#F type) T, VData(F, v) is a new "vector" diagonal function with domain N/v,
#F and range (=element type) T.vectorType(v).
#F
#F It uses vpack(_, _, ...) pseudo-instruction to pack scalars into vectors.
#F This instrucion has to be converted into architecture specific sequence
#F by the backend.
#F
Class(VData, Function, rec(
   __call__ := (self, func, v) >> WithBases(self, rec(
       func := Checked(IsFunction(func) or IsFuncExp(func), func),
       v := v,
       operations := PrintOps)),

   print := self >> Print(self.name, "(", self.func, ", ", self.v, ")"),
   rChildren := self >> [ self.func, self.v ],
   rSetChild := rSetChildFields("func", "v"),

   at := (self, i) >> let(
       n := self.func.domain(), v := self.v, f := self.func.lambda(),
       ApplyFunc(vpack, List([0..v-1], j->f.at(v*i + j)))),

   lambda := self >> let(
       n := self.func.domain(), v := self.v, f := self.func.lambda(),
       i := Ind(n/v),
       Lambda(i, ApplyFunc(vpack, List([0..v-1], j->f.at(v*i + j))))),

   domain := self >> div(self.func.domain(), self.v),
   range := self >> TVect(self.func.range(), self.v),
   free := self >> self.func.free()
));

#################################################################################
#F VDup( <diagfunc>, <vlen> ) - new diagonal gen. function returning vectors
#F
#F
#F Given a diagonal generating function F with domain N, and range (=element
#F type) T, VData(F, v) is a new "vector" diagonal function with domain N,
#F and range (=element type) T.vectorType(v). Every vector element is the
#F replicated original scalar.
#F
#F It uses vdup(<elt>, <vlen>) pseudo-instruction to duplicate scalars
#F to form vectors.  This instrucion has to be converted into
#F architecture specific sequence by the backend.
#F
Class(VDup, Function, rec(
   __call__ := (self, func, v) >> WithBases(self, rec(
       func := Checked(IsFunction(func) or IsFuncExp(func), func),
       v := v,
       operations := PrintOps)),

   print := self >> Print(self.name, "(", self.func, ", ", self.v, ")"),
   rChildren := self >> [ self.func, self.v ],
   rSetChild := rSetChildFields("func", "v"),

   at := (self, i) >> vdup(self.func.at(i), self.v),

   lambda := self >> let(
       i := Ind(self.func.domain()),
       Lambda(i, vdup(self.func.at(i), self.v))),

   domain := self >> self.func.domain(),
   range := self >> TVect(self.func.range(), self.v),
   free := self >> self.func.free()
));

#################################################################################
#F VDupOnline( <diagfunc>, <vlen> ) - is an online version of VDup.
#F
#F VDupOnline semantically the same as VDup(), but unlike VDup() it is pulled out
#F of fPrecompute and executed at compute time.
#F

Class(VDupOnline, VDup, rec(
    range := (self) >> let( t := self.func.range(),
        Cond( IsVecT(t), TVect(t.t, t.size*self.v), TVect(t, self.v)))
));

#################################################################################
#F RCVData( <cplx_vdiagfunc> ) - func returning real vecs given complex vec func
#F
#F Given a diagonal generating function F with domain N, and range of
#F of complex vectors TVect(CT, v), RCVData(F) is a new real vector diagonal
#F function with domain N and range TVect(CT.realType(), 2*v)
#F
Class(RCVData, Function, rec(
   __call__ := (self, func) >> WithBases(self, rec(
       func := Checked((IsFunction(func) or IsFuncExp(func)) and IsVecT(func.range()), func),
       operations := PrintOps)),

   print := self >> Print(self.name, "(", self.func, ")"),
   rChildren := self >> [ self.func ],
   rSetChild := rSetChildFields("func"),

   at := (self, i) >> let(
       t := self.func.range(),
       f := self.func.lambda(),
       Cond(IsVecT(t),
	    ApplyFunc(vpack, ConcatList([0..t.size-1], j -> [ re(velem(f.at(i), j)), im(velem(f.at(i), j)) ])),
	    ApplyFunc(vpack, [ re(f.at(i)), im(f.at(i)) ]))
   ),

   lambda := self >> let(j := Ind(self.domain()), Lambda(j, self.at(j))),
   domain := self >> self.func.domain(),
   range := self >> let(t := self.func.range(),
       Cond(IsVecT(t), TVect(t.t.realType(), 2*t.size), TVect(t.realType(), 2))),
   free := self >> self.func.free()
));

#################################################################################
#F RCVDataSplit(<cplx_func>, <v>)
#F
#F Given a complex function f RCVDataSplit(f, v)
#F is a new function that returns v-way vectors of the form
#F [ re0, re0, re1, re1, ... ], [-im0, im0, -im1, im1, ...]
#F
Class(RCVDataSplit, Function, rec(
   __call__ := (self, cplx_func, v) >> WithBases(self, rec(
       func := Checked(IsFunction(cplx_func) or IsFuncExp(cplx_func), cplx_func),
       v := v,
       operations := PrintOps)),

   print := self >> Print(self.name, "(", self.func, ", ", self.v, ")"),
   rChildren := self >> [ self.func, self.v ],
   rSetChild := rSetChildFields("func", "v"),

   lambda := self >> let(
       i := Ind(self.domain()),
       n := self.func.domain(),
       v := self.v,
       f := self.func.lambda(),
       vars := DropLast(f.mkVars(), 1),
       ip := idiv(i, 2),
       Lambda(vars :: [i],
       cond(imod(i, 2),
            ApplyFunc(vpack, ConcatList([0..v/2-1], j ->
	        [ -im(f.at(vars :: [(v/2)*ip+j])), im(f.at(vars :: [(v/2)*ip+j])) ])),
            ApplyFunc(vpack, ConcatList([0..v/2-1], j ->
		[  re(f.at(vars :: [(v/2)*ip+j])), re(f.at(vars :: [(v/2)*ip+j])) ]))))
   ),

   at := (self, i) >> self.lambda().at(i),

   # each cpx element -> 2 duplicated reals (r r -i i) = 4 reals -> packed into vectors of length <v>
   domain := self >> let(d := self.func.domain(),
       Cond(IsType(d), d, 4 * self.func.domain() / self.v)),

   range := self >> let(t := self.func.range(), TVect(t.realType(), self.v)),
   free := self >> self.func.free()
));

#################################################################################
#F RCVDataSplitVec( <cplx_vdiagfunc> )
#F
#F Given a v/2-way complex vector function f RCVDataSplitVec(f, v)
#F is a new function that returns v-way real vectors of the form
#F [ re0, re0, re1, re1, ... ], [-im0, im0, -im1, im1, ...]
#F
Class(RCVDataSplitVec, RCVData, rec(
   # lambda := <inherited from RCVData, implemented via at>
   #
   at := (self, i) >> let(
       v := self.func.range().size,
       f := self.func.lambda(),
       ip := idiv(i, 2),
       cond(imod(i, 2),
            ApplyFunc(vpack, ConcatList([0..v-1], j ->
		        [ -im(velem(f.at(ip), j)), im(velem(f.at(ip), j)) ])),
            ApplyFunc(vpack, ConcatList([0..v-1], j ->
		        [  re(velem(f.at(ip), j)), re(velem(f.at(ip), j)) ])))
   ),

   # each cpx element -> 2 duplicated reals (r r -i i) = 4 reals
   #  -> packed into vectors of length <v>
   domain := self >> let(d:=self.func.domain(),
       Cond(IsType(d), d, 2 * self.func.domain())),
));

#################################################################################
#F VDataRDup(<real-vector-func>)
#F
#F Returns new function of twice the vector length, with every element
#F duplicated.
#F
#F
Class(VDataRDup, RCVData, rec(
   # lambda := <inherited from RCVData, implemented via at>
   #
   at := (self, i) >> let(
       v := self.func.range().size,
       f := self.func.lambda(),
       ApplyFunc(vpack, ConcatList([0..v-1],
	           j -> [ velem(f.at(i), j), velem(f.at(i), j) ]))
   ),

   # domain doesn't change, but output vectors get longer (ie range is different)
   domain := self >> self.func.domain(),
));

#################################################################################
#F fStretch( <diagfunc>, num, den ) - new diagonal gen. function stretching the old one
#F
Class(fStretch, Function, rec(
   __call__ := (self, func, num, den) >> WithBases(self, rec(
       func := Checked(IsFunction(func) or IsFuncExp(func), func),
       den := den,
       num := num,
       operations := PrintOps)),

   print := self >> Print(self.__name__, "(", self.func, ", ", self.num, ", ", self.den, ")"),
   rChildren := self >> [ self.func, self.num, self.den ],
   rSetChild := rSetChildFields("func", "num", "den"),

   at := (self, i) >> self.lambda().at(i).eval(),

   lambda := self >> let(
       n := self.func.domain(),
       f := self.func.lambda(),
       i := Ind(n * self.num / self.den),
       i_old := self.den * idiv(i, self.num) + imod(i, self.num),
       Lambda(i, cond(leq(imod(i, self.num), self.den - 1), f.at(i_old), 0))
   ),

   domain := self >> self.func.domain() * self.num/self.den,
   range := self >> self.func.range(),
   free := self >> self.func.free()
));

#F ==========================================================================
#F VDiag(<func>, <v>)
#F
#F Vectorized diagonal. <func> must be a plain scalar function.
#F
#F Diag(f) == VDiag(f, v)
#F
Class(VDiag, Diag, rec(
    abbrevs := [
	(f, v) -> Checked(IsInt(v),
	    Cond(IsList(f), let(ff := List(f, toExpArg),
		    [ FList(UnifyTypes(List(f, x -> x.t)), f), v]),
		 IsFunction(f),
		    [ f, v ]))
    ],

    new := (self, f, v) >> SPL(WithBases(self, rec(element := f, v := v))).setDims(),

    print := (self, indent, indentStep) >> Print(self.name, "(", self.element, ", ", self.v, ")"),
    needInterleavedLeft := False,
    needInterleavedRight := False,
    transpose := self >> self,
    conjTranspose := self >> ObjId(self)(FConj(self.element), self.v),
    # self.v is not exposed in rChildren
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.v),
));

#################################################################################
#F VDiag_x_I(<func>, <v>)
#F
#F VTensor(Diag(f), v) == VDiag_x_I(f, v)
#F
Class(VDiag_x_I, VDiag, rec(
    toAMat := self >> Tensor(Diag(self.element), I(self.v)).toAMat(),
    dims := self >> self.v * Diag(self.element).dims(),
    needInterleavedLeft := False,
    needInterleavedRight := False,
    transpose := self >> self
));

#################################################################################
#F VRCDiag(<func>, <v>)
#F
#F Semantically equivalent to (I tensor L) * RCDiag * (I tensor L), and vectorized.
#F
#F Expects <func> and Input/output data to be in vector-interleaved complex format.
#F I.e. that packing is r0,r1,r2,...,r_(v-1); i0,i1,..,i_(v-1); ... into vectors.
#F
#F Example, given a complex function <f> (range = TComplex), and vector length <v>:
#F (different SIMD slots perform subsequent complex mults)
#F
#F    VRCDiag(VData(fCompose(RCData(f.domain()), fTensor(fId(d.domain()/v), L(2*v, 2))), v), v)
#F
#F Example, given a complex function <f> (range = TComplex), and vector length <v>:
#F (different SIMD slots perform SAME complex mults)
#F
#F    VRCDiag(VDup(RCData(f), v), v)
#F
Declare(VRCDiag);
Class(VRCDiag, VDiag, rec(
    print := (self, indent, indentStep) >> Print(self.name, "(", self.element, ", ", self.v, ")"),

    toAMat := meth(self)
        local elts, k, res, r, i;
        elts := self.element.tolist();
        Constraint(Length(elts) mod 2 = 0);
        res := [];
        for k in [ 1 .. Length(elts)/2 ] do
            r := EvalScalar(elts[2 * k - 1]);
            i := EvalScalar(elts[2 * k]);
            Add(res, BlockMat([[Diag(r), Diag(-i)], [Diag(i), Diag(r)]]).toAMat());
        od;
        return DirectSumAMat(res);
    end,

    toloop := self >> let(
        func := self.element, # r,i,r,i,...
        j    := Ind(idiv(func.domain(),2)),
        re   := func.at(2*j),
        im   := func.at(2*j+1),
        ISum(j, j.range,
                VScat(H(2*j.range, 2, 2*j, 1), self.v) *
                VBlk([[re, -im],
                      [im,  re]], self.v) *
                VGath(H(2*j.range, 2, 2*j, 1), self.v))),

    dims := self >> let(n:=self.element.domain() * self.v, [n,n]),  # data stores r,i pairs
    needInterleavedLeft := False,
    needInterleavedRight := False,
    transpose := self >> VRCDiag(FRConj(self.element), self.v)
));

#################################################################################
#F RCVDiag(<func>, <v>)
#F
#F Semantically equivalent to RCDiag, but vectorized.
#F Expects <func> to pack r0,i0,r1,i1,.. into vectors. Input/Output data format is same.
#F
#F Example, given a complex function <f> (range = TComplex), and vector length <v>:
#F (different SIMD slots perform subsequent complex mults)
#F
#F    RCDiag(RCData(f)) == RCVDiag(VData(RCData(f), v), v)
#F
#F Example, given a complex function <f> (range = TComplex), and vector length <v>:
#F (different SIMD slots perform SAME complex mults)
#F
#F    RCDiag(RCData(f)) == RCVDiag(RCVData(VDup(f, v/2)), v)
#F
#F Example, given a complex vector function <f> (range = TVect(TComplex, v/2))
#F
#F    RCVDiag(RCVData(f), v),
#F
Class(RCVDiag, VDiag, rec(
    print := (self, indent, indentStep) >> Print(self.name, "(", self.element, ", ", self.v, ")"),

    toAMat := self >> RCDiagonalAMat(
	ConcatList(self.element.tolist(), x -> Cond(IsValue(x), x.v, List([1..self.v], i->velem(x,i)))),
	IdentityPermAMat(2)),

    dims := self >> let(n:=self.element.domain() * self.v, [n,n]),  # data stores r,i pairs
    needInterleavedLeft := True,
    needInterleavedRight := True,

    transpose := self >> self  # this is invalid! need to conjugate complex function
));

#################################################################################
#F RCVDiagSplit(<func>, <v>)
#F
#F Semantically equivalent to RCDiag, but vectorized, and expects different <func>
#F data format than RCVDiag.
#F
#F Expects <func> to pack r0,r0,r1,r1,...; -i0,i0,-i1,i1,...; ... into vectors.
#F Expects Input/Output data to pack r0,i0,r1,i1,.. into vectors.
#F
#F Example, given a complex function <f> (range = TComplex), and vector length <v/2>:
#F (different SIMD slots perform SAME complex mults)
#F
#F  RCVDiagSplit(RCVDataSplitVec(VDup(f, v/2)), v);
#F
#F Example, given a complex function <f> (range = TComplex), and vector length <v/2>:
#F (different SIMD slots perform different complex mults)
#F
#F  RCVDiagSplit(RCVDataSplit(f, v), v);
#F
Class(RCVDiagSplit, VDiag, rec(
    print := (self, indent, indentStep) >> Print(self.name, "(", self.element, ", ", self.v, ")"),

    toAMat := meth(self)
        local elts, k, res, r, i, p;
        elts := self.element.tolist();
        Constraint(Length(elts) mod 2 = 0);
        res := [];
        p := Tensor(I(self.v/2), J(2));
        for k in [ 1 .. Length(elts)/2 ] do
            r := EvalScalar(elts[2 * k - 1]);
            i := EvalScalar(elts[2 * k]);
            Add(res, SUM(Diag(r), Diag(i)*p).toAMat());
        od;
        return DirectSumAMat(res);
    end,

    dims := self >> let(n:=self.element.domain() * (self.v/2), [n,n]),  # data stores r,i pairs
    needInterleavedLeft := True,
    needInterleavedRight := True,
    transpose := self >> self   # this is invalid! need to conjugate complex function
));

#################################################################################
#F VScale(<spl>, <constant>, <v>)
#F
#F Vecorized version of Scale
#F
Declare(VScale);
Class(VScale, BaseMat, SumsBase, rec(
    new := (self, spl, s, v) >> SPL(WithBases(self,
        rec(dimensions := spl.dims(), _children:=[spl], s:=s, v := v))),
    #-----------------------------------------------------------------------
    dims := self >> self.dimensions,
    child := (self, i) >> self._children[i],
    rChildren := self >> [self._children[1], self.s, self.v],
    rSetChild := meth(self, n, what)
        if n=1 then self._children[1] := what;
        elif n=2 then self.s := what;
        elif n=3 then self.v := what;
        else Error("<n> must be in [1..3]");
        fi;
    end,
    #-----------------------------------------------------------------------
    print := (self, i, is) >> Print(self.name, "(", self.child(1).print(i+is,is), ", ", self.s, ", ", self.v, ")"),
    unroll := self >> self,
    transpose := self >> VScale(self.child(1).transpose(), self.s, self.v),
    isReal := self >> self.child(1).isReal(),
    #-----------------------------------------------------------------------
    toAMat := self >> AMatMat(MatSPL(Scale(self.s, self.child(1))))
));

#################################################################################
# Below are the version of above objects hacked up for online twiddle computation
# The difference is in the codegen, as it uses a different style of underlying
# function, that is recursive.
#
# NOTE: Find a cleaner and more general solution
# NOTE: document this
# NOTE: .toAMat() does not work for the below object, because self.func is different
#
Class(RCVOnlineDiag, RCVDiag);
Class(VOnlineDiag, VDiag);
Class(VOnlineDiag_x_I, VDiag_x_I);


text:
DotProduct := function(u, v)
	return u*v;
end;

CrossProduct := function(u, v)
	return [
		u[2]*v[3] - u[3]*v[2],
		u[3]*v[1] - u[1]*v[3],
		u[1]*v[2] - u[2]*v[1] ];
end;

ScalarTripleProduct := function(u, v, w)
	return DotProduct(u, CrossProduct(v, w));
end;

VectorTripleProduct := function(u, v, w)
	return CrossProduct(u, CrossProduct(v, w));
end;

a := [3, 4, 5];
b := [4, 3, 5];
c := [-5, -12, -13];

DotProduct(a, b);
# 49

CrossProduct(a, b);
# [ 5, 5, -7 ]

ScalarTripleProduct(a, b, c);
# 6

# Another way to get it
Determinant([a, b, c]);
# 6

VectorTripleProduct(a, b, c);
# [ -267, 204, -3 ]


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

#F AllocBuffersSMP(pfx, code, map, opts)
#F
#F Returns [newbuffers, allocation_code, code]
#F and fills the <map> with the remapping from original array-typed variables
#F to the new dynamically allocated pointer-typed variables.
#F
#F if opts.zeroallocate is true, the buffers are zeroed out
#F if opts.buffPrefixExtra is set, it's appended to pfx
AllocBuffersSMP := function(pfx, code, map, opts)
    local nameprefix, newbuffers, buffers, b, bufmap, v, nthreads, tid, newv, alloc_code;
    ## Pull out the buffer declarations
    [buffers, code] := PullBuffersSMP(code, IsArrayT);
    alloc_code := [];
    newbuffers := [];
	if IsBound(opts.buffPrefixExtra) then
		nameprefix := Concat(pfx, String(opts.buffPrefixExtra));
	else
		nameprefix := pfx;
	fi;
    for b in buffers do
        [v, nthreads, tid] := b;
        if nthreads = 1 then
            newv := var.fresh_t(nameprefix, TPtr(v.t.t));
        else
            newv := var.fresh_t(Concat("thr", nameprefix), TPtr(v.t.t));
        fi;
        if IsBound(v.value) then
            newv.value := v.value;
        fi;
        Add(newbuffers, newv);
        if IsBound(opts.zeroallocate) and opts.zeroallocate then
            Add(alloc_code, zallocate(newv, TArray(v.t.t, v.t.size*nthreads)));
        else
            Add(alloc_code, allocate(newv, TArray(v.t.t, v.t.size*nthreads)));
        fi;
        map.(v.id) := When(nthreads=1, newv, newv + tid*v.t.size);
    od;
    return [newbuffers, alloc_code, code];
end;

#F UnifyBuffersSMP(pfx, code, map, ubuf)
#F
#F Returns [ubuf_size, ptr_init_code, code]
#F and fills the <map> with the remapping from original array-typed variables
#F to the new dynamically allocated pointer-typed variables.
#F
#F ptr_init_code initializes the pointer-typed variables as pointers to
#F sections of <ubuf>
#F
UnifyBuffersSMP := function(pfx, code, map, ubuf, opts)
    local map, vars, alloc, ofs, usize, ptr_init_code;
    [vars, alloc, code] := AllocBuffersSMP(pfx, code, map,opts);
    ofs   := ScanL(alloc, (prev, c) -> prev + c.exp.size * sizeof(c.exp.t), 0);
    usize := Sum(alloc, a -> a.exp.size * sizeof(a.exp.t));

    ptr_init_code := List([1..Length(vars)], i -> assign(vars[i], tcast(vars[i].t, (ubuf + ofs[i]))));
    return [ usize, ptr_init_code, code ];
end;

Class(RecCodegen, RecCodegenMixin, SMPCodegenMixin, DefaultCodegen, rec(
    Formula := meth(self, o, y, x, opts)
        local code, datas, prog, params, init_code, destroy_code, codelet_codes, codelet_recs,
          datvars, dalloc, bufvars, bufalloc, map, ignore, num_threads, smp, io, _data, v;

        [x, y] := self.initXY(x, y, opts);
    o := o.child(1);
    params := Set(Collect(o, param));
    datas := Collect(o, @(1, var, e->IsBound(e.init)));
    smp := Collect(o, SMPSum);
    num_threads := When(smp=[], 1, smp[1].nthreads);
#    if not ForAll(smp, x->x.nthreads=num_threads) then Error("Non-uniform num_threads in SMPSum's"); fi;
    io := When(x=y, [x], [y, x]);

    ## Generating code : codelets
    codelet_recs := CompileCodelets(o, opts);
    codelet_codes := List(codelet_recs, clrec ->
        func(TVoid, clrec.name, Concatenation([Y, X], clrec.params), clrec.code));

    map := tab();
    ## Generating code : main body
    code := SReduce(self(o, y, x, opts), opts);
    code := BlockUnroll(code, opts);
    code := DeclareHidden(code); # NOTE: do I need this?
    code := func(TVoid, "transform", Concatenation(
        When(IsBound(opts.subParams), opts.subParams, []), params, io),
                When(num_threads=1, code,
                                    smp_fork(num_threads, code)));

    [bufvars, bufalloc, code] := AllocBuffersSMP("buf", code, map, opts);

    ## Generating code : initialization
    [datvars, dalloc, ignore] := AllocBuffersSMP("dat", decl(datas, skip()), map, opts);
    init_code := func(TVoid, "init", [],
        chain(bufalloc, dalloc, List(datas, x -> SReduce(x.init, opts))));
	destroy_code := func(TVoid, "destroy", [], skip());

    for v in bufvars do Add(v.t.qualifiers, "static"); od;
    for v in datvars do Add(v.t.qualifiers, "static"); od;
    _data := When(not IsBound(opts.smp) or (IsBound(opts.smp) and IsBound(opts.smp.OmpMode) and opts.smp.OmpMode = "for"), (i,j,k)->k, (i,j,k)->data(i, j, k));

    prog := program(
        codelet_codes,
        decl(Concatenation(datvars, bufvars),
        _data(var("NUM_THREADS", TInt), V(num_threads),
            SubstVars(chain(init_code, code, destroy_code), map))));
    prog.dimensions := o.dimensions;
    return prog;
    end,
));


text:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

Class(JoinDirectSums, RuleSet);
RewriteRules(JoinDirectSums, rec(
    combine := ARule(Compose, [@(1, DelayedDirectSum), @(2, DelayedDirectSum,
        e -> let(c1:=@(1).val.children(), c2:=e.children(), Length(c1) = Length(c2) and ForAll([1..Length(c1)], i->c1[i].dims()[2] = c2[i].dims()[1]))) ],
        e -> let(c1:=@(1).val.children(), c2:=@(2).val.children(), [ DelayedDirectSum(List([1..Length(c1)], i -> c1[i] * c2[i])) ])
    )
));

Class(TerminateDirectSums, RuleSet);
RewriteRules(TerminateDirectSums, rec(
    terminate := Rule(@(1, DelayedDirectSum), e-> DirectSum(@(1).val.children()).sums())
));