######################################################################
 ## Trigonometric Polynomials
 ######################################################################
 ## Jamie Mulholland & Michael Monagan
 ## August 17,2000
 ## trigpolylib.txt  (for Maple)
 ## ***Required library ***  setpartition.txt  (included below)
 ######################################################################
 ##  This library contains procedures for trigonometric polynomials.
 ##  There are five procedures: trig_degree, trig_divide,
 ##  trig_irreduc, trig_simplify, trig_GCD, trig_factors.
 ######################################################################

 phi := proc(a,s,c) algsubs(s^2=1-c^2,expand(a)) end:

 ##======================================================================
 ## Main procedures for trigonometric polynomials:
 ##======================================================================

printf("trig_degree(p,s,c) - compute TD(p)\n");
 trig_degree:=proc(p,s,c)
  ##  input:  A polynomial in s and c.
  ## output:  The trigonometric degree.
  RETURN(degree(phi(p,s,c),{s,c}));
 end:


printf("trig_divide(a,b,s,c) - test if b divides a\n");
printf("trig_divide(a,b,s,c,q) - and compute the quotient\n");
 trig_divide:=proc(a,b,s,c,q)
    ## input:  Two polynomials a,b in s and c.
    ##output:  true/false depending on whether b divides a or not.  An
    ##         optional 5th argument can be used to return the quotient
    ##         if b divides a.
   local a1,b1,da,db,q1,dq;
   a1:=convert_to_t(a,s,c,t,'da');
   b1:=convert_to_t(b,s,c,t,'db');
   a1:=a1*(1+t^2)^da; b1:=b1*(1+t^2)^db;
   if divide(a1,b1,'q1') and degree(a1,t)-degree(b1,t)<=2*(da-db) then
     if nargs = 5 then q:=custom_convert_to_sc(q1,da-db,t,s,c); fi;
     RETURN(true);
   else RETURN(false); fi;
 end:

printf( "trig_irreduc(p,s,c) - test if p is irreducible\n");
 trig_irreduc:=proc(p,s,c)
    ## input:  A trigonometric polynomial p is s,c.
    ##  output: 'true' if p is irreducible in Q[s,c]/<
    ##  s^2+c^2-1 > , otherwise 'false'.
   local p1,dp,t,fp1;
   # check if p is constant
   if diff(p,s)=0 and diff(q,c)= 0 then RETURN(true); fi;
   p1:=algsubs(s^2=1-c^2,p);
   dp:=degree(p1);
   p1:=normal(subs(s=2*t/(1+t^2),c=(1-t^2)/(1+t^2),p1));
   p1:=p1*(1+t^2)^dp;
   if degree(p1) < 2*dp-1 then RETURN(false) fi;
   if dp=1 or irreduc(p1) then RETURN(true) fi;

   # p1 is reducible by this point
   fp1:=factors(p1); fp1 := fp1[2];
   nops(fp1)=2 and degree(fp1[1][1],t) mod 2=1 and degree(fp1[2][1],t) mod 2=1
 end:

printf("trig_simplify(r,s,c) - simplify a ratio of trig polynomials\n");
 trig_simplify:=proc(r,s,c)
 # input: A rational trigonometric polynomial a/b.
 #output: A rational trigonometric polynomial abar/bbar whose
 #        numerator and denominator are of smallest trig degree
 #        and such that abar/bbar = a/b.
   local a,b,a1,b1,da,db,g1,ca,cb,l;
   a:=numer(r); b:=denom(r);
   a1:=convert_to_t(a,s,c,t,'da');
   a1:=a1*(1+t^2)^da;
   b1:=convert_to_t(b,s,c,t,'db');
   b1:=b1*(1+t^2)^db;
   g1:=gcd(a1,b1,'ca','cb');
   if da<=db then ca:=ca*(1+t^2)^(db-da);
   else cb:=cb*(1+t^2)^(da-db); fi;
   l:=max(ceil(degree(ca)/2),ceil(degree(cb)/2));
   RETURN(normal(custom_convert_to_sc(ca,l,t,s,c)/
                 custom_convert_to_sc(cb,l,t,s,c)));
 end:


printf("trig_gcd(p,q,s,c) - compute the trig GCD g of p and q\n");
printf("trig_gcd(p,q,s,c,pbar,qbar) - and the cofactors p/g and q/g\n");
trig_gcd:=proc(p,q,s,c,pbar,qbar)
##  input:  Two trigonometric polynomials p and q.
## output:  The trigonometric GCD of p and q.
   local p1,q1,g1,dpbar,dqbar,GCD,leftover,p1bar,q1bar,ofactors,f,ic;
   p1:=convert_to_t(p,s,c,t,'dp');
   q1:=convert_to_t(q,s,c,t,'dq');
   p1:=p1*(1+t^2)^dp; q1:=q1*(1+t^2)^dq;
   g1:=gcd(p1,q1,'p1bar','q1bar');
   dpbar:=degree(p1bar); dqbar:=degree(q1bar);
   leftover:=min(dp-ceil(dpbar/2),dq-ceil(dqbar/2));
   if degree(g1)>2*leftover then
     ofactors:=sort(eofactors(g1,t)[2],proc(x,y) evalb(degree(x[1])<=degree(y[1])) end);
     f:=ofactors[1];
     divide(g1,f[1],'g1'); p1bar:=p1bar*(f[1]); q1bar:=q1bar*(f[1]);
     dpbar:=degree(p1bar); dqbar:=degree(q1bar);
     leftover:=min(dp-ceil(dpbar/2),dq-ceil(dqbar/2));
   fi;
   GCD:=custom_convert_to_sc(g1,leftover,t,s,c);
   # Make the GCD monic over Z
   GCD := monic(GCD,[s,c],'ic');
   if nargs = 6 then
     pbar:=ic*custom_convert_to_sc(p1bar,dp-leftover,t,s,c);
     qbar:=ic*custom_convert_to_sc(q1bar,dq-leftover,t,s,c);
   fi;
   RETURN(GCD);
 end:



printf("trig_factors(p,s,c) - output all trig factorizations of p\n");
trig_factors:=proc(p,s,c)
##  input:  A trigonometric polynomial p (polynomial in s and c).
## output:  All factorizations of p.
   local t,p1,nf,f,dp,eflist,oflist,numsing,const,tempf,etrigfactors,trigfactors,
         allfactorizations,numodd,K,bijec,intlist,otrigfactorizations,i,
         leftovers,choice_singles,j,singles,trigfactorizations,nc;
   p1:=convert_to_t(p,s,c,t,'dp');
   p1:=p1*(1+t^2)^dp;
   nf:=num_freefactors(degree(p1),dp);
   eflist,oflist:=eofactors(p1,t,'const');
   numodd:=add(f[2],f=oflist);
   #######################################################################
   ### numsing is the maximum possible number of odd degree factors there
   ### can be which are left single (not paired up).  Whereas K is the
   ### exact maximum number of odd degree factors that can be left single.
   #######################################################################
   if numodd mod 2 = 0 then numsing:=2*nf; else numsing:=2*nf+1; fi;
   K:=min(numodd,numsing);

   allfactorizations:={};
   etrigfactors:= convert_list_to_sc(eflist,t,s,c);
   #etrigfactors,const:= makemonic(etrigfactors,const);
   #etrigfactors:=op(etrigfactors);

   ##########################
   ### Simple case first.
   ##########################
   if oflist=[] then
     etrigfactors,const:=makemonic(etrigfactors,const);
     RETURN([const,etrigfactors]);
   fi;

   bijec,intlist:=natmap(oflist);
   otrigfactorizations:=NULL;

   if numodd mod 2 = 0 then
       for i from 0 to K by 2 do
          leftovers:=nf-i/2;
          if i=0 then
            if leftovers<>0 then
               otrigfactorizations:=
                 op(map(proc(x) [op(x),[1/2+c/2,leftovers]]end,to_trig(intlist,[],bijec,t,s,c)));
            else
               otrigfactorizations:=op(to_trig(intlist,[],bijec,t,s,c));
            fi;
          else
            choice_singles:=combinat[choose](intlist,i);
            for j from 1 to nops(choice_singles) do
              singles:=choice_singles[j];
              if leftovers<>0 then
               otrigfactorizations:=otrigfactorizations,
                 op(map(proc(x) [op(x),[1/2+c/2,leftovers]]end,to_trig(intlist,singles,bijec,t,s,c)));
              else
               otrigfactorizations:=otrigfactorizations,
op(to_trig(intlist,singles,bijec,t,s,c));
              fi;
            od;
          fi;
       od;
       #trigfactorizations:=map(proc(x) [const,[etrigfactors,op(x)]] end,[otrigfactorizations]);
       otrigfactorizations := [otrigfactorizations];
print(numberofoddfactorizations=nops(otrigfactorizations));
       trigfactorizations:=map(proc(odd) local c,of,ef;
                of,c := makemonic(odd,const);
                ef,c := makemonic(etrigfactors,c);
                [c,[op(ef),op(of)]]
                               end,otrigfactorizations);
    else
      for i from 1 to K by 2 do
          leftovers:=nf-(i-1)/2;
          choice_singles:=combinat[choose](intlist,i);
          for j from 1 to nops(choice_singles) do
              singles:=choice_singles[j];
              if leftovers<>0 then
               otrigfactorizations:=otrigfactorizations,
                 op(map(proc(x) [op(x),[1/2+c/2,leftovers]] end,to_trig(intlist,singles,bijec,t,s,c)));
              else
               otrigfactorizations:=otrigfactorizations,
op(to_trig(intlist,singles,bijec,t,s,c));
              fi;
          od;
       od;
       #trigfactorizations:=map(proc(x) [const,[etrigfactors,op(x)]] end,otrigfactorizations);
       otrigfactorizations := [otrigfactorizations];
print(numberofoddfactorizations=nops(otrigfactorizations));
       trigfactorizations:=map(proc(odd) local c,of,ef;
                ef,c := makemonic(etrigfactors,const);
                of,c := makemonic(odd,c);
                [c,[op(ef),op(of)]]
                               end,otrigfactorizations);
   fi;

                                 #        end,trigfactorizations);
   RETURN({op(trigfactorizations)});
end:

##===============================================================
##===============================================================
##Sub procedures:
##===============================================================
 ## PROCEDURE: num_freefactors
 ##  input:   An integer dp (representing the degree of a poly p in Q[t])
 ##           and an integer n.
 ## output:  (dp-2*n)/2; if dp; even , (dp-2*n-1)/2; if dp; odd.
 ##           This number corresponds to the number of (1+t^2) factors in
 ##           the denominator of psi[t](p) which are free to be moved
 ##           around.
num_freefactors:=proc(dp,n)
   RETURN(floor(1/2*(2*n-dp)));
end:


## PROCEDURE: natmap (natural mapping)
## input:  A list of n irreducible pairwise relatively prime polynomials with
## multiplicities.
## output:  A bijection from the set or irreducible polynomials and [n], also
##         a multiset (list)of integers between 1 and n where then number or
##         times the integer i appears is the multiplicity of the polynomial
##         to which it corresponds.
natmap:=proc(L)
   RETURN([seq(i=L[i][1],i=1..nops(L))],[seq(i$L[i][2],i=1..nops(L))]);
end:

##__________________________________________________________________________
 ##Operations for converting between t and s,c:
 ##
 ## PROCEDURE: to_trig
 ## input: Two lists of integers, intlist, singles, a list of
 ##        [polynomial,integer],and an integer n.
 ## output: A list of lists of [trigonometric polynomial, multiplicity]
 ##         where a trig poly corresponds to each of the integers in the
 ##         list 'singles' or a pair of integers in 'intlist' .
to_trig:=proc(intlist::list(integer),singles::list(integer),bijec,t,s,c)
   local pairings,polypairing,mpairs,mpairstrig,singlestrig,
         alldifftrigfact,const;
   pairings:=setpartition(minuslist(intlist,singles),2);
   singlestrig:=convert_list_to_sc(findmultiplicity(subs(bijec,singles)),t,s,c);
   if pairings=[] then RETURN([singlestrig]);
   else
     polypairing:=subs(bijec,pairings);
     mpairs:=map(multpairs,polypairing);
     mpairstrig:=map(proc(x) convert_list_to_sc(x,t,s,c) end,mpairs);
     alldifftrigfact:=map(proc(x) [op(x),op(singlestrig)] end,mpairstrig);
     #alldifftrigfact:=map(proc(x) local c,l; l,c:=makemonic(x,1);
          #RETURN([c,[l]]); end,alldifftrigfact);
     RETURN(alldifftrigfact);
   fi;
end:

## PROCEDURE: convert_to_t
printf("psi[t](a,s,c) - psi_t(a(s,c))\n");
psi := proc(a,s,c) local t;
   if not type(procname,indexed) then ERROR("invalid input") fi;
   t := op(procname);
   convert_to_t(a,s,c,t)
end:

convert_to_t:=proc(a,s,c,t,d)
    local a1,da;
    a1:=algsubs(s^2=1-c^2,expand(a));  da:=degree(a1);
    if nargs = 5 then d:=da; fi;
    RETURN(normal(subs(s=2*t/(1+t^2),c=(1-t^2)/(1+t^2),a)));
end:


## PROCEDURE: convert_list_to_sc
##
## **Works only for lists of even degree polynomials.
convert_list_to_sc:=proc(plist,t,s,c)
 map(proc(x) [convert_to_sc(op(1,x),t,s,c),op(2,x)] end,plist);
end:


##PROCEDURE: convert_to_sc
## input:  A polynomial p in t.
##output:  The tightest polynomial in s,c which corresponds to p.
##         ( By tightest I mean the polynomial in s,c whose image
##         under psi[t] is  p/(1+t^2)^n where n=ceil(deg(p)/2)  ).

printf("convert_to_sc(a,t,s,c,d) - convert a(t)/(1+t^2)^d to s/c\n");
convert_to_sc:=proc(p,t,s,c,deg::nonnegint) option
remember,system;
 local d,r,A,ans;
   if nargs=4 then d:=ceil(degree(p)/2); else d := deg fi;
   r:=resultant(p-A*(1+t^2)^d,s*t-1+c,t);
   r:=algsubs(s^2=1-c^2,r);
   ans:=solve(r,A);
   RETURN(ans);
end:

## PROCEDURE: custom_convert_to_sc
## input: A polynomial p in Q[t] and a positive integer d.
##output: The polynomial in s and c corresponding to p/(1+t^2)^d.
custom_convert_to_sc:=proc(p,d::integer,t,s,c) option
remember,system;
   local n,r,A,ans;
   if 2*d<degree(p,t) then ERROR(`in converting to trig`); fi;
   r:=resultant(p-A*(1+t^2)^d,s*t-1+c,t);
   r:=algsubs(s^2=1-c^2,r);
   ans:=solve(r,A);
   RETURN(ans);
end:

##_______________________________________________________________
##List operation:
##
##PROCEDURE: minuslist
minuslist:=proc(L1::list(integer),L2::list(integer))
   local i,j,L;
   L:=L1;
   for i from 1 to nops(L2) do
       if member(L2[i],L,'j') then
          L:=subsop(j=NULL,L);
       fi;
   od;
 RETURN(L);
end:

##_______________________________________________________________
##Polynomial list operations:
##
##PROCEDURE: eofactors
##  input: A polynomial p in Q[t].
## output:  Two maple-type multi-lists: one containig all nonconstant
##          even degree factors of p and one containing all odd
##          degree factors p.
eofactors:=proc(p,t,const)
    local p1,fl,elist,olist;
    p1:=factors(p);
    fl:=op(2,p1);
    if nargs = 3 then const:=op(1,p1); fi;
    elist:=select(proc(x) degree(op(1,x),t) mod 2 =0 end,fl);
    olist:=select(proc(x) degree(op(1,x),t) mod 2 =1 end,fl);
    RETURN(elist,olist);
end:

## PROCEDURE: multpairs
## input:  A list of  pairs of polynomials.
## output: A list of pairs , [polynomial, integer], where the first entry
##        is a polynomial corresponding to the product of a pair of polynomials
##        in the input list and the second entry is the number of pairs in the
##        input list which have this product.
multpairs:=proc(L::list(list(polynom)))
   local plist,multip,templist,n,finallist,i,j;
   plist:=map(proc(x) x[1]*x[2] end,L);
   ##############################
   ## now find multiplicities
   ##############################
   finallist:=findmultiplicity(plist);
   RETURN(finallist);
end:

## PROCEDURE: findmultiplicity
##  input:  A list (multiset) of polynomials.
## output: A list of pairs [polynomial,multiplicity].
findmultiplicity:=proc(L)
   local mlist,templist,n,i,j;
   templist:=L; mlist:=NULL;
   for i from 1 to nops(L) do
    n:=0;
    for j from 1 to nops(templist) do
      if normal(L[i]-templist[j])=0 then
         n:=n+1
      fi;
    od;
    templist:=remove(proc(x) simplify(L[i]-x)=0 end,templist);
    if n<>0 then mlist:=mlist,[L[i],n];fi;
   od;
   RETURN([mlist]);
end:

##________________________________________________________________
##Operations to convert to monic polynomials:
##
## PROCEDURE: makemonic (over Z)
## input:  A list of [trigpoly,multiplicity] and a constant ,c.
##output:  A list of [monicpolynomial,multiplicity] and the product of all
##         lcoeffs and oldconst.
makemonic:=proc(plist::list,oldconst)
   local n,newconst,newplist,i;
   n:=nops(plist);
   newconst:=oldconst;
   newplist:= map(proc(x) [monic(op(1,x),[s,c],'lc'),op(2,x),lc^op(2,x)] end, plist);
   for i from 1 to n do
     newconst:=newconst*op(3,op(i,newplist));
     newplist:=subsop(i=[op(1,op(i,newplist)),op(2,op(i,newplist))],newplist);
   od;
   RETURN(newplist,newconst);
end:

## PROCEDURE: monic (over Z)
##  input:  A polynomial p with indeterminants in L.
## output:  A monic polynomial corresponding to p.  An optional 3rd argument
##          will return the lcoeff of p.
monic:=proc(p::polynom,L::list,ic)
   local r;
   r:=icontent(p)*sign(p,L);
   if nargs=3 then ic:=r; fi;
   RETURN(p/r);
end:

 #################################################
 ##  Procedure:  setpartition
 #################################################
 ##  Mike Monagan & Jamie Mulholland
 ##  August 16, 2000
 ##  setpartition.txt  (for Maple)
 #################################################
 ## Examples:
 ## (i)  setpartition([1,1,2,3],2) -> [[[1,1],[2,3]],[[1,2],[1,3]]]
 ## (ii) setpartition({1,2,3,4},2) -> {{{1,2},{3,4}},{{1,3},{2,4}},
 ##                                    {{1,4},{2.3}}}
 #################################################

setpartition := proc(S::{list,set},k::posint)
##setpartition([1,1,2,3],2) -> [ [[1,1],[2,3]], [[1,2],[1,3]] ]
   local P,N,T,FI,BI,F,G,i,s,t,u;
   #Default case
   if nops(S)=0 then RETURN(S); fi;
   if type(S,set) then P := setpartition([op(S)],k);
      {seq( {seq( {seq(u,u=t)},t=s)}, s = P )};
   else
      N := nops(S);
      if N mod k <> 0 then ERROR("wrong number of elements in list") fi;
      T := {op(S)};
      T := sort([op(T)]);
      FI := {seq(T[i]=i,i=1..nops(T))};
      BI := {seq(i=T[i],i=1..nops(T))};
      F := sort(subs(FI,S));
      P := NULL;
      while F <> FAIL do G := [seq(F[(i-1)*k+1..i*k],i=1..N/k)];
         P := P, G; F := kpart(F,k); od;
      subs(BI,[P]);
   fi
end:

kpart := proc(L,k)
   local F,R,S,i,j,m,x,y;
   if nops(L)=k then RETURN(FAIL) fi;
   F := L[1..k];
   R := L[k+1..-1];
   S := kpart(R,k);
   if S<>FAIL then RETURN( [op(F),op(S)] ) fi;
   R := sort(R);
   for i from k by -1 to 1 do
       x := F[i];
       m,y := getnext(x,R);
       if m=0 then next fi;
       F := subsop(i=y,F);
       R := subsop(m=NULL,R);
       R := insert(x,R);
       for j from i+1 to k do
           x := F[j];
           R := insert(x,R);
           m,y := getnext(y-1,R);
           if m=0 then break; fi;
           F := subsop(j=y,F);
           R := subsop(m=NULL,R);
       od;
       if m=0 then next fi;
       S := nexttuple(F,R,k);
       if S=FAIL then RETURN(FAIL); fi;
       RETURN( [op(F),op(S)] );
   od;
   FAIL;
end: insert := proc(x,R) local k;
    for k to nops(R) while x>R[k] do od;
    if k>nops(R) then [op(R),x] else [op(1..k-1,R),x,op(k..-1,R)] fi;
end: getnext := proc(x,R) local k;
    for k to nops(R) do if R[k]>x then RETURN(k,R[k]) fi; od;
    RETURN(0,x);
end: nexttuple := proc(f,r,k)
    local T,R,i,x,y,m;
    if nops(r)=0 then RETURN([]) fi;
    if r[1]<f[1] then RETURN(FAIL) fi;
    T := [];
    R := r;
    for i to k do
        x := f[i];
        m,y := getnext(x-1,R);
        if m=0 then RETURN(FAIL) fi;
        T := [op(T),y];
        R := subsop(m=NULL,R);
        if y>x then break fi;
    od;
    for i from i+1 to k do
        m,y := getnext(y-1,R);
        if m=0 then RETURN(FAIL) fi;
        T := [op(T),y];
        R := subsop(m=NULL,R);
    od;
    R := nexttuple(T,R,k);
    if R=FAIL then RETURN(FAIL) fi;
    [op(T),op(R)];
end: