######################################################################
##ALG.txt: Save this file as  ALG.txt                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read ALG.txt<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: 

read `AlgDataFile.txt`:
 
print(`Created:  Nov. 2015`):
print(` This is ALG.txt `):
print(`It is one of the packages that accompany the article `):
print(` "Factorization of C-finite Sequences" `):
print(`by Manuel Kauers and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):


print(`---------------------------------------`):
 print(`For a list of the Checking procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Data procedures type ezraD();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezraD:=proc()

if args=NULL then
 print(` The Data procedures are:   DimerData, IsingData  `):
 print(``):

else
ezra(args):
fi:

end:

ezraC:=proc()

if args=NULL then
 print(` The Checking procedures are:  CheckeTp, CheckpTe  `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  Ep, eTp, FactorP1, Kfol1, KfolOld, Khaber1, Khezka1, MulReci, Profile1, Profile11, pTe, TP  `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Cond,   FactorP, Khaber, Khezka,  Mul, ProfileS, TestFact  `):
 print(` `):



elif nops([args])=1 and op(1,[args])=CheckeTp then
print(`CheckeTp(n,N) : checks eTp(e,n,N); Try:`):
print(` CheckeTp(4,6); `):

elif nops([args])=1 and op(1,[args])=CheckpTe then
print(`CheckpTe(n) : checks pTe(p,n); Try:`):
print(` CheckpTe(4); `):

elif nops([args])=1 and op(1,[args])=Cond then
print(`Cond(P,m,n): inputs a symbol P, and positive integers m and n both larger than 1, and outputs the condition(s) for`):
print(`the power sums of a polynomial of degree m*n to be the tensor product of polynomials of degree m and n. Try:`):
print(`Cond(P,2,2);`):

elif nops([args])=1 and op(1,[args])=DimerData then
print(`DimerData(z,x): The denominators of the generating functions for the dimer-problem in an i by infinity strip`):
print(`for i from 1 to 8`):
print(`NOTE: THis procere, due to its size, is in a separate file called AlgDataFile.txt, that is automatically read`):
print(`and should reside in the same directory in your computer`):
print(`Try: `):
print(`DimerData(z,x); ` ):

elif nops([args])=1 and op(1,[args])=Ep then
print(`Ep(p,n,N): expresses the power functions from 1 to N in terms of p_1(x[1], ..., x[n]), ... p_n(x[1], ..., x[n])`):
print(`Try: `):
print(`Ep(p,2,4); `):

elif nops([args])=1 and op(1,[args])=eTp then
print(`eTp(e,n,N): inputs a symbol e and a positive integer n, and outputs polynomial expressions for the power functions`):
print(`p_i(x[1], ..., x[n]) in terms of the elementary symmetric functions e[1], .., e[n], using Newton's relations.`):
print(`for i from 1 to N.`):
print(`Try: `):
print(` eTp(e,3,3); `):


elif nops([args])=1 and op(1,[args])=FactorP then
print(`FactorP(p,x,X): The profile of the tensor-product of p(x) and its inverse x^d*p(1/x).`):
print(`Try: `):
print(`FactorP(1+5*x+101*x^2,x,X); `):

elif nops([args])=1 and op(1,[args])=FactorP1 then
print(`FactorP1(P,x,X): the repetition factorization profile of the polynomial P of x. `):
print(`It first factors P over the complex, and then outputs the sum of X[repeats]`):
print(` Try: `):
print(` FactorP1((x^2+16*x+101)^11*(4*x^3-101)^5,x,X); `):

elif nops([args])=1 and op(1,[args])=FactorP1old then
print(`FactorP1old(P,x,X): the factorization profile of the polynomial P of x. `):
print(`It first factors P over the rationals, and then outputs the set`):
print(` [deg,power] for all the factors that have power higher than 1. Try: `):
print(` FactorP1old((x^2+16*x+101)^11*(4*x^3-101)^5,x,X); `):

elif nops([args])=1 and op(1,[args])=IsingData then
print(`IsingData(z,x): a list of length 4 whose (i-1)-th entry is the characteristic polynomial`):
print(`of the i by infinity transfer, of degree 2^i, in x, for i from 2 to 5. z is the parameter`):
print(`NOTE: THis procere, due to its size, is in a separate file called AlgDataFile.txt, that is automatically read`):
print(`and should reside in the same directory in your computer`):
print(`Try:`):
print(`IsingData(z,x); `):

elif nops([args])=1 and op(1,[args])=KfolOld then
print(`VERY SLOW, DO NOT USE, FOR CHECKING PURPOSES ONLY`):
print(`KfolOld(P1,P2,x): inputs two polynomials, P1, and P2, of the variable x, and outputs a polynomial of degree<= deg(P1)*deg(P2)`):
print(`satisfyied by the product of any root of P1 with and any root of P2. For example, try:`):
print(`KfolOld(x^2+5*x+11,x^2+17*x-13,x);`):
print(``):


elif nops([args])=1 and op(1,[args])=Mul then
print(`Mul(a1,a2,x): the "product" of polynomials a and b of x. The new way, using symmetric functions. Try:`):
print(`Mul((x-2)*(x-3),(x-5)*(x-11),x);`):

elif nops([args])=1 and op(1,[args])=Kfol1 then
print(`Kfol1(P1,P2,x,A): inputs two polynomials, P1, and P2, of the variable x, and a positive integer A.`):
print(` Outputs a polynomial of degree A, if it exists, or FAIL, `):
print(`satisfyied by the product of any root of P1 with and any root of P2. For example, try:`):
print(`Kfol1(x^2+5*x+11,x^2+17*x-13,x,4);`):
print(``):

elif nops([args])=1 and op(1,[args])=Khaber then
print(`Khaber(P1,P2,x): inputs two polynomials, P1, and P2, of the variable x, and outputs a polynomial of degree <=deg(P1)*deg(P2)`):
print(`satisfyied by the sum of any root of P1 with and any root of P2. For example, try:`):
print(`Khaber(x^2+5*x+11,x^2+17*x-13,x);`):
print(``):

elif nops([args])=1 and op(1,[args])=Khaber1 then
print(`Khaber1(P1,P2,x,A): inputs two polynomials, P1, and P2, of the variable x, and a positive integer A. Outputs a polynomial of degree A`):
print(`satisfyied by the sum of any root of P1 with and any root of P2, or FAIL. For example, try:`):
print(`Khaber1(x^2+5*x+11,x^2+17*x-13,x,4);`):
print(``):


elif nops([args])=1 and op(1,[args])=Khezka then
print(`Khezka(P1,x,p): inputs a polynomial, P1, of the variable x, and a positive integer p,`):
print(`and outputs a polynomial `):
print(` satisfyied by the p^th power of and any root of P1. For example, try: `):
print(` Khezka(x^2+5*x+11,x,2); `):

elif nops([args])=1 and op(1,[args])=Khezka1 then
print(`Khezka1(P1,x,p,A): inputs a polynomial, P1, of the variable x, and a positive integer p,`):
print(`and outputs a polynomial of degree A`):
print(` satisfyied by the p^th power of and any root of P1. For example, try: `):
print(` Khezka1(x^2+5*x+11,x,2,2); `):


elif nops([args])=1 and op(1,[args])=MulReci then
print(`MulReci(p,x): inputs a polynomial p in the variable x, of degree d, say, outputs the polynomial of degree d^2-d that is`):
print(` Mul(p,p(1/x)*x^d)/(1-x)^d. Try: `):
print(` MulReci(x^2+11*x+5,x); `):

elif nops([args])=1 and op(1,[args])=pTe then
print(`pTe(p,n): inputs a symbol p and a positive integer n, and outputs polynomial expressions for the elementary symmetric functions`):
print(`e_i(x[1], ..., x[n]) in terms of the power functions p[1], .., p[n], using Newton's relations.`):
print(`for i from 1 to n.`):
print(`Try:`):
print(`pTe(p,3);`):

elif nops([args])=1 and op(1,[args])=Profile1 then
print(`Profile1(L,X,k): inputs a list of positive integers (larger than 1),L, of size k, say, outputs the profile`):
print(`of the tensor product of the tensor-product p[1]x...xp[k] where p[i] is a generic polynomial of degree L[i].`):
print(` It does it for k random choices. If they don't all agree, it returns FAIL `):
print(`Warning: very slow, for testing only. Please use ProfileS(L,X)`):
print(` Try: `):
print(` Profile1([2,2],X,5); `):
print(` Profile1([2,4],X,3): `):


elif nops([args])=1 and op(1,[args])=Profile11 then
print(`Profile11(L,X): inputs a list of positive integers (larger than 1),L, of size k, say, outputs the profile`):
print(`of the tensor product of the tensor-product p[1]x...xp[k] where p[i] is a generic polynomial of degree L[i].`):
print(` It does it for one random choice `):
print(` Try: `):
print(` Profile11([2,2],X); `):
print(` Profile11([2,4],X); `):


elif nops([args])=1 and op(1,[args])=ProfileS then
print(`ProfileS(L,X): inputs a list of positive integers (larger than 1),L, of size k, say, outputs the profile`):
print(`of the tensor product of the tensor-product p[1]x...xp[k] where p[i] is a generic polynomial of degree L[i].`):
print(` It is fully rigorous.` ):
print(` ProfileS([2,2],X); `):
print(` ProfileS([2,4],X): `):

elif nops([args])=1 and op(1,[args])=TestFact then

print(`TestFact(p,x,L): inputs a polynomial p of degree mul(L[i],i=1..nops(L)) of x, and outputs true iff it is  a tensor product of`):
print(` a polynomials of degrees L[1], .., L[nops(L)]. Try:`):
print(`TestFact(Mul(101+24*x+511*x^2,41+162*x+153*x^2,x),x,[2,2]);`):

elif nops([args])=1 and op(1,[args])=TestFactOld then
print(`TestFactOld(p,x,m,n): inputs a polynomial p of degree m*n of x, and outputs true iff it is  a tensor product of`):
print(` a polynomial of degree m and a polynomial of degree n `):
print(` TestFactOld(Mul(1+4*x+5*x^2,4+16*x+15*x^2,x),x,2,2); `):


elif nops([args])=1 and op(1,[args])=TP then
print(`TP(L1,L2): The tensor product of lists L1 and L2`):
print(`Try:`):
print(`TP([a1,a2],[b1,b2]);`):


else
print(`There is no ezra for`,args):
fi:
 
end:



###start data



#DimerData(z,x): The denominators of the generating functions for the dimer-problem in an i by infinity strip
#for i from 1 to 8
DimerData:=proc(z,x)
[x^2-1, x^2-1-z*x, x^4+1-(2*z^2+2)*x^2, x^4+1-z^2*x^3-(3*z^2+2)*x^2-z^2*x, x^8+
1-(3*z^4+8*z^2+4)*x^6-(-10*z^4-16*z^2-6)*x^4-(3*z^4+8*z^2+4)*x^2, x^8+1-z^3*x^7
-(6*z^4+10*z^2+4)*x^6-(5*z^5+5*z^3)*x^5-(z^6-13*z^4-20*z^2-6)*x^4-(-5*z^5-5*z^3
)*x^3-(6*z^4+10*z^2+4)*x^2+z^3*x, x^16+1-(4*z^6+20*z^4+24*z^2+8)*x^14-(-52*z^8-\
192*z^6-256*z^4-144*z^2-28)*x^12-(64*z^10+416*z^8+892*z^6+844*z^4+360*z^2+56)*x
^10-(-16*z^12-160*z^10-744*z^8-1408*z^6-1216*z^4-480*z^2-70)*x^8-(64*z^10+416*z
^8+892*z^6+844*z^4+360*z^2+56)*x^6-(-52*z^8-192*z^6-256*z^4-144*z^2-28)*x^4-(4*
z^6+20*z^4+24*z^2+8)*x^2, x^16+1-z^4*x^15-(10*z^6+30*z^4+28*z^2+8)*x^14-(15*z^8
+35*z^6+19*z^4)*x^13-(7*z^10-78*z^8-300*z^6-354*z^4-168*z^2-28)*x^12-(z^12-75*z
^10-230*z^8-217*z^6-63*z^4)*x^11-(-15*z^12+148*z^10+822*z^8+1442*z^6+1146*z^4+
420*z^2+56)*x^10-(69*z^12+232*z^10+303*z^8+182*z^6+43*z^4)*x^9-(-119*z^12-642*z
^10-1673*z^8-2304*z^6-1644*z^4-560*z^2-70)*x^8-(69*z^12+232*z^10+303*z^8+182*z^
6+43*z^4)*x^7-(-15*z^12+148*z^10+822*z^8+1442*z^6+1146*z^4+420*z^2+56)*x^6-(z^
12-75*z^10-230*z^8-217*z^6-63*z^4)*x^5-(7*z^10-78*z^8-300*z^6-354*z^4-168*z^2-\
28)*x^4-(15*z^8+35*z^6+19*z^4)*x^3-(10*z^6+30*z^4+28*z^2+8)*x^2-z^4*x, x^32+1-(
5*z^8+40*z^6+84*z^4+64*z^2+16)*x^30-(-190*z^12-1200*z^10-3026*z^8-3872*z^6-2632
*z^4-896*z^2-120)*x^28-(655*z^16+7400*z^14+30734*z^12+65392*z^10+80039*z^8+
58488*z^6+25116*z^4+5824*z^2+560)*x^26-(-550*z^20-9600*z^18-75506*z^16-291824*z
^14-635896*z^12-847104*z^10-715572*z^8-384192*z^6-126672*z^4-23296*z^2-1820)*x^
24-(125*z^24+3000*z^22+41860*z^20+312560*z^18+1269923*z^16+3077832*z^14+4746678
*z^12+4822192*z^10+3260653*z^8+1448360*z^6+404404*z^4+64064*z^2+4368)*x^22-(-\
3275*z^24-63200*z^22-546320*z^20-2580000*z^18-7423128*z^16-13887552*z^14-\
17518258*z^12-15134672*z^10-8937678*z^8-3532000*z^6-888888*z^4-128128*z^2-8008)
*x^20-(20775*z^24+326600*z^22+2165340*z^20+8191760*z^18+19832526*z^16+32436304*
z^14+36765756*z^12+29101024*z^10+15961703*z^8+5915064*z^6+1405404*z^4+192192*z^
2+11440)*x^18-(-41650*z^24-558400*z^22-3359060*z^20-11855040*z^18-27215340*z^16
-42684320*z^14-46777648*z^12-36011264*z^10-19292248*z^8-7003776*z^6-1633632*z^4
-219648*z^2-12870)*x^16-(20775*z^24+326600*z^22+2165340*z^20+8191760*z^18+
19832526*z^16+32436304*z^14+36765756*z^12+29101024*z^10+15961703*z^8+5915064*z^
6+1405404*z^4+192192*z^2+11440)*x^14-(-3275*z^24-63200*z^22-546320*z^20-2580000
*z^18-7423128*z^16-13887552*z^14-17518258*z^12-15134672*z^10-8937678*z^8-\
3532000*z^6-888888*z^4-128128*z^2-8008)*x^12-(125*z^24+3000*z^22+41860*z^20+
312560*z^18+1269923*z^16+3077832*z^14+4746678*z^12+4822192*z^10+3260653*z^8+
1448360*z^6+404404*z^4+64064*z^2+4368)*x^10-(-550*z^20-9600*z^18-75506*z^16-\
291824*z^14-635896*z^12-847104*z^10-715572*z^8-384192*z^6-126672*z^4-23296*z^2-\
1820)*x^8-(655*z^16+7400*z^14+30734*z^12+65392*z^10+80039*z^8+58488*z^6+25116*z
^4+5824*z^2+560)*x^6-(-190*z^12-1200*z^10-3026*z^8-3872*z^6-2632*z^4-896*z^2-\
120)*x^4-(5*z^8+40*z^6+84*z^4+64*z^2+16)*x^2, 1-z^5*x^31-(15*z^8+70*z^6+112*z^4
+72*z^2+16)*x^30-(35*z^11+140*z^9+171*z^7+65*z^5)*x^29-(28*z^14-322*z^12-2265*z
^10-5184*z^8-5768*z^6-3388*z^4-1008*z^2-120)*x^28-(9*z^17-566*z^15-3215*z^13-\
6692*z^11-6587*z^9-3087*z^7-551*z^5)*x^27-(z^20-255*z^18+1353*z^16+19012*z^14+
68971*z^12+125890*z^10+133221*z^8+84938*z^6+32032*z^4+6552*z^2+560)*x^26-(-35*z
^21+1984*z^19+15158*z^17+46708*z^15+77312*z^13+74348*z^11+41517*z^9+12474*z^7+
1561*z^5)*x^25-(411*z^22-3406*z^20-55574*z^18-291977*z^16-829849*z^14-1442558*z
^12-1611098*z^10-1174278*z^8-552608*z^6-160888*z^4-26208*z^2-1820)*x^24-(-2164*
z^23-15597*z^21-52968*z^19-105361*z^17-129494*z^15-99521*z^13-47920*z^11-14785*
z^9-3114*z^7-395*z^5)*x^23-(5527*z^24+61547*z^22+386154*z^20+1608919*z^18+
4462986*z^16+8297046*z^14+10472463*z^12+9038642*z^10+5307687*z^8+2073470*z^6+
512512*z^4+72072*z^2+4368)*x^22-(-6907*z^25-62987*z^23-295050*z^21-909216*z^19-\
1898966*z^17-2657435*z^15-2455558*z^13-1463330*z^11-535745*z^9-108423*z^7-9163*
z^5)*x^21-(4508*z^26-4591*z^24-378657*z^22-2676224*z^20-9843835*z^18-22728088*z
^16-35204576*z^14-37654642*z^12-28066567*z^10-14480232*z^8-5043640*z^6-1125124*
z^4-144144*z^2-8008)*x^20-(-1591*z^27+41809*z^25+525554*z^23+2582552*z^21+
7213271*z^19+12799308*z^17+15047773*z^15+11830501*z^13+6121965*z^11+1989750*z^9
+365715*z^7+28765*z^5)*x^19-(300*z^28-25078*z^26-95157*z^24+783883*z^22+7101270
*z^20+25736222*z^18+55416381*z^16+78964798*z^14+77734646*z^12+53630348*z^10+
25795557*z^8+8435826*z^6+1777776*z^4+216216*z^2+11440)*x^18+x^32-(-28*z^29+6237
*z^27-90066*z^25-1297784*z^23-6358401*z^21-17032260*z^19-28498360*z^17-31413624
*z^15-23177509*z^13-11309780*z^11-3490046*z^9-613764*z^7-46563*z^5)*x^17-(z^30-\
694*z^28+42087*z^26+192164*z^24-927407*z^22-9606494*z^20-34884094*z^18-73719263
*z^16-102492918*z^14-98357116*z^12-66229900*z^10-31153572*z^8-9984576*z^6-\
2066064*z^4-247104*z^2-12870)*x^16-(28*z^29-6237*z^27+90066*z^25+1297784*z^23+
6358401*z^21+17032260*z^19+28498360*z^17+31413624*z^15+23177509*z^13+11309780*z
^11+3490046*z^9+613764*z^7+46563*z^5)*x^15-(300*z^28-25078*z^26-95157*z^24+
783883*z^22+7101270*z^20+25736222*z^18+55416381*z^16+78964798*z^14+77734646*z^
12+53630348*z^10+25795557*z^8+8435826*z^6+1777776*z^4+216216*z^2+11440)*x^14-(
1591*z^27-41809*z^25-525554*z^23-2582552*z^21-7213271*z^19-12799308*z^17-\
15047773*z^15-11830501*z^13-6121965*z^11-1989750*z^9-365715*z^7-28765*z^5)*x^13
-(4508*z^26-4591*z^24-378657*z^22-2676224*z^20-9843835*z^18-22728088*z^16-\
35204576*z^14-37654642*z^12-28066567*z^10-14480232*z^8-5043640*z^6-1125124*z^4-\
144144*z^2-8008)*x^12-(6907*z^25+62987*z^23+295050*z^21+909216*z^19+1898966*z^
17+2657435*z^15+2455558*z^13+1463330*z^11+535745*z^9+108423*z^7+9163*z^5)*x^11-
(5527*z^24+61547*z^22+386154*z^20+1608919*z^18+4462986*z^16+8297046*z^14+
10472463*z^12+9038642*z^10+5307687*z^8+2073470*z^6+512512*z^4+72072*z^2+4368)*x
^10-(2164*z^23+15597*z^21+52968*z^19+105361*z^17+129494*z^15+99521*z^13+47920*z
^11+14785*z^9+3114*z^7+395*z^5)*x^9-(411*z^22-3406*z^20-55574*z^18-291977*z^16-\
829849*z^14-1442558*z^12-1611098*z^10-1174278*z^8-552608*z^6-160888*z^4-26208*z
^2-1820)*x^8-(35*z^21-1984*z^19-15158*z^17-46708*z^15-77312*z^13-74348*z^11-\
41517*z^9-12474*z^7-1561*z^5)*x^7-(z^20-255*z^18+1353*z^16+19012*z^14+68971*z^
12+125890*z^10+133221*z^8+84938*z^6+32032*z^4+6552*z^2+560)*x^6-(-9*z^17+566*z^
15+3215*z^13+6692*z^11+6587*z^9+3087*z^7+551*z^5)*x^5-(28*z^14-322*z^12-2265*z^
10-5184*z^8-5768*z^6-3388*z^4-1008*z^2-120)*x^4-(-35*z^11-140*z^9-171*z^7-65*z^
5)*x^3-(15*z^8+70*z^6+112*z^4+72*z^2+16)*x^2+z^5*x]:
end:


###end data
#Khaber1(P1,P2,x,A): inputs two polynomials, P1, and P2, of the variable x, and a positive integer A,
#and outputs a polynomial of degree A
#satisfyied by the sum of any root of P1 with and any root of P2. For example, try:
#Khaber1(x^2+5*x+11,x^2+17*x-13,x,4);
Khaber1:=proc(P1,P2,x,A) local d1,d2,i,gu,eq,var,a,T1,T2,x1,x2,Pa,P,i1,i2:

d1:=degree(P1,x):

d2:=degree(P2,x):


Pa:=-(P1-coeff(P1,x,d1)*x^d1)/coeff(P1,x,d1):
T1[d1]:=-(P1-coeff(P1,x,d1)*x^d1)/coeff(P1,x,d1):

for i from d1+1 to A do
T1[i]:=expand(T1[i-1]*x):
T1[i]:=expand(subs(x^d1=Pa,T1[i])):
od:


Pa:=-(P2-coeff(P2,x,d2)*x^d2)/coeff(P2,x,d2):
T2[d2]:=-(P2-coeff(P2,x,d2)*x^d2)/coeff(P2,x,d2):

for i from d2+1 to A do
T2[i]:=expand(T2[i-1]*x):
T2[i]:=expand(subs(x^d2=Pa,T2[i])):
od:

P:=x^A+add(a[i]*x^i,i=0..A-1):
var:={seq(a[i],i=0..A-1)}:

gu:=expand(subs(x=x1+x2,P)):

for i1 from d1 to A do
 gu:=expand(subs(x1^i1=subs(x=x1,T1[i1]),gu)):
od:

for i2 from d2 to A do
 gu:=expand(subs(x2^i2=subs(x=x2,T2[i2]),gu)):
od:

eq:={seq(seq(coeff(coeff(gu,x1,i1),x2,i2),i1=0..d1-1),i2=0..d2-1)}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
fi:

subs(var,P):

end:

#Khaber(P1,P2,x): inputs two polynomials, P1, and P2, of the variable x, and outputs a polynomial of degree deg(P1)*deg(P2)
#satisfyied by the sum of any root of P1 with and any root of P2. For example, try:
#Khaber(x^2+5*x+11,x^2+17*x-13,x);
Khaber:=proc(P1,P2,x) local gu,i:

for i from 1 to degree(P1,x)*degree(P2,x) do
 gu:=Khaber1(P1,P2,x,i):
  if gu<>FAIL then
    RETURN(gu):
  fi:
od:
FAIL:

end:




#Kfol1(P1,P2,x,A): inputs two polynomials, P1, and P2, of the variable x, and outputs a polynomial of degree A
#satisfyied by the product of a root of P1 with and any root of P2. For example, try:
#Kfol1(x^2+5*x+11,x^2+17*x-13,x);
Kfol1:=proc(P1,P2,x,A) local d1,d2,i,gu,eq,var,a,T1,T2,x1,x2,Pa,P,i1,i2:

d1:=degree(P1,x):

d2:=degree(P2,x):


Pa:=-(P1-coeff(P1,x,d1)*x^d1)/coeff(P1,x,d1):
T1[d1]:=-(P1-coeff(P1,x,d1)*x^d1)/coeff(P1,x,d1):

for i from d1+1 to A do
T1[i]:=expand(T1[i-1]*x):
T1[i]:=expand(subs(x^d1=Pa,T1[i])):
od:


Pa:=-(P2-coeff(P2,x,d2)*x^d2)/coeff(P2,x,d2):
T2[d2]:=-(P2-coeff(P2,x,d2)*x^d2)/coeff(P2,x,d2):

for i from d2+1 to A do
T2[i]:=expand(T2[i-1]*x):
T2[i]:=expand(subs(x^d2=Pa,T2[i])):
od:

P:=x^A+add(a[i]*x^i,i=0..A-1):
var:={seq(a[i],i=0..A-1)}:

gu:=expand(subs(x=x1*x2,P)):

for i1 from d1 to A do
 gu:=expand(subs(x1^i1=subs(x=x1,T1[i1]),gu)):
od:

for i2 from d2 to A do
 gu:=expand(subs(x2^i2=subs(x=x2,T2[i2]),gu)):
od:

eq:={seq(seq(coeff(coeff(gu,x1,i1),x2,i2),i1=0..d1-1),i2=0..d2-1)}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
fi:

subs(var,P):

end:


#KfolOld(P1,P2,x): inputs two polynomials, P1, and P2, of the variable x, and outputs a polynomial 
#satisfyied by the sum of any root of P1 with and any root of P2. For example, try:
#KfolOld(x^2+5*x+11,x^2+17*x-13,x);
KfolOld:=proc(P1,P2,x) local gu,i:

for i from 1 to degree(P1,x)*degree(P2,x) do
 gu:=Kfol1(P1,P2,x,i):
  if gu<>FAIL then
    RETURN(sort(gu)):
  fi:
od:
FAIL:

end:



#Khezka1(P1,x,p,A): inputs a polynomial, P1, of the variable x, and a positive integer p,
#and outputs a polynomial of degree A
#satisfyied by the p^th power of and any root of P1. For example, try:
#Khezka1(x^2+5*x+11,x,2,2);
Khezka1:=proc(P1,x,p,A) local d1,i,gu,eq,var,a,T1,x1,Pa,P,i1:

d1:=degree(P1,x):

T1[0]:=1:

Pa:=-(P1-coeff(P1,x,d1)*x^d1)/coeff(P1,x,d1):
T1[d1]:=-(P1-coeff(P1,x,d1)*x^d1)/coeff(P1,x,d1):

for i from d1+1 to A*p do
T1[i]:=expand(T1[i-1]*x):
T1[i]:=expand(subs(x^d1=Pa,T1[i])):
od:

P:=x^A+add(a[i]*x^i,i=0..A-1):
var:={seq(a[i],i=0..A-1)}:

gu:=expand(subs(x=x1^p,P)):

for i1 from 0 to A do
 gu:=expand(subs(x1^(i1*p)=subs(x=x1,T1[i1*p]),gu)):
od:

eq:={seq(coeff(gu,x1,i1),i1=0..d1-1)}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
fi:

subs(var,P):

end:

#Khezka(P1,x,p): inputs a polynomial, P1 of the variable x, and a positive integer p
#outputs a polynomial
#satisfyied by the p-th powers of the roots P1 . Try:
#Khezka(x^2+5*x+11,x^2+17*x-13,x,3);
Khezka:=proc(P1,x,p) local i,gu:

for i from 1 to degree(P1,x) do
 gu:=Khezka1(P1,x,p,i):
  if gu<>FAIL then
    RETURN(gu):
  fi:
od:
FAIL:

end:




#eTp(e,n,N): inputs a symbol e and a positive integer n, and outputs polynomial expressions for the power functions
#p_i(x[1], ..., x[n]) in terms of the elementary symmetric functions e[1], .., e[n], using Newton's relations.
#for i from 1 to N.
#Try:
#eTp(e,3,3);
eTp:=proc(e,n,N) local gu,k,gu1,r:
gu:=[e[1]]:

if N<=n then
for k from 2 to n do
 gu1:=expand(-((-1)^k*k*e[k]+add((-1)^r*e[r]*gu[k-r],r=1..k-1))):
 gu:=[op(gu),gu1]:
od:

RETURN(gu):
fi:

if N>n then
for k from 2 to n do
 gu1:=expand(-((-1)^k*k*e[k]+add((-1)^r*e[r]*gu[k-r],r=1..k-1))):
 gu:=[op(gu),gu1]:
od:

for k from n+1 to N do
 gu1:=expand(add((-1)^(r-1)*e[r]*gu[k-r],r=1..n)):
 gu:=[op(gu),gu1]:
od:

RETURN(gu):
fi:

end:

#CheckeTp(n,N) : checks eTp(e,n,N); Try:
#CheckeTp(4,6);
CheckeTp:=proc(n,N) local e,i,P,t,gu,x,j:
P:=expand(mul(1+t*x[i],i=1..n)):
gu:=eTp(e,n,N):
gu:=expand(subs({seq(e[i]=coeff(P,t,i),i=1..n)},gu)):
evalb(expand({seq(gu[i]-add(x[j]^i,j=1..n),i=1..nops(gu))})={0}):
end:



#pTe(p,n): inputs a symbol p and a positive integer n, and outputs polynomial expressions for the elementary symmetric functions
#e_i(x[1], ..., x[n]) in terms of the power functions p[1], .., p[n], using Newton's relations.
#for i from 1 to n.
#Try:
#pTe(p,3);
pTe:=proc(p,n) local gu,k,gu1,r:
gu:=[p[1]]:

for k from 2 to n do
 gu1:=expand((-1)^(k-1)*(p[k]+add((-1)^r*gu[r]*p[k-r],r=1..k-1))/k):
 gu:=[op(gu),gu1]:
od:

gu:



end:

#CheckpTe(n) : checks pTe(p,n); Try:
#CheckpTe(4);
CheckpTe:=proc(n) local p,i,P,t,gu,x,j:
P:=expand(mul(1+t*x[i],i=1..n)):
gu:=pTe(p,n):
gu:=expand(subs({seq(p[i]=add(x[j]^i,j=1..n),i=1..n)},gu)):
evalb(expand({seq(gu[i]-coeff(P,t,i),i=1..n)})={0}):
end:

#Ep(p,n,N): expresses the power functions from 1 to N in terms of p_1(x[1], ..., x[n]), ... p_n(x[1], ..., x[n])
#Try:
#Ep(p,2,4);
Ep:=proc(p,n,N) local gu,mu,e,j,i:
if N<=n then
 RETURN([seq(p[i],i=1..N)]):
fi:
mu:=pTe(p,n):
gu:=eTp(e,n,N):

[seq(p[i],i=1..n), seq(expand(subs({seq(e[j]=mu[j],j=1..n)},gu[i])),i=n+1..N)]:
end:


#Cond(P,m,n): inputs a symbol P, and positive integers m and n both larger than 1, and outputs the condition(s) for
#the power sums of a polynomial of degree m*n to be the tensor product of polynomials of degree m and n. Try:
#Cond(P,2,2);
Cond:=proc(P,m,n) local gu1,gu2,gu,p,q,i:
 gu1:=Ep(p,m,m*n):
 gu2:=Ep(q,n,m*n):

gu:=[seq(P[i]-gu1[i]*gu2[i],i=1..m*n)]:

Groebner[Basis](gu,plex(seq(p[i],i=1..m),seq(q[i],i=1..n), seq(P[i],i=1..m*n)));
end:





#Mul(a1,a2,x): the "product" of polynomials a and b of x. Try:
#Mul((x-2)*(x-3),(x-5)*(x-11),x);
Mul:=proc(a1,a2,x) local a11,a21,d1,d2,i,e1,e2,p1,p2,P,EL:
d1:=degree(a1,x):
d2:=degree(a2,x):

a11:=a1/coeff(a1,x,d1):
a21:=a2/coeff(a2,x,d2):
e1:=[seq((-1)^i*coeff(a11,x,d1-i),i=1..d1)]:
e2:=[seq((-1)^i*coeff(a21,x,d2-i),i=1..d2)]:
p1:=eTp(e1,d1,d1*d2):
p2:=eTp(e2,d2,d1*d2):
P:=[seq(p1[i]*p2[i],i=1..d1*d2)]:
EL:=pTe(P,d1*d2):
sort(x^(d1*d2)+add((-1)^i*EL[i]*x^(d1*d2-i),i=1..d1*d2)):
end:

#MulReci(p,x): inputs a polynomial p in the variable x, of degree d, say, outputs the polynomial of degree d^2 that is
#Mul(p,p(1/x)*x^d). Try:
#MulReci(x^2+11*x+5,x);
MulReci:=proc(p,x) local d,p1:
d:=degree(p,x):
p1:=normal(x^d*subs(x=1/x,p)):
factor(normal(Mul(p,p1,x)/(1-x)^d)):
end:

#TestFactOld(p,x,m,n): inputs a polynomial p of degree m*n of x, and outputs true iff it is  a tensor product of
#a polynomial of degree m and a polynomial of degree n
#TestFactOld(Mul(1+4*x+5*x^2,4+16*x+15*x^2,x),x,2,2);
TestFactOld:=proc(p,x,m,n) local gu,lu1,lu2,gu1,i:

if degree(p,x)<>m*n then
 print(p, `should have been of degree`, m*n):
 RETURN(FAIL):
fi:

gu:=factor(normal(MulReci(p,x)/(x-1)^(m*n))):


if m<>n then

lu1:=1:
lu2:=1:

for i from 1 to nops(gu) do

gu1:=op(i,gu):


if type(gu1,`^`) then
 if op(2,gu1)=n then
   lu1:=lu1*op(1,gu1):
 elif op(2,gu1)=m then
   lu2:=lu2*op(1,gu1):
fi:
fi:
od:

if degree(lu1,x)=m*(m-1) and degree(lu2,x)=n*(n-1) then
  RETURN(true):
else
 RETURN(false):
fi:


else

lu1:=1:
lu2:=1:

for i from 1 to nops(gu) do

gu1:=op(i,gu):


if type(gu1,`^`) then
 if op(2,gu1)=n then
   lu1:=lu1*op(1,gu1):
 fi:
fi:

od:

if degree(lu1,x)=2*m*(m-1) then
  RETURN(true):
else
 RETURN(false):
fi:



fi:



end:



 

#FactorP1old(P,x,X): the factorization profile of the polynomial P of x. 
#It first factors P over the rationals, and then outputs the set
#[deg,power] for all the factors that have power higher than 1. Try:
#FactorP1old((11*x^3+101)^11*(x^2+6*x+1)^11,x,X);
FactorP1old:=proc(P,x,X) local P1,P1a, gu,i:
P1:=factor(P):

gu:=[]:

for i from 1 to nops(P1) do
 P1a:=op(i,P1):
 if type(P1a,`^`) then
  gu:=[op(gu),[degree(op(1,P1a),x),op(2,P1a)]]:
 fi:
od:


add(X[gu[i]],i=1..nops(gu)):


end:


#FactorP1(P,x,X): the factorization profile of the polynomial P of x. 
#It first factors P over the complex, and then outputs the set
#of all repetitins  Try:
#FactorP1((11*x^3+101)^11*(x^2+6*x+1)^11,x,X);
FactorP1:=proc(P,x,X) local P1,P1a, gu,i:
P1:=factor(P):

if type(P1,`+`) then
 RETURN(0):
fi:

if type(P1,`^`) then
 RETURN([op(2,P1)$degree(op(1,P1),x) ]):
fi:
gu:=[]:

for i from 1 to nops(P1) do
 P1a:=op(i,P1):
 if type(P1a,`^`) then
  gu:=[op(gu),op(2,P1a)$degree(op(1,P1a),x) ]:
 fi:
od:


add(X[gu[i]],i=1..nops(gu)):


end:




#Profile11(L,X): inputs a list of positive integers (larger than 1),L, of size k, say, outputs the profile
#of the tensor product of the tensor-product p[1]x...xp[k] where p[i] is a generic polynomial of degree L[i].
#It does it for one random choice
#Try:
#Profile11([2,2],X);
#Profile11([2,4],X):
Profile11:=proc(L,X) local L1,ra,x,i,i1,r,p,gu:

L1:=sort(L):
ra:=rand(10..1000):

r:=add(ra()*x^i1,i1=0..L1[1]):
for i from 2 to nops(L1) do
 p:=add(ra()*x^i1,i1=0..L1[i]):
 r:=Mul(r,p,x):
od:

gu:=FactorP1(MulReci(r,x),x,X);


end:

#Profile1(L,X,k): inputs a list of positive integers (larger than 1),L, of size k, say, outputs the profile
#of the tensor product of the tensor-product p[1]x...xp[k] where p[i] is a generic polynomial of degree L[i].
#It does it for k random choice, and returns FAIL if they disagree
#Profile11([2,2],X,10);
#Profile11([2,4],X,3);
Profile1:=proc(L,X,k) local gu,i:
option remember:
gu:={Profile11(L,X)}:

for i from 2 to k do
 gu:=gu union {Profile11(L,X)}:
od:

if nops(gu)>1 then
 RETURN(FAIL):
else
 RETURN(gu[1]):
fi:

end:




#TestFact(p,x,L): inputs a polynomial p of degree mul(L[i],i=1..nops(L)) of x, and outputs true iff it is  a tensor product of
#a polynomials of degrees L[1], .., L[nops(L)]. Try:
#TestFact(Mul(1+4*x+5*x^2,4+16*x+15*x^2,x),x,[2,2]);
TestFact:=proc(p,x,L) local Pr,gu,X,i1:

if degree(p,x)<>mul(L[i1],i1=1..nops(L)) then
 print(`the degree of the input polynomial should be`,  mul(L[i1],i1=1..nops(L)) ):
 RETURN(FAIL):
fi:

Pr:=ProfileS(L,X):

if Pr=FAIL then
 RETURN(FAIL):
fi:

gu:=FactorP1(MulReci(p,x),x,X):

evalb(gu=Pr):

end:


#FactorP(p,x,X): The profile of the tensor-product of p(x) and its inverse x^d*p(1/x).
#Try:
#FactorP(1+5*x+101*x^2,x,X)
FactorP:=proc(p,x,X)
FactorP1(MulReci(p,x),x,X):
end:




#TP(L1,L2): The tensor product of lists L1 and L2
#Try:
#TP([a1,a2],[b1,b2]);
TP:=proc(L1,L2) local i,j:
[seq(seq(L1[i]*L2[j],j=1..nops(L2)),i=1..nops(L1))]:
end:



#ProfileS(L,X): inputs a list of positive integers (larger than 1),L, of size k, say, outputs the profile
#of the tensor product of the tensor-product p[1]x...xp[k] where p[i] is a generic polynomial of degree L[i].
#it is fully rigorous. Try:
#ProfileS([2,2],X);
#ProfileS([2,2,2],X);
ProfileS:=proc(L,X) local gu,a,i,j,mu,lu,T:
gu:=[seq(a[1,j],j=1..L[1])]:

for i from 2 to nops(L) do
 gu:=TP(gu,[seq(a[i,j],j=1..L[i])]):
od:

mu:=TP(gu,[seq(1/gu[i],i=1..nops(gu))]):

lu:=convert(mu,set):

for i from 1 to nops(lu) do
T[lu[i]]:=0:
od:

for i from 1 to nops(mu) do
T[mu[i]]:=T[mu[i]]+1:
od:

gu:=add(X[T[lu[i]]],i=1..nops(lu)) -X[convert(L,`*`)]:
gu:=gu-coeff(gu,X[1],1)*X[1]:
gu:
end: