//Research supported in part by OTKA PD75264 and János Bolyai Research Scholarship of the 
//Hungarian Academy of Sciences
//
//This MAGMA code implements the algorithm given by Fuchs and Pethő. It was used to compute
//numerical examples in the paper "On Composite Rational Functions" by Pethő and Tengely.
//
//The function works as follows:
//t denotes the number of roots/poles of g
//n denotes the number of roots/poles of f
//tipus is in {0,1}, it is 0 if S_infty is empty, otherwise it is 1.
//So one way to call the function is CFunc(2,3,1);
//we obtain a set of cardinality 18, an element of this set is:
//<[
//        { 1 },
//        { 2 },
//        { 3 }
//    ], <1, 2, 2>, [
//        X[1] - X[3] - 1/4*X[4] + 1/4*X[5],
//        X[2] - X[3] - 1/2*X[4] + 1/2*X[5]
//    ]>
//
//this element is a record having 3 parts
//1st: a list contaning the partitions of [1,2,..n]
//2nd: the exponent vector <l_1,l_2,..,l_n>
//3rd: the system of equations defining the variety
//
//Other way to use the function is to use the optional
//parameters PSet and exptup.
//PSet is a list contaning an appropriate partition of [1,2,..n],
//exptup is a list containing exponents l_1 up to l_n.
//In this way we get a system of equation corresponding to
//the given partition and exponent vector e.g.:
//CFunc(2,3,1:PSet:=[{1},{2},{3}]);
//the output is:
//{
//    <[
//        { 1 },
//        { 2 },
//        { 3 }
//    ], <1, 2, 2>, [
//        X[1] - X[3] - 1/4*X[4] + 1/4*X[5],
//        X[2] - X[3] - 1/2*X[4] + 1/2*X[5]
//    ]>,
//    <[
//        { 1 },
//        { 2 },
//        { 3 }
//    ], <2, 2, 1>, [
//        X[1] - X[3] - 1/4,
//        X[2] - X[3] + 1/4,
//        X[4] - X[5] + 1
//    ]>,
//    <[
//        { 1 },
//        { 2 },
//        { 3 }
//    ], <2, 1, 2>, [
//        X[1] - X[3] - 1/2,
//        X[2] - X[3] - 1/4,
//        X[4] - X[5] - 1
//    ]>
//}
//That is we obtained all systems with the given partition [{1},{2},{3}].
//In a similar way one can compute all systems with exponent vector [[1,2,2]]:
//CFunc(2,3,1:exptup:=[[1,2,2]]);
//
//If we have a system of equations, then sometimes it is possible to compute
//parameterized solutions. To do so we may use Maple or Sage, using Maple it goes
//as follows:
//solve({X[1] - X[3] - 1/2,X[2] - X[3] - 1/4,X[4] - X[5] - 1});
//the solution is {X[4] = 1 + X[5], X[3] = X[2] + 1/2, X[1] = 1/4 + X[2], X[2] = X[2], X[5] = X[5]},
//that is we have two parameters X[2] and X[5].
//

CFunc:=function(t,n,tipus:PSet:=[],exptup:=[])
/////////////////////////////////////////////////////
particio:=function(s)
N:=&+s;
sset:=SequenceToSet(s);
H:=[];

for k in sset do
 for l in Subsets({1..N},k) do
  Include(~H,[l]);
  s1:=Exclude(s,k);
  for t in [1..#s1] do Hnew:=[];
   for h in H do
    for l1 in Subsets({1..N} diff &join h,s1[t]) do
     Include(~H,Include(h,l1));
     Exclude(~H,h);
    end for;
   end for;
  end for;
 end for;
end for;

return H;
end function;
/////////////////////////////////////////////////////

/////////////////////////////////////////////////////
lf:=function(nn)
nset:=[1..nn-1];
cA:=car<nset,nset>;
for n1 in [3..nn] do
 cA:=Flat(car<cA,nset>);
end for;
return cA;
end function;
/////////////////////////////////////////////////////
if tipus eq 1 then

P<[X]>:=PolynomialRing(Rationals(),n+t);
Pol1:=&*[X[a1]-X[a2]: a1,a2 in [1..n]|a1 lt a2];
Pol2:=&*[X[n+a1]-X[n+a2]: a1,a2 in [1..t]|a1 lt a2];
Pol_AB:=(Pol1*Pol2);
R<Y>:=PolynomialRing(P);
ParN:=Partitions(n,t+1);
Rendszer:={};
for k in ParN do
if #PSet eq 0 then
Npar:=particio(k); print Npar;
else
Npar:=PSet;
end if;
if #exptup eq 0 then
LFN:=lf(n);
else
LFN:=exptup;
end if;
 for kk in Npar do
  for l in LFN do
   PRend:={};
   for j in [1..t] do
    PRend:=PRend join {X[n+j]*(&*[(Y-X[m])^l[m]: m in kk[1]])+&*[(Y-X[m])^l[m]: m in kk[j+1]]};
   end for; //j
   PRendseq:=SetToSequence(PRend);
   PR:=&cat[Coefficients(PRendseq[s1]-PRendseq[s2]): s1,s2 in [1..#PRendseq]|s1 lt s2];
   PR_id:=ideal<Parent(PR[1])|PR>;
   PR_GB:=GroebnerBasis(PR_id);
   if not (&or[IsDivisibleBy(Pol_AB,gbpol): gbpol in PR_GB]) then
   fac_GB:=[Factorization(gbpol): gbpol in PR_GB];
   PR_GBreduced:=[&*[cpol[1]: cpol in bpol]: bpol in fac_GB];
   reducedsystem:=[(TT div Gcd(Pol_AB,TT)) div Gcd(TT div Gcd(Pol_AB,TT),Gcd(Pol_AB,TT)): TT in PR_GBreduced];
   if not 1 in reducedsystem then
   xixj:=&or [(X[i]-X[j]) in ideal<P|reducedsystem>: i,j in [1..n]|i ne j];
   yiyj:=&or [(X[i]-X[j]) in ideal<P|reducedsystem>: i,j in [n+1..n+t]|i ne j];
   if (not xixj) and (not yiyj) then Rendszer:=Rendszer join {<kk,l,reducedsystem>}; end if;
   end if;
   end if;
  end for; //lf(n)
 end for; //Npar
end for; //ParN

else

P<[X]>:=PolynomialRing(Rationals(),n+t);
Pol1:=&*[X[a1]-X[a2]: a1,a2 in [1..n]|a1 lt a2];
Pol2:=&*[X[n+a1]-X[n+a2]: a1,a2 in [1..t]|a1 lt a2];
Pol_AB:=(Pol1*Pol2);
R<Y>:=PolynomialRing(P);
ParN:=Partitions(n,t);
Rendszer:={};
for k in ParN do
if #PSet eq 0 then
Npar:=particio(k); print Npar;
else
Npar:=PSet;
end if;
if #exptup eq 0 then
LFN:=lf(n);
else
LFN:=exptup;
end if;
 for kk in Npar do
  for l in LFN do
   PRend:={};
   for j in [[j1,j2,j3]: j1,j2,j3 in [1..t]|j1 lt j2 and j2 lt j3] do
    PRend:=PRend join {(X[n+j[1]]-X[n+j[2]])*(&*[(Y-X[m])^l[m]: m in kk[j[3]]]) + (X[n+j[3]]-X[n+j[1]])*(&*[(Y-X[m])^l[m]: m in kk[j[2]]]) + 
           (X[n+j[2]]-X[n+j[3]])*(&*[(Y-X[m])^l[m]: m in kk[j[1]]])};
   end for; //j
   PRendseq:=SetToSequence(PRend);
   PR:=&cat[Coefficients(PRendseq[s1]): s1 in [1..#PRendseq]];
   if #PR eq 0 then Rendszer:=Rendszer join {<kk,l,[]>};
   else
   PR_id:=ideal<Parent(PR[1])|PR>;
   PR_GB:=GroebnerBasis(PR_id);
   if not (&or[IsDivisibleBy(Pol_AB,gbpol): gbpol in PR_GB]) then
   fac_GB:=[Factorization(gbpol): gbpol in PR_GB];
   PR_GBreduced:=[&*[cpol[1]: cpol in bpol]: bpol in fac_GB];
   reducedsystem:=[(TT div Gcd(Pol_AB,TT)) div Gcd(TT div Gcd(Pol_AB,TT),Gcd(Pol_AB,TT)): TT in PR_GBreduced];
   if not 1 in reducedsystem then
   xixj:=&or [(X[i]-X[j]) in ideal<P|reducedsystem>: i,j in [1..n]|i ne j];
   yiyj:=&or [(X[i]-X[j]) in ideal<P|reducedsystem>: i,j in [n+1..n+t]|i ne j];
   if (not xixj) and (not yiyj) then Rendszer:=Rendszer join {<kk,l,[(TT div Gcd(Pol_AB,TT)) div Gcd(TT div Gcd(Pol_AB,TT),Gcd(Pol_AB,TT)): TT in PR_GBreduced]>}; end if;
   end if;
   end if;
   end if;
  end for; //lf(n)
 end for; //Npar
end for; //ParN

end if;
return Rendszer;
end function;
