##########################################################################
## RamanujanCubic.txt: Save this file as   RamanujanCubic.txt to use it, #
# stay in the                                                            #
## same directory, get into Maple (by typing: maple <Enter> )            #
## and then type:  read  `RamanujanCubic.txt`:  <Enter>                  #	
## Then follow the instructions given there                              #
##                                                                       #
## Written by Doron Zeilberger, Rutgers University ,                     #
##  DoronZeil at gmail dot com                                           # 
#########################################################################


print(`First Written: July 2020: tested for Maple 2018 `):
print(`Version  : July 2020:  `):
print():
print(`This is  RamanujanCubic.txt A Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(` Automatic Solving of Cubic Diophantine Equations Inspired by Ramanujan`):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/RamanujanCubic.txt `):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():

print(`-------------------------------`):
print(`For general help, and a list of the MAIN functions,`):

print(` type "ezra();". For specific help type "ezra(procedure_name);" `):

print(`-------------------------------`):


print(`-------------------------------`):
print(`For general help, and a list of the STORY functions,`):

print(` type "ezraSt();". For specific help type "ezra(procedure_name);" `):

print(`-------------------------------`):

print(`For a list of the supporting functions type: ezra1();`):
print(` For specific help type "ezra(procedure_name);" `):

print(`-------------------------------`):


print(`-------------------------------`):

print(`For a list of the  functions from RgfDio.txt type: ezraR();`):
print(` For specific help type "ezraR(procedure_name);" `):

print(`-------------------------------`):


with(combinat):



ezraSt:=proc()
if args=NULL then

print(`The Story procedures are`):
print(`  Paper, PaperV, TheoremABs, TheoremABsV `):

else
ezra(args):

fi:


end:

ezraOld:=proc()
if args=NULL then

print(`The Depretated procedures are`):
print(`  ExtractCoeffsOld, FindEqsPold, SolDEPold `):

else
ezra(args):

fi:


end:


ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(`  AllSols, AllSolsP, BFd1, Box1, Cands, CleanUp, ExtractCoeffs, Find43, FindEqsP,  FindGoodR, GenPol1, IsTrivial, Katan, MakeRama1, `):
print(`MakeRamaSol, RandPolH,  `):
print(` Reject, Rootk, `):
 print(` SolBF, SolBF1, SolBF1a,  , SolDEP, SortLE  `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then


print(`The MAIN procedures are:`):
print(` BFd, BFdG, EriJ, EriJs,  FindDE, MakeRama, Morph43, ParAB3, TheoremAB  `):



elif nargs=1 and args[1]=AllSols then
print(`AllSols(P,m,n,K): inputs a polynomial P in m and n and a positive integer K outputs all`):
print(`solutions in -K<=m,n<=K of P(m,n)=0. Try:`):
print(`AllSols(m^2-2*n^2-1,m,n,30);`):


elif nargs=1 and args[1]=AllSolsP then
print(`AllSolsP(P,m,n,K): inputs a polynomial P in m and n and a positive integer K outputs all`):
print(`solutions in 0<=m,n<=K of P(m,n)=0. Try:`):
print(`AllSolsP(m^2-2*n^2-1,m,n,30);`):

elif nargs=1 and args[1]=BFd1 then
print(`BFd1(A,d,K): inputs a list A of integers of length k say, and a positive integer d`):
print(`finds all solutions of`):
print(`add(A[i]*X[i]^d,i=1..k)=0 with X[i], i from 1 to k-1 from 1 to K. Try:`):
print(`BFd1([1,1,1,1],3,5);`):

elif nargs=1 and args[1]=BFd then
print(`BFd(A,d,K): inputs a list A of integers of length k say, and a positive integer d`):
print(`finds all good solutions of`):
print(`add(A[i]*X[i]^d,i=1..k)=0 with X[i], i from 1 to k-1 from 1 to K. Try:`):
print(`BFd([1,1,1,1],3,5);`):


elif nargs=1 and args[1]=BFdG then
print(`BFdG(A,d,K): inputs a list A of integers of length k say, and a positive integer d`):
print(`finds all solutions of`):
print(`add(A[i]*X[i]^d,i=1..k)=0 with X[i], i from 1 to k-1 from -K to K (except 0). Try:`):
print(`BFdG([1,1,1,1],3,5);`):

elif nargs=1 and args[1]=Box1 then
print(`Box1(L): All the vectors [x[1], .., x[nops(L)]] such that  x[i] belongs to L[i]. Try: `):
print(`Box1([{-1,1},{-2,2}]);`):


elif nargs=1 and args[1]=Cands then
print(`Cands(A,m,n,K,r): inputs a list of integers A of length 4 and symbols m,n,and a pos. integer K, and a pos. integer r.`):
print(`Outputs a set of quadruple of quadratic polynomials of size <=r `):
print(`such that for each of its members, add(A[i]*L[i]^3,i=1..4)=0. It will return the empty set if none found. K is the trial parameter. Try:`):
print(`Cands([1,1,1,1],m,n,5,1);`):

elif nargs=1 and args[1]=CleanUp then
print(`CleanUp(S,var): Given a list of tuples of polynomials in the list of variables var, kicks out redundanies, Try:`):
print(`CleanUp([m,-m],[m]); `):

elif nargs=1 and args[1]=EriJ then
print(`EriJ(E,L1,L2): uses Eri Jabotinsky's method to get yet-anothe solution of`):
print(`E[1]*X1^3+E[2]*X2^3+E[3]*X3^3+E[4]*X4^3=0 from two given solutions.`):
print(`Try:`):
print(`EriJ([1,1,1,1],[3,4,5,-6],[m,-m,n,-n]);`):


elif nargs=1 and args[1]=EriJs then
print(`EriJs(E,L1,L2,m,n): uses Eri Jabotinsky's method to get yet-anothe sultion of`):
print(` E[1]*X1^3+E[2]*X2^3+E[3]*X3^3+E[4]*X4^3=0 from two given solutions. `):
print(` here L1 and L2 are two known solutions possibly in terms of symbols m and n `):
print(` Try: `):
print(` EriJs([1,1,1,1],[3,4,5,-6],[m,-m,n,-n],m,n); `):

elif nargs=1 and args[1]=ExtractCoeffs then
print(`ExtractCoeffs(P,v): inputs a polynomial P in the list of variables v, outputs the set of coefficients of all monomials. Try:`):
print(`ExtractCoeffs(5*m1*m2+3*m1^2,[m1,m2]);`):

elif nargs=1 and args[1]=ExtractCoeffsOld then
print(`ExtractCoeffsOld(P,m,k): inputs a polynomial P in m[1], ..., m[k] outputs the set of coefficients of all monomials. Try:`):
print(`ExtractCoeffsOld(5*m[1]*m[2]+3*m[1]^2,m,2);`):


elif nargs=1 and args[1]=Find43 then
print(`Find43(A,L,m,n,M): Inputs a list of intgeres A of  length 4 and two numerical solution to`):
print(`add(A[i]*X[i]^3,i=1..4)=0 and a pos. integer M, finds all the good Morph43`):
print(`Morph43(A,L,M,m,n,M) (q.d.) where M is  a permutation of -L. Try:`):
print(`Find43([1,1,1,1],[3,4,5,-6],m,n,5);`):

elif nargs=1 and args[1]=FindDE then
print(`FindDE(L,m,X): inputs a list L of polynomials in m[1], ..., m[k-1] (where k=nops(L),  outputs the minimal polynomial`):
print(`in X[1], ... X[k] such that`):
print(`P(X[1],..., X[k]]) such that P(L[1], ..., L[k]) is identically zero. Try:`):
print(`FindDE([m[1]^2+1,m[1]^3+2],m,X);`):
print(`FindDE([m[1]^2-m[2]^2,2*m[1]*m[2],m[1]^2+m[2]^2],m,X);`):
print(`FindDE([m[1]^2-m[2]^2,2*m[1]*m[2],m[1]^2+m[2]^2],m,X);`):


elif nargs=1 and args[1]=FindEqsP then
print(`FindEqsP(P,x,L,v,Pars): Inputs`):
print(` (i)a polynomial P in the list of variables x `):
print(` (ii) a list L of the same length as x  of the variables v the coefficients may involve parameters. `):
print(` (iii) the set of parameters. `):
print(`Outputs: `):
print(` The set of coefficients in the variables in v of all monomials. Try: `):
print(`FindEqsP(X^2+Y^2-Z^2,[X,Y,Z],[a20*m^2+a11*m*n+a02*n^2,b20*m^2+b11*m*n+b02*n^2,c20*m^2+c11*m*n+c02*n^2],[m,n],{a20,a11,a02,b20,b11,b02,c20,c11,c02});`):


elif nargs=1 and args[1]=FindEqsPold then
print(`FindEqsPold(P,X,r,m,k,d,c): Inputs: `):
print(`(i)a polynomial P in X[1], .., X[r], (where X is a symbol)`):
print(` (ii) a symbol m `):
print(` (iii) pos. integers k and d `):
print(` (iv) a symbol c. Outputs `):
print(` (i) a list of length r with tenative expressions of generic polynomials with indexed variables given by c[i]'s `):
print(` let's call it L `):
print(` (ii) The set of equations in the c[i]'s such that P(L[1], ..., L[r]) is identically zero `):
print(` (iii) the set of c[i]'2. Try: `):
print( `FindEqsPold(X[1]^2+X[2]^2-X[3]^2,X,3,m,2,2,c); `):
print(` FindEqsPold(X[1]^3+X[2]^3+X[3]^3+X[4]^3,X,4,m,2,2,c); `):



elif nargs=1 and args[1]=FindGoodR then
print(`FindGoodR(m,n,K): finds all the solutions to X1^3+X2^3+X3^3+X4^3 in terms of (m,n) that`):
print(`yields one solution in positive integers in m,n<=K to  equalling 1 or -1 for one of its componets.`):
print(`Try:`):
print(`FindGoodR(m,n,100);`):

elif nargs=1 and args[1]=GenPol1 then
print(`GenPol1(m,k,d,c,st): Inputs a symbol m, a pos. integer k, a pos. integer d, and a symbol c, outputs`):
print(`A generic homog. polynomial in m[1], ..,m[k] with indexed coefficients c[i] starting at i=st. `):
print(`It also outputs the set of coefficients. Try:`):
print(`GenPol1(m,2,2,c,1);`):


elif nargs=1 and args[1]=IsTrivial then
print(`IsTrivial(L): decides whether L is trivial. Try:`):
print(` IsTrivial([2*n^2,3*n^2,5*n^2]); `):

elif nargs=1 and args[1]=Katan then
print(`Katan(P,m,n,K): inputs a polynomial P in m and n finds the minimum value (except 0),  in absolute value`):
print(`in [-K,K]^2 and where it is and whether it is posive or negative. Try:`):
print(`Katan(m^2-2*n^2,m,n,10);`):


elif nargs=1 and args[1]=MakeRama then
print(`MakeRama(L,m,n,t,K): inputs a list L of homog. quadratic polynomials in m and n, a variable t and a large pos. integer K`):
print(`outputs what MakeRama1(L,m,n,t,K) (q.v.)  let's call it [i,c,Pair]`):
print(`together with a tuple of the same length as L where`):
print(`the i-th entry is replaced by c and the other ones, way L[j], is the generating functions of the expressions for L[j]`):
print(`if m and n are replaced by Pair[1],Pair[2]. Try:`):
print(`MakeRama([m^2-7*m*n-9*n^2, 2*m^2+4*m*n+12*n^2, -2*m^2-10*n^2, -m^2-9*m*n+n^2],m,n,t,100);`):
print(`MakeRama([3*a^2+5*a*b-5*b^2, 4*a^2-4*a*b+6*b^2,5*a^2-5*a*b-3*b^2,-6*a^2+4*a*b-4*b^2],a,b,t,100);`):
print(`To reproduce the original Ramanujan result as explained in the article "Another look at an Amazing identity of Ramanujan"`):
print(`by  Jung Hun Han and Michael D. Hirschhorn, Mathematics Magazine Vol. 79, Nov. 4 (Oct. 2006), pp. 302-304; Try: `):
print(`MakeRama([u^2+7*u*v-9*v^2,2*u^2-4*u*v+12*v^2,-2*u^2-10*v^2,-u^2+9*u*v+v^2],u,v,t,100);`):

elif nargs=1 and args[1]=MakeRama1 then
print(`MakeRama1(L,m,n,t,K): inputs a list L of homog. quadratic polynomials in m and n`):
print(`finds the best entry i if it exists together with the constant. The output is [i,c,PairOfRationalFunctions of L[i]=c. Try:`):
print(`MakeRama1([m^2-7*m*n-9*n^2, 2*m^2+4*m*n+12*n^2, -2*m^2-10*n^2, -m^2-9*m*n+n^2],m,n,t,100);`):


elif nargs=1 and args[1]=MakeRamaSol then
print(`MakeRamaSol(L,m,n,t,K): Same as MakeRama(L,m,n,t,K) (q.v.) but using SolDE2 rather than SolQuad `):
print(`To reproduce the original Ramanujan result as explained in the article "Another look at an Amazing identity of Ramanujan"`):
print(`by  Jung Hun Han and Michael D. Hirschhorn, Mathematics Magazine Vol. 79, Nov. 4 (Oct. 2006), pp. 302-304; Try: `):
print(`MakeRamaSol([u^2+7*u*v-9*v^2,2*u^2-4*u*v+12*v^2,-2*u^2-10*v^2,-u^2+9*u*v+v^2],u,v,t,100);`):

elif nargs=1 and args[1]=Morph43 then
print(`Morph43(A,L1,L2,m,n,M): Given two specific solutions of add(A[i]*X[i]^3,i=1..4)=0, tries to find`):
print(`a doubly-infinite set of solutions where X[i] are quadratic polynomials in m and n with middle coefficients from -M to M. Try:`):
print(`Morph43([1,1,1,1],[3,4,5,-6],[-5,6,-3,-4],m,n,5); `):
print(` Morph43([1,1,1,1],[1,2,-2,-1],[-9,12,-10,1],m,n,10);`):


elif nargs=1 and args[1]=Paper then
print(`Paper(C,t,X,Y,Z,K1,K2): Inputs a positive integer B, variables t, X,Y,Z, and positive integers K1 and K2 (guessing parameters)`):
print(`outputs an article about Diophantine equations of the form`):
print(`A*X^3+B*Y^3+C*Z^3=Constant, for some constant for A<B<=C where A and B are relatively prime.`):
print(` Try: `):
print(` Paper(3,t,X,Y,Z,10,100): `):


elif nargs=1 and args[1]=PaperV then
print(`PaperV(C,m,n,t,X,Y,Z,K1,K2): Inputs a positive integer B, variables t, X,Y,Z, and positive integers K1 and K2 (guessing parameters)`):
print(`outputs a MOTIVATED article about Diophantine equations of the form`):
print(`A*X^3+B*Y^3+C*Z^3=Constant, for some constant for A<B<=C where A and B are relatively prime.`):
print(` Try: `):
print(` PaperV(3,m,n,t,X,Y,Z,10,100): `):

elif nargs=1 and args[1]=ParAB3 then
print(`ParAB3(A,B,m,n,K): inputs positive integers A and B and outputs a parametric equation, in terms of quadratic polynomials in m and n`):
print(`to the Diphantine equation `):
print(`A*X^3+A*Y^3+B*Z^3+B*W^3=0 using Eri Jabotinsky's trick by first guessing a particular solution of size <=K. If`):
print(`none found it returns FAIL. It returns a list of length 2. The first is the numerical solution found empircally usinf BFd, and the second`):
print(`the parametric equation using the trick. Try:`):
print(`ParAB3(2,3,m,n,8);`):

elif nargs=1 and args[1]=RandPolH then
print(`RandPolH(m,k,d,K): A random polynomial in m[1], ..., m[k], of homeg. degree k with coefficients from -K to K. Try:`):
print(`RandPolH(m,2,2,3);`):


elif nargs=1 and args[1]=Reject then
print(`Reject(P,var,N): Given a  polynomial with integer coefficients in the set of variables var and a positive integer N`):
print(`decides whether there exists a prime power less than N such that if you solve it mod this power there is no solution of`):
print(`P=0 and hence there is no solution. try:`):
print(`Reject(19*X^2-Y^2-2,[X,Y],30);`):

elif nargs=1 and args[1]=Rootk then
print(`Rootk(n,k): the kth-root of n or FAIL. Try:`):
print(`Rootk(-8,3);`):

elif nargs=1 and args[1]=SolBF then
print(`SolBF(L,Var,M): inputs a list L of equations in the variables Var, and  apos. integer M, finds all solutions in integers from -M to M. Try`):
print(`SolBF([x^2-1,x+y-3,x+y+z-10],{x,y,z},5);`):

elif nargs=1 and args[1]=SolBF1 then
print(`SolBF1(P,Var,M): Give a polynomial P in the list of variables Var finds all solutions of P(Var)=0`):
print(`with values from -M to M by brute force. Try:`):
print(`SolBF1(x^2-2*y^2-1,[x,y],20);`):


elif nargs=1 and args[1]=SolBF1a then
print(`SolBF1a(P,Var,M): Give a polynomial P in the list of variables Var finds all solutions of P(Var)=0`);
print(`with values from -M to M by brute force. Try:`):
print(`SolBF1a(x^2-2*y^2-1,[x,y],3);`):


elif nargs=1 and args[1]=SolDEP then
print(`SolDEP(P,x,L,v,Pars,M): Inputs `):
print(` (i)a polynomial P in the list of variables x `):
print(` (ii) a list L of the same length as x  of the variables v the coefficients may involve parameters. `):
print(` (iii) the set of parameters `):
print(` (iv) a positive integer M `):
print(`Outputs:`):
print(`The set of specializations of L with integers between -M and M such that P(L[1], ..., L[nops(L)])=0. Try:`):
print(`SolDEP(X^2+Y^2-Z^2,[X,Y,Z],[m^2+a11*m*n+a02*n^2,b20*m^2+b11*m*n+b02*n^2,c20*m^2+c11*m*n+c02*n^2],[m,n],{a11,a02,b20,b11,b02,c20,c11,c02},2);`):

print(`L:=[a02*n^2+a11*m*n+3*m^2, b02*n^2+b11*m*n+4*m^2, c02*n^2+c11*m*n+5*m^2, d02*n^2+d11*m*n-6*m^2]:`):
print(`SolDEP(X^3+Y^3+Z^3+W^3,[X,Y,Z,W],L,[m,n],{a02,a11,b02,b11,c02,c11,d02,d11},6);`):
print(`L:=[a02*n^2+a11*m*n+a20*m^2, b02*n^2+b11*m*n+b20*m^2, c02*n^2+c11*m*n+c20*m^2, d02*n^2+d11*m*n+d20*m^2]:`):
print(` var:={a02,a11,a20,b02,b11,b20,c02,c11,c20,d02,d11,d20}: `):
print(`L1:=subs({a02=3,b02=4,c02=5,d02=-6},L): var1:=var minus {a02,b02,c02,d02}:`):
print(`SolDEP(X^3+Y^3+Z^3+W^3,[X,Y,Z,W],L1,[m,n],var1,6);`):


elif nargs=1 and args[1]=SolDEPold then
print(`SolDEPold(P,X,r,m,k,d,M): Inputs`):
print(` (i)a polynomial P in X[1], .., X[r], (where X is a symbol) `):
print(` (ii) a symbol m `):
print(` (iii) pos. integers k and d. Outputs a set of parametric solutions to `):
print(` the diophantine equation P(X[1], ..., X[r])=0 where each solution is a list `):
print(` of length r whose entries are  homog. polynomials in m[1],...,m[k] of total degree d, and the `):
print(`coefficients are between -M and M, excluding trivial solutions. Try: `):
print(`SolDEPold(X[1]^2+X[2]^2-X[3]^2,X,3,m,2,2,2);`):
print(`SolDEPold(X[1]^3+X[2]^3+X[3]^3+X[4]^3,X,4,m,2,2,6);`):

elif nargs=1 and args[1]=SortLE then
print(`SortLE(L,Vars): inputs a list of polynomials L in the set of variables Vars, outputs the sorted`):
print(`version in increasing number of variables. Try:`):
print(`SortLE([x+y+z,z,y+z,x+z],{x,y,z});`):

elif nargs=1 and args[1]=TheoremAB then
print(`TheoremAB(A,B,m,n,t,X,Y,Z,K1,K2): Inputs positive integers A and B, letters m,n, t and X and guessing parameters K1 and K2`):
print(`(fairly large positive integers, start K1 around 10 and K2 around 100), outputs`):
print(`(i) A list of length 4, let's call it L, consisting of parametric equation for the solution of A*X[1]^3+A*X[2]^3+B*X[3]^3+B*X[4]^3=0 in terms of`):
print(` quadratic polynomials in m, n`):
print(`(ii) A list of length 4 of the form [Place,constant,f1,f2] where 1<=Place<=4, constant is a rational number and f1 and f2 are rational functions of `):
print(` the variable t such that if a(i) is the coefficient of t^i  in the power-series expansion of f1`):
print(`and b(i) is the coefficient of t^i  in the power-series expansion of f2 then L[Place](a(i),b(i))=consant for all i. In other words `):
print(`it is the solution of the Diphantine equation  L[Place](m,n)=constant in the C-finite ansatz`):
print(`(iii) A list of length 4 there of whose entries are rational functions in t and one of them is a constant such that`):
print(`they form solutions of A*(X[1]^3+X[2]^3)+B*(X[3]^3+X[4]^3)=0`):
print(`(iv) A Diophantine equation of the from A*X^3+B*Y^3+B*Z^3=CONSTANT `):
print(` (v): A triple of rational functions in t whose MacLaurin coefficients form solutions of that Diophantine equation. Try: `):
print(`TheoremAB(1,2,m,n,t,X,Y,Z,10,100);`):


elif nargs=1 and args[1]=TheoremABs then
print(`TheoremABs(A,B,t,X,Y,Z,K1,K2): States a theorem about infinitely many  solutions  to the`):
print(`diophantine equation`):
print(` A*X^3+B*Y^3+B*Z^3=Constant, for some constant, Constant. `):
print(`Try: `):
print(` TheoremABs(1,2,t,X,Y,Z,10,100): `):

elif nargs=1 and args[1]=TheoremABsV then
print(`TheoremABsV(A,B,m,n,t,X,Y,Z,K1,K2): States and EXPLAINS a  theorem about infinitely many  solutions  to the`):
print(`diophantine equation`):
print(`A*X^3+B*Y^3+B*Z^3=Constant, for some constant, Constant.`):
print(`Try:`):
print(`TheoremABsV(1,2,m,n,t,X,Y,Z,10,100):`):

else
 print(`There is no such thing as`, args):

fi:


end:


##start stuff from Pell.txt
ezraP1:=proc()
if args=NULL then

print(`The supporing procedures are:`):
print(` CheckDE,  DE,DEm, GFkbN, GuessRec,  PellC,  Seqk , UP `):

else
ezra(args):

fi:


end:

ezraP:=proc()
if args=NULL then

print(` Pell.txt: A Maple package for  `):

print(`The MAIN procedures are:`):
print(` AllPellGF, DErat, PellFt, PellGF, PellLike, PellR, PellRg , SolQuad `):


elif nargs=1 and args[1]=AllPellGF then
print(`AllPellGF(N,x,K): All the solutions to Pell equations up to n=N^2-1,using PellR (i.e. continuted fractions). Try:`):
print(`AllPellGF(10,x,300);`):

elif nargs=1 and args[1]=CheckDE then
print(`AllPellGF(N,x,K): All the solutions to Pell equations (using PellGF(n,x,K)), up to n=N^2-1. try:`):
print(`AllPellGF(10,x,1000);`):

elif nargs=1 and args[1]=CheckDE then
print(`CheckDE(k,a0,a1,b0,b1,N): checks DE(k,a0,a1,b0,b1,X,Y) up to N terms. Try:`):
print(`CheckDE(5,a0,a1,b0,b1,10);`):


elif nargs=1 and args[1]=DE then
print(`DE(k,a0,a1,b0,b1,X,Y,t): Finds a polynomial P(X,Y) such that if`):
print(`X=a(n) Y=b(n) if a(n) b(n) are solutions of the recurrence`):
print(`x(n)=k*x(n-1)-x(n-2) with a(0)=a0, a(1)=a(1), b(0)=b0, b(1)=b1.`):
print(`Then we have P(a(n),b(n))=0 `):
print(`It also returns the generating functions for a(n) and b(n) in terms of t`):
print(`DE(6,-1,1,1,1,X,Y,t);`):
print(` DE(6,1,3,0,2,X,Y,t);`):
print(`DE(3,1,2,3,4,X,Y,t); `):



elif nargs=1 and args[1]=DEm then
print(`DEm(k,a0,a1,b0,b1,X,Y,t): Finds a polynomial P in the variables X and Y such that`):
print(`if a(n) b(n) are solutions of the recurrence`):
print(`x(n)=k*x(n-1)+x(n-2) with a(0)=a0, a(1)=a(1), b(0)=b0, b(1)=b1.  `):
print(` then P(a(n),b(n))=(-1)^n `):
print(`It also returns the generating functions for a(n) and b(n). Try:`):
print(`DEm(9,0,1,1,9,X,Y,t);`):

elif nargs=1 and args[1]=DErat then
print(`DErat(Zug,t, X,Y): inputs a pair of rational functions Zug and tries to find`):
print(` a quadratic polynomil P in the variables X and Y such that`):
print(`if a(n) b(n) are the Taylor coefficients (about t=0) or Zug[1],Zug[2] `):
print(`then P(a(n),b(n))=0`):
print(`DErat([1/(1-3*t+t^2),t/(1-3*t+t^2)],t, X,Y);`):

elif nargs=1 and args[1]=GFkbN then
print(`GFkbN(k,b,N,t): Given integers or symbols k,b,N, outputs the generating function of`):
print(`the sequence a(n)^2-N*b(n)^2 if Sum(a(n)*x^n,n=0..infinity)= (1-k*x)/(1-2*k*x+x^2)`):
print(`Sum(b(n)*x^n,n=0..infinity)= (b*x)/(1-2*k*x+x^2). Try:`):
print(`GFkbN(3,2,2,t);`):

elif nargs=1 and args[1]=GuessRec then
print(`GuessRec(L): inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant cofficients `):
print(`satisfying it. It returns the initial  values and the operator`):
print(` as a list. For example try:`):
print(`GuessRec([1,1,1,1,1,1]);`):


elif nargs=1 and args[1]=PellC then
print(`PellC(N,M): inputs a non-perfect-square pos. integer N, and a large pos. integer M, and outputs the data (k,a0,a1,b0,b1) such that`):
print(`DE(k,a0,a1,b0,b1,X,Y,t)[1] is Pell's equation X^2-N*Y^2=1.`):
print(`It uses the contunued fraction of sqrt(N). Try: `):
print(`PellC(19,100);`):

elif nargs=1 and args[1]=PellFt then
print(`PellFt(N1,N2,M,t): inputs a non-perfect-square rel. prime pos. integer N1, N2,and a large pos. integer M, and outputs the list such that`):
print(`of pairs of rational functions whose coefficients satisfy N1*X^2-N2*Y^2=c for small constant. Try:`):
print(`PellFt(1,19,100,t);`):


elif nargs=1 and args[1]=PellCg then
print(`PellCg(N,M): inputs a non-perfect-square pos. integer N, and a large pos. integer M, and outputs the `):
print(`list of data [j,(k,a0,a1,b0,b1)] such that`):
print(`DE(k,a0,a1,b0,b1,X,Y,t)[1] is the  equation X^2-N*Y^2=j, for some CONST.`):
print(`It uses the contunued fraction of sqrt(N) `):
print(`PellCg(19,100);`):


elif nargs=1 and args[1]=PellGF then
print(`PellGF(n,x,K): solves Pell's equation X^2-N*Y^2=1, tries all y from 1 to K until 1+N*y^2 is a perfect square. `):
print(`Ooutputs the generating functions for a(n),b(n) such that a(n)^2-N*b(n)^2=1`):
print(`if none found, returns FAIL and K should be made bigger. `):
print(`It uses the pre-set generating functions (NOT using contiuned fractions explicitly)`):
print(`(1-k*x)/(1-2*k*x+x^2), b*x/(1-2*k*x+x^2)]):`):
print(`PellGF(2,x,100)`):

elif nargs=1 and args[1]=PellLike then
print(`PellLike(a0,a1,b0,t,X,Y): Gives the most general Pell-like eqation satisfied by`):
print(`satisfied by`):
print(`The coefficients of [(a0+a1*t)/(1-k*t+t^2),(b0+b1*t)/(1-k*t+t^2)] for the appropriate k`):
print(`Together with the equation of the from A1*X^2+A2*Y^2=A3`):
print(`Try: `):
print(`PellLike(a0,a1,b0,b1,t,X,Y); `):

elif nargs=1 and args[1]=PellR then
print(`PellR(N,M,t): inputs a non-perfect-square pos. integer N, and a large pos. integer M,  and a variable name t.`):
print(`Outputs the pair of rational functions A(t),B(t) such that their respective coefficients satisfy Pell's equation X^2-N*Y^2=1. `):
print(`It uses the contunued fraction of sqrt(N). Try: `):
print(` PellR(19,100,t); `):

elif nargs=1 and args[1]=PellRg then
print(`PellRg(N,M,t): inputs a non-perfect-square pos. integer N, and a large pos. integer M,  and a variable name t.`):
print(`Outputs the `):
print(`list of of triples [A(t),B(t),m] such that their respective coefficients satisfy Pell's equation X^2-N*Y^2=m. Try`):
print(`It uses the contunued fraction of sqrt(N). Try: `):
print(`PellRg(19,100,t);`):

elif nargs=1 and args[1]=Seqk then
print(`Seqk(k,a0,a1,N): The first N terms of the C-finite sequence satisying a(n)=k*a(n-1)-a(n-2)`):
print(` with a(0)=a0, a(1)=a0, Try: `):
print(` Seqk(5,1,2,20); `):

elif nargs=1 and args[1]=SolQuad then
print(`SolQuad(Q,m,n,t,K): Given a quadratic homog. polynomial Q in m and n finds solutions,  in terms of generating functions in t`):
print(`of Q(m,n)=D for select D. K is a guessing parameter. Try:`):
print(`SolQuad(m^2-2*n^2,m,n,t,100);`):


elif nargs=1 and args[1]=UP then
print(`UP(L): Given a sequence L, decides whether it is`):
print(`ultimately periodic, and`):
print(`returns the beginning and the periodic part`):
print(`For example, try: UP([3,1,4,2,4,2,4,2,4,2,4,2,4,2]);`):

else
 print(`There is no such thing as`, args):

fi:


end:





##############Start Detecting ultimate periodicity


#Shrink1(L): normalizes L=[L1,L2]
Shrink1:=proc(L) local L1,L2:
L1:=L[1]: L2:=L[2]:
if L1=[] then
  RETURN(L):
fi:

if L1[nops(L1)]=L2[nops(L2)] then
  RETURN([[op(1..nops(L1)-1,L1)],[L2[nops(L2)],op(1..nops(L2)-1,L2)]]):
else
  RETURN(L):
fi:
end:

Shrink:=proc(L) local L1:

L1:=L:

while L1<>Shrink1(L1) do
L1:=Shrink1(L1):
od:
L1:
end:

#UP1(L,p): is the sequence L ultimately periodic of period p?
#returns the beginning and the periodic part
#For example, try: UP1([3,1,4,2,4,2,4,2,4,2,4,2],2);
UP1:=proc(L,p) local i,m,gu,mu,i1,gu1,gu2:


m:=nops(L):
 if m-2*p+1<1 then
 print(`Too short`):
 RETURN(FAIL):
fi:

m:=nops(L):

if L[m]<>L[m-p] then 
 RETURN(FAIL):
fi:

if L[m-1]<>L[m-1-p] then 
 RETURN(FAIL):
fi:

if L[m-2]<>L[m-2-p] then 
 RETURN(FAIL):
fi:

if L[m-3]<>L[m-3-p] then 
 RETURN(FAIL):
fi:


if not [op(m-p+1..m,L)]=[op(m-2*p+1..m-p,L)] then
 RETURN(FAIL):
fi:


gu2:=[op(m-p+1..m,L)]:


for i from 1 while [op(i..i+p-1,L)]<>gu2 do od:

gu1:=[op(1..i-1,L)]:

gu:=[gu1,gu2]:
mu:=[op(gu1),
 seq(op(gu2),i1=1..trunc((m-nops(gu1)) /p) )
    ]:

     if mu<>[op(1..nops(mu),L)] then

       FAIL:
      else
       Shrink(gu):
      fi:

end:





#UP(L): Given a sequence L, decides whether it is
#ultimately periodic, and
#returns the beginning and the periodic part
#For example, try: UP([3,1,4,2,4,2,4,2,4,2,4,2]);
UP:=proc(L) local p,m,gu:

m:=nops(L):

for p from 1 to trunc((m-10)/2) do
 gu:=UP1(L,p):

  if gu<>FAIL then
    RETURN(gu):
  fi:
od:


FAIL:

end:
####end detecting ultimate periodicity


###start from Cfinite.txt

#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:

eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:




#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

for d from 1 to trunc(nops(L)/2)-2 do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:


#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:


#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
#   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):
L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
#print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(f):
fi:
end:


#RtoC(R,t): Given a rational function R(t) finds the
#C-finite description of its sequence of coefficients
#For example, try:
#RtoC(1/(1-t-t^2),t);
RtoC:=proc(R,t) local S1,S2,D1,R1,d,f1,i,a:
R1:=normal(R):
D1:=denom(R1):
d:=degree(D1,t):
if degree(numer(R1),t)>=d then
print(`The degree of the top must be less than the degree of the bottom`):
 RETURN(FAIL):
fi:

a:=coeff(D1,t,0):
S2:=[seq(-coeff(D1,t,i)/a,i=1..degree(D1,t))]:
f1:=taylor(R,t=0,d+1):
S1:=[seq(coeff(f1,t,i),i=0..d-1)]:
[S1,S2]:
end:

###end from Cfinite.txt

#Seqk(k,a0,a1,N): The first N terms of the C-finite sequence satisying a(n)=k*a(n-1)-a(n-2)
#with a(0)=a0, a(1)=a0, Try:
#Seqk(5,1,2,20);
Seqk:=proc(k,a0,a1,N) local gu,x,i:
gu:=[a0,a1]:
for i from 3 to N do
gu:=[op(gu),k*gu[-1]-gu[-2]]:
od:
gu:
end:



#DEold(k,a0,a1,b0,b1,t): Finds a polynomial equation satisfied by
#X=a(n) Y=b(n) if a(n) b(n) are solutions of the recurrence
#x(n)=k*x(n-1)-x(n-2) with a(0)=a0, a(1)=a(1), b(0)=b0, b(1)=b1. 
#It also returns the generating functions for a(n) and b(n)
#Try:
#DEold(3,1,2,3,4,X,Y,t) 
DEold:=proc(k,a0,a1,b0,b1,X,Y,t) local x,A1,A2,B1,B2,gu,Z,eq,P,A,B:
gu:=solve({A1+A2=a0,A1*x+A2/x=a1},{A1,A2}):
A1:=subs(gu,A1):A2:=subs(gu,A2):
gu:=solve({B1+B2=b0,B1*x+B2/x=b1},{B1,B2}):
B1:=subs(gu,B1):B2:=subs(gu,B2):

eq:={numer(X-A1*Z-A2/Z), numer(Y-B1*Z-B2/Z),numer(x+1/x-k)}:

P:=Groebner[Basis](eq,plex(x,Z,X,Y))[1]:

A:=CtoR([[a0,a1],[k,-1]],t):
B:= CtoR([[b0,b1],[k,-1]],t):

[P,[A,B]]:
end:

#CheckDE(k,a0,a1,b0,b1,N): checks DE(k,a0,a1,b0,b1,X,Y) up to N terms. Try:
#CheckDE(5,a0,a1,b0,b1,10);
CheckDE:=proc(k,a0,a1,b0,b1,N) local gu,mu1,mu2,X,Y,i,t:
gu:=DE(k,a0,a1,b0,b1,X,Y,t)[1]:
mu1:=Seqk(k,a0,a1,N):
mu2:=Seqk(k,b0,b1,N):
evalb({seq(expand(subs({X=mu1[i],Y=mu2[i]},gu)),i=1..N)}={0}):
end:



#PellC(N,M): inputs a non-perfect-square pos. integer N, and a large pos. integer M, and outputs the data (k,a0,a1,b0,b1) such that
#DE(k,a0,a1,b0,b1,X,Y,t)[1] is Pell's equation X^2-N*Y^2=CONST, for some CONST. Try:
#PellC(19,100);
PellC:=proc(N,M) local lu,mu,mu1,i,p,L1,L2,C1,C2,a0,a1,b0,b1,k,X,Y,t:

if not (type(N,integer) and N>1 and not type (sqrt(N),integer) ) then
# print(`Bad input`):
 RETURN(FAIL):
fi:

convert(sqrt(N),confrac,M,lu):
mu:=[seq(numer(lu[i])^2-N*denom(lu[i])^2,i=1..nops(lu))]:
mu1:=UP(mu):

if mu1[1]<>[] then
# print(`Not yet implemented`):
 RETURN(FAIL):
fi:

mu1:=mu1[2]:

if mu1[nops(mu1)]<>1 then
# print(`Not yet implemented`):
 RETURN(FAIL):
fi:

p:=nops(mu1):

L1:=[seq(numer(lu[i*p]),i=1..trunc(nops(lu)/p) )]:
C1:=GuessRec(L1):
if C1=FAIL then
# print(`Make `, M, `bigger `):
 RETURN(FAIL):
fi:


L2:=[seq(denom(lu[i*p]),i=1..trunc(nops(lu)/p) )]:
C2:=GuessRec(L2):
if C2=FAIL then
# print(`Make `, M, `bigger `):
 RETURN(FAIL):
fi:


if C1[2][2]<>-1 then
 RETURN(FAIL):
fi:

if C2[2][2]<>-1 then
 RETURN(FAIL):
fi:

k:=C1[2][1]:

if C2[2][1]<>k then
 RETURN(FAIL):
fi:


a0:=C1[1][1]:a1:=C1[1][2]:
b0:=C2[1][1]:b1:=C2[1][2]:


lu:=[k,k*a0-a1,a0,k*b0-b1,b0]:

if not type(normal(DE(op(lu),X,Y,t)[1]/(X^2-N*Y^2-1)),numeric) then
print(lu, `gave `):
print( DE(op(lu),X,Y,t)[1]):
RETURN(FAIL):
fi:

lu:

end:




#PellR(N,M,t): inputs a non-perfect-square pos. integer N, and a large pos. integer M,  and a variable name t.
#Outputs the pair of rational functions A(t),B(t) such that their respective coefficients satisfy Pell's equation X^2-N*Y^2=1. Try
#PellR(19,100,t);
PellR:=proc(N,M,t) local gu,A,B,A1,B1,i:
gu:=PellC(N,M):

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

A:=CtoR([[gu[2],gu[3]],[gu[1],-1]],t):
B:= CtoR([[gu[4],gu[5]],[gu[1],-1]],t):
A1:=taylor(A,t=0,101):
B1:=taylor(B,t=0,101):
if {seq(coeff(A1,t,i)^2-N*coeff(B1,t,i)^2,i=1..100)}<>{1} then
 RETURN(FAIL):
fi:

[A,B]:
end:



#GFkbN(k,b,N,t): Given integers or symbols k,b,N, outputs the generating function of
#the sequence a(n)^2-N*b(n)^2 if Sum(a(n)*x^n,n=0..infinity)= (1-k*x)/(1-2*k*x+x^2)
#Sum(b(n)*x^n,n=0..infinity)= (b*x)/(1-2*k*x+x^2). Try:
#GFkbN(3,2,2,t);
GFkbN:=proc(k,b,N,t) local gu1,gu2,gu,i,x:

gu1:=taylor((1-k*x)/(1-2*k*x+x^2),x=0,11):
gu2:=taylor(b*x/(1-2*k*x+x^2),x=0,11):

gu:=[seq(coeff(gu1,x,i)^2- N*coeff(gu2,x,i)^2,i=0..10)]:

CtoR(GuessRec(gu),t):
end:





#PellCg(N,M): inputs a non-perfect-square pos. integer N, and a large pos. integer M, and outputs the 
#list of data [j,(k,a0,a1,b0,b1)] such that
#DE(k,a0,a1,b0,b1,X,Y,t)[1] is the  equation X^2-N*Y^2=j, for some CONST. Try:
#PellCg(19,100);
PellCg:=proc(N,M) local lu,mu,mu1,i,p,L1,L2,C1,C2,a0,a1,b0,b1,k,X,Y,j,gu,vu,kak,t:


if not (type(N,integer) and N>1 and not type (sqrt(N),integer) ) then
# print(`Bad input`):
 RETURN(FAIL):
fi:

convert(sqrt(N),confrac,M,lu):
mu:=[seq(numer(lu[i])^2-N*denom(lu[i])^2,i=1..nops(lu))]:
mu1:=UP(mu):

if mu1[1]<>[] then
# print(`Not yet implemented`):
 RETURN(FAIL):
fi:

mu1:=mu1[2]:


p:=nops(mu1):

gu:=[]:


for j from 1 to p do
vu:=mu1[j]:

L1:=[seq(numer(lu[i*p+j-p]),i=1..trunc(nops(lu)/p)-1 )]:


C1:=GuessRec(L1):
if C1=FAIL then
# print(`Make `, M, `bigger `):
 RETURN(FAIL):
fi:


L2:=[seq(denom(lu[i*p+j-p]),i=1..trunc(nops(lu)/p)-1 )]:


C2:=GuessRec(L2):
if C2=FAIL then
# print(`Make `, M, `bigger `):
 RETURN(FAIL):
fi:


if C1[2][2]<>-1 then
 RETURN(FAIL):
fi:

if C2[2][2]<>-1 then
 RETURN(FAIL):
fi:

k:=C1[2][1]:

if C2[2][1]<>k then
 RETURN(FAIL):
fi:


a0:=C1[1][1]:a1:=C1[1][2]:
b0:=C2[1][1]:b1:=C2[1][2]:


kak:=[k,k*a0-a1,a0,k*b0-b1,b0]:

if not type(normal(DE(op(lu),X,Y,t)[1]/(X^2-N*Y^2-vu)),numeric) then
print( DE(op(kak),X,Y,t)[1]):
print(`did not work out`):
RETURN(FAIL):
fi:

gu:=[op(gu),[vu,kak]]:

od:

gu:

end:



#PellRg(N,M,t): inputs a non-perfect-square pos. integer N, and a large pos. integer M,  and a variable name t.
#Outputs the 
#list of of triples [A(t),B(t),m] such that their respective coefficients satisfy Pell's equation X^2-N*Y^2=m. Try
#PellRg(19,100,t);
PellRg:=proc(N,M,t) local gu,A,B,A1,B1,i,j,mu:
gu:=PellCg(N,M):

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

mu:=[]:

for j from 1 to nops(gu) do

A:=CtoR([[gu[j][2][2],gu[j][2][3]],[gu[j][2][1],-1]],t):
B:= CtoR([[gu[j][2][4],gu[j][2][5]],[gu[j][2][1],-1]],t):
A1:=taylor(A,t=0,101):
B1:=taylor(B,t=0,101):
if {seq(coeff(A1,t,i)^2-N*coeff(B1,t,i)^2,i=1..100)}<>{gu[j][1]} then
 RETURN(FAIL):
fi:

mu:=[op(mu), [A,B,gu[j][1]] ]:
od:

mu:

end:



#PellGF(n,x,K): solves Pell's equation X^2-N*Y^2=1, tries all b from 1 to K, outputs the generating functions for a(n),b(n)
#if none found, returns FAIL. It uses the pre-set generating functions (NOT using contiuned fractions explicitly). Try:
#PellGF(n,x,K)
PellGF:=proc(n,x,K) local k,b:

if type(sqrt(n),integer) then
 RETURN(FAIL):
fi:

   
for b from 1 to K while not type(sqrt(1+n*b^2),integer) do od:

if b=K+1 then
 RETURN(FAIL):
else
 k:=sqrt(1+n*b^2):
RETURN([(1-k*x)/(1-2*k*x+x^2), b*x/(1-2*k*x+x^2)]):
fi:

end:


#AllPellGFold(N,x,K): All the solutions to Pell equations up to n=N^2-1,using PellR (i.e. continuted fractions). Try:
#AllPellGFold(10,x,300);
AllPellGFold:=proc(N,x,K) local gu,i:
gu:={seq(i,i=1..N^2-1)} minus {seq(i^2,i=1..N)};
gu:=sort(gu):
[seq([gu[i],PellR(gu[i],K,x)],i=1..nops(gu))]:
end:


#AllPellGF(N,x,K): All the solutions to Pell equations up to n=N^2-1,using PellR (i.e. continuted fractions). Try:
#AllPellGF(10,x,300);
AllPellGF:=proc(N,x,K) local gu,i,mu,lu,vu,f1,f2,k,b:
gu:={seq(i,i=1..N^2-1)} minus {seq(i^2,i=1..N)};
gu:=sort(gu):

vu:=[]:

for i from 1 to nops(gu) do

lu:=PellFt(1,gu[i],K,x)[1]:
if lu[1]<>1 then
 RETURN(FAIL):
else

f1:=lu[2]:
f2:=lu[3]:

k:=-coeff(denom(f1),x,1)/2:

b:=-coeff(numer(f2),x):

if k^2-gu[i]*b^2<>1 then
 RETURN(FAIL):
else
f1:=(1-k*x)/(1-2*k*x+x^2):
f2:=b*x/(1-2*k*x+x^2):
vu:=[op(vu),[gu[i],[f1,f2] ]]:
fi:

fi:
od:
vu:

end:



#DE(k,a0,a1,b0,b1,X,Y,t): Finds a polynomial P in the variables X and Y such that
#if a(n) b(n) are solutions of the recurrence
#x(n)=k*x(n-1)-x(n-2) with a(0)=a0, a(1)=a(1), b(0)=b0, b(1)=b1. 
#then P(a(n),b(n))=0
#It also returns the generating functions for a(n) and b(n)
#Try:
#DE(3,1,2,3,4,X,Y,t) 
DE:=proc(k,a0,a1,b0,b1,X,Y,t) local gu1,gu2,P,c1,c2,c3,eq,var,i:

gu1:=SeqFromRec([[a0,a1],[k,-1]],15):

gu2:=SeqFromRec([[b0,b1],[k,-1]],15):

P:=c1*X^2+c2*X*Y+c3*Y^2+1:
var:={c1,c2,c3}:

eq:={seq(subs({X=gu1[i],Y=gu2[i]},P),i=1..6)}:

var:=solve(eq,var):

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


P:=numer(subs(var,P)):


if {seq(subs({X=gu1[i],Y=gu2[i]},P),i=1..15)}<>{0} then
 RETURN(FAIL):
fi:

[P,[CtoR([[a0,a1],[k,-1]],t),CtoR([[b0,b1],[k,-1]],t)]]:

end:





#DEm(k,a0,a1,b0,b1,X,Y,t): Finds a polynomial P in the variables X and Y such that
#if a(n) b(n) are solutions of the recurrence
#x(n)=k*x(n-1)+x(n-2) with a(0)=a0, a(1)=a(1), b(0)=b0, b(1)=b1. 
#then P(a(n),b(n))=(-1)^n
#It also returns the generating functions for a(n) and b(n)
#Try:
#DEm(9,0,1,1,9,X,Y,t);
DEm:=proc(k,a0,a1,b0,b1,X,Y,t) local gu1,gu2,P,c1,c2,c3,eq,var,i:

gu1:=SeqFromRec([[a0,a1],[k,1]],15):

gu2:=SeqFromRec([[b0,b1],[k,1]],15):

P:=c1*X^2+c2*X*Y+c3*Y^2:
var:={c1,c2,c3}:

eq:={seq(subs({X=gu1[i],Y=gu2[i]},P)-(-1)^(i-1),i=1..6)}:

var:=solve(eq,var):

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

P:=subs(var,P):


if {seq(subs({X=gu1[i],Y=gu2[i]},P)-(-1)^(i-1),i=1..15)}<>{0} then
 RETURN(FAIL):
fi:

[P,[CtoR([[a0,a1],[k,1]],t),CtoR([[b0,b1],[k,1]],t)]]:

end:






#DErat(Zug,t, X,Y): inputs a pair of rational functions Zug and tries to find
#Finds a quadratic polynomil P in the variables X and Y such that
#if a(n) b(n) are the Taylor coefficients (about t=0) or Zug[1],Zug[2] then
#then P(a(n),b(n))=0
#DErat([1/(1-3*t+t^2),t/(1-3*t+t^2)],t, X,Y);
DErat:=proc(Zug,t,X,Y) local gu1,gu2,P,c1,c2,c3,eq,var,i,gu:

gu1:=taylor(Zug[1],t=0,17):
gu1:=normal([seq(coeff(gu1,t,i),i=1..16)]):

gu2:=taylor(Zug[2],t=0,17):
gu2:=normal([seq(coeff(gu2,t,i),i=1..16)]):

P:=c1*X^2+c2*X*Y+c3*Y^2+1:
var:={c1,c2,c3}:

eq:={seq(subs({X=gu1[i],Y=gu2[i]},P),i=1..6)}:

var:=solve(eq,var):

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


P:=numer(subs(var,P)):


if normal({seq(subs({X=gu1[i],Y=gu2[i]},P),i=1..16)})<>{0} then
 RETURN(FAIL):
fi:

P:=factor(P):

gu:=P:

if type(P,`*`) then

for i from 1 to nops(P) do

if (normal(diff(op(i,P),X))=0 and normal(diff(op(i,P),Y))=0) then
 gu:=factor(normal(gu/op(i,P))):
fi:
od:
fi:

factor(gu):

end:


#PellLike(a0,a1,b0,t,X,Y): Gives the most general pair of rational functions and the resulting Pell-Like equation
#satisfied by its coefficients
#Try:
#PellLike(a0,a1,b0,b1,t,X,Y);
PellLike:=proc(a0,a1,b0,b1,t,X,Y) local k,Zug,gu,lu,k0:

Zug:=[(a0+a1*t)/(1-k*t+t^2),(b0+b1*t)/(1-k*t+t^2)]:

gu:=DErat(Zug,t,X,Y):

lu:=coeff(coeff(gu,X,1),Y,1):


k0:=solve(lu,k):


Zug:=normal(subs(k=k0,Zug)):

gu:=DErat(Zug,t,X,Y):

if normal(coeff(coeff(gu,X,1),Y,1))<>0 then
 RETURN(FAIL):
fi:

gu:=factor(coeff(gu,X,2))*X^2+factor(coeff(gu,Y,2))*Y^2+ factor(coeff(coeff(gu,X,0),Y,0)):
Zug:=[collect(Zug[1],t),collect(Zug[2],t)]:

[Zug,gu]:

end:


###start N1*X^2-N2*Y^2=1


#PellFt1(N1,N2,M,t): inputs a non-perfect-square rel. prime pos. integer N1, N2,and a large pos. integer M, and outputs the data (k,a0,a1,b0,b1) such that
#DE(k,a0,a1,b0,b1,X,Y,t)[1] is the equation N1*X^2-N2*Y^2=CONST, for some CONST. Try:
#PellFt1(1,19,100,t);
PellFt1:=proc(N1,N2,M,t) local lu,mu,mu1,i,p,L1,L2,C1,C2,a0,a1,b0,b1,k,X,Y,j,ku,f1,f2,la,T,vu:

if not (type(N1,integer) and N1>0 ) then
 RETURN(FAIL):
fi:

if not (type(N2,integer) and N2>1 and not type (sqrt(N2),integer) ) then
 RETURN(FAIL):
fi:

if igcd(N1,N2)<>1 then
 RETURN(FAIL):
fi:

convert(sqrt(N2/N1),confrac,M,lu):
mu:=[seq(N1*numer(lu[i])^2-N2*denom(lu[i])^2,i=1..nops(lu))]:
mu1:=UP(mu):

if mu1[1]<>[] then
 RETURN(FAIL):
fi:

mu1:=mu1[2]:


p:=nops(mu1):

ku:=[]:

for j from 1 to nops(mu1) do

L1:=[seq(numer(lu[i*p+j]),i=1..trunc(nops(lu)/p)-1 )]:
C1:=GuessRec(L1):
if C1=FAIL then
 RETURN(FAIL):
fi:


L2:=[seq(denom(lu[i*p+j]),i=1..trunc(nops(lu)/p)-1 )]:
C2:=GuessRec(L2):
if C2=FAIL then
 RETURN(FAIL):
fi:

f1:=CtoR(C1,t): f2:=CtoR(C2,t):

if {seq(N1*coeff(taylor(f1,t=0,41),t,i)^2-N2*coeff(taylor(f2,t=0,41),t,i)^2,i=0..40)}<>{mu1[j]} then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

la:=[mu1[j],f1,f2]:

ku:=[op(ku),la]:
od:

lu:=convert(sort([seq(abs(ku[i][1]),i=1..nops(ku))]),set):

lu:=[seq(op([lu[i],-lu[i]]),i=1..nops(lu))]:

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

for i from 1 to nops(ku) do
 T[ku[i][1]]:=T[ku[i][1]] union {ku[i]}:
od:

vu:=[]:

for i from 1 to nops(lu) do

if T[lu[i]]<>{} then
 vu:=[op(vu),T[lu[i]][1] ] :
fi:

od:

vu:

end:


#PellFt(N1,N2,M,t): inputs a non-perfect-square rel. prime pos. integer N1, N2,and a large pos. integer M, and outputs the data (k,a0,a1,b0,b1) such that
#DE(k,a0,a1,b0,b1,X,Y,t)[1] is the equation N1*X^2-N2*Y^2=CONST, for some CONST. Try:
#PellFt(1,19,100,t);
PellFt:=proc(N1,N2,M,t) local  M1,gu:

M1:=10:

while M1<=M do
 gu:=PellFt1(N1,N2,M1,t):
 if gu<>FAIL then
   RETURN(gu):
 else
 M1:=M1*2:
 fi:
od:

FAIL:
end:



#SolQuad(Q,m,n,t,K): Given a quadratic homog. polynomial Q in m and n finds solutions,  in terms of generating functions in t
#of Q(m,n)=D for some constant D. 
#It reuturns a pair of rational functions [f1,f2] in t togetger with the constant. Try:
#SolQuad(m^2-2*n^2,m,n,t,K);
SolQuad:=proc(Q,m,n,t,K) local A,B,C,D1,gu,f1,f2,f1t,f2t,lu,i:

A:=coeff(Q,m,2):
B:=coeff(coeff(Q,m,1),n,1):
C:=coeff(Q,n,2):

if not (type(A,integer) and type(B,integer) and type(C,integer)) then
 print(A,B,C, `should have been integers `):
 RETURN(FAIL):
fi:

if expand(Q-A*m^2-B*m*n-C*n^2)<>0 then
 print(Q, `is not a homog.quadratic form in`, m,n ):
 RETURN(FAIL):
fi:

D1:= B^2-4*A*C:

if D1<=0 then
# print(`The discriminat is not positive`):
 RETURN(FAIL):
fi:

gu:=PellR(D1,K,t):

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

f2:=gu[2]:

f1:=normal((gu[1]-B*gu[2])/(2*A)):


f1t:=taylor(f1,t=0,41):
f2t:=taylor(f2,t=0,41):

lu:={seq(subs({m=coeff(f1t,t,i), n=coeff(f2t,t,i)},Q),i=0..40)}:



if nops(lu)<>1 then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

[f1,f2,lu[1]]:

end:

#####End  stuff from Pell.txt

#######################Start  stuff from RatDio.txt
ezra1R:=proc()
if args=NULL then

print(`The SUPPORTING procedures for Rational functions are`):
print(`   DE,DEm, FindEqs, FindEqs1 , FindEqsM, FindVars,  GenPol, GenPolHst, SolBF, SolBF1, SortLE, SR1  `):

else
ezra(args):

fi:

end:


ezraR:=proc()
if args=NULL then

print(`The MAIN procedures for rational functions are: DEg, DEgM, DEgO, SolDE2, SolDE2m , SR `):
print(``):


elif nargs=1 and args[1]=DE then
print(`DE(k,a0,a1,b0,b1,X,Y,t): Finds a polynomial P in the variables X and Y such that`):
print(`if a(n) b(n) are solutions of the recurrence`):
print(`x(n)=k*x(n-1)-x(n-2) with a(0)=a0, a(1)=a(1), b(0)=b0, b(1)=b1. `):
print(`then P(a(n),b(n))=0`):
print(`It also returns the generating functions for a(n) and b(n)`):
print(`Try:`):
print(`DE(3,1,2,3,4,X,Y,t);`):


elif nargs=1 and args[1]=DEg then
print(`DEg(L,t,X,d): inputs a list L of rational functions in t of length k say in the varible t,`):
print(`a letter X and a pos. integer d,outputs a polynomial in X[1], ..., X[k] of degree<= d `):
print(`let's call it P(X[1], ..., X[k]), such`):
print(`that if you plug-in the respective Taylor coefficients of L[i] are solutions of the`):
print(`Diphantine equation P(X[1], ..., X[k])=0. Try:`):
print(`DEg([(1+t)/(1-3*t+t^2) , 4/(1-3*t+t^2) ],t,X,2);`):
print(`DEg([(c0+c1*t)/(1-k*t+t^2) , d0*(1-(c1*k+2*c0)/(c0*k+2*c1)*t)/(1-k*t+t^2) ],t,X,2);`):
print(`DEg( subs(c1=(1-k*c0)/2,[(c0+c1*t)/(1-k*t+t^2) , d0*(1-(c1*k+2*c0)/(c0*k+2*c1)*t)/(1-k*t+t^2) ]),t,X,2);`):
print(`DEg([1/2*(c0*k*t-2*c0-t)/(k*t-t^2-1), -1/2*d0*(c0*k^2*t-4*c0*t-k*t+2)/(k*t-t^2-1)],t,X,2);`):
#print(`DEg([(1+4273*t-791*t^2)/(1-t)/(1-6887*t+t^2) , 2*(1-1154*t+505*t^2)/(1-t)/(1-6887*t+t^2) , (2+482*t+812*t^2)/(1-t)/(1-6887*t+t^2) ],t,X,3);`):



elif nargs=1 and args[1]=DEgM then
print(`DEgM(L,t,X,d): inputs a list L of rational functions in t of length k say in the varible t,`):
print(`a letter X and a pos. integer d,outputs a polynomial in X[1], ..., X[k] of degree d plus a constant term `):
print(`let's call it P(X[1], ..., X[k]), such`):
print(`that if you plug-in the respective Taylor coefficients of L[i] are solutions of the`):
print(`Diphantine equation P(X[1], ..., X[k])=(-1)^k. Try:`):
print(`DEgM([1/(1-3*t-t^2) , t/(1-3*t-t^2) ],t,X,2);`):
print(`DEgM([1/(1-3*t-t^3) , t/(1-3*t-t^3), t^2/(1-3*t-t^3) ],t,X,3);`):

elif nargs=1 and args[1]=DEgO then
print(`DEgO(L,t,X,d): inputs a list L of rational functions in t of length k say in the varible t,`):
print(`a letter X and a pos. integer d,outputs a polynomial in X[1], ..., X[k] of degree d plus a constant term `):
print(`let's call it P(X[1], ..., X[k]), such`):
print(`that if you plug-in the respective Taylor coefficients of L[i] are solutions of the`):
print(`Diphantine equation P(X[1], ..., X[k])=1. Try:`):
print(`DEgO([(1+t)/(1-3*t+t^2) , 4/(1-3*t+t^2) ],t,X,2);`):
print(`DEgO([(c0+c1*t)/(1-k*t+t^2) , d0*(1-(c1*k+2*c0)/(c0*k+2*c1)*t)/(1-k*t+t^2) ],t,X,2);`):
print(`DEgO([(1+4273*t-791*t^2)/(1-t)/(1-6887*t+t^2) , 2*(1-1154*t+505*t^2)/(1-t)/(1-6887*t+t^2) , (2+482*t+812*t^2)/(1-t)/(1-6887*t+t^2) ],t,X,3);`):


elif nargs=1 and args[1]=DEm then
print(`DEm(k,a0,a1,b0,b1,X,Y,t): Finds a polynomial P in the variables X and Y such that`):
print(`if a(n) b(n) are solutions of the recurrence`):
print(`x(n)=k*x(n-1)+x(n-2) with a(0)=a0, a(1)=a(1), b(0)=b0, b(1)=b1. `):
print(`then P(a(n),b(n))=(-1)^n`):
print(`It also returns the generating functions for a(n) and b(n)`):
print(`Try: `):
print(`DEm(9,0,1,1,9,X,Y,t);`):

elif nargs=1 and args[1]=FindEqs then
print(`FindEqs(L,t,Vars,P,X,Y): inputs a pair of rational functions in t phrased in terms of parameters in Vars`):
print(`and a  polynomial P in the variables X and Y, finds a Groebner basis `):
print(`such that the  if a(n) and b(n) are the coefficients of t^n of L[1], L[2] respectively then`):
print(`P(a(n),b(n))=0. If there is no solution it returns FAIL. Try:`):
print(`FindEqs([(a0+a1*t)/(1-k*t+t^2),(b0+b1*t)/(1-k*t+t^2)],t,[a0,a1,b0,b1,k],X^2-2*Y^2-1,X,Y);`):


elif nargs=1 and args[1]=FindEqs1 then

print(`FindEqs1(L,t,Vars,P,X,Y,K): inputs a pair of rational functions in t phrased in terms of parameters in Vars`):
print(`and a  polynomial P in the variables X and Y, finds a Groebner basis for the set of  nops(Vars)+K equations`):
print(`such that the  if a(n) and b(n) are the coefficients of t^n of L[1], L[2] respectively then`):
print(`P(a(n),b(n))=0. If there is no solution it returns FAIL. Try:`):
print(`FindEqs1([(a0+a1*t)/(1-6*t+t^2),(b0+b1*t)/(1-6*t+t^2)],t,[a0,a1,b0,b1],X^2-2*Y^2-1,X,Y,4);`):



elif nargs=1 and args[1]=FindEqsM then
print(`FindEqsM(L,t,Vars,P,X,Y): inputs a pair of rational functions in t phrased in terms of parameters in Vars`):
print(`and a  polynomial P in the variables X and Y, finds a Groebner basis for the set of nops(Vars)+K equations `):
print(`such that the  if a(n) and b(n) are the coefficients of t^n of L[1], L[2] respectively then`):
print(`P(a(n),b(n))=(-1)^n. If there is no solution it returns FAIL. Try:`):
print(`FindEqsM([(a0+a1*t)/(1-9*t-t^2),(b0+b1*t)/(1-9*t-t^2)],t,[a0,a1,b0,b1],X^2-9*X*Y-Y^2,X,Y);`):

elif nargs=1 and args[1]=FindVars then
print(`FindVars(P,Vars): inputs a polynomial in a subset of the set of variables Vars,`):
print(`finds the subest that actually participates in it. Try:`):
print(`FindVars(x^3+y*z,{x,y,z,w});`):




elif nargs=1 and args[1]=GenPol then
print(`GenPol(m,k,d,c): Inputs a symbol m, a pos. integer k, a pos. integer d, and a symbol c, outputs`):
print(`A generic polynomial in m[1], ..,m[k] of degree d with indexed coefficients c[i] and starting at i=st. `):
print(`It also outputs the set of coefficients. Try:`):
print(`GenPol(m,2,2,c);`):


elif nargs=1 and args[1]=GenPolHst then
print(`GenPolHst(m,k,d,c,st): Inputs a symbol m, a pos. integer k, a pos. integer d, and a symbol c, outputs`):
print(`A generic homog. polynomial in m[1], ..,m[k] with indexed coefficients c[i] starting at i=st. `):
print(`It also outputs the set of coefficients. Try:`):
print(`GenPolHst(m,2,2,c,1);`):

elif nargs=1 and args[1]=SolBF then
print(`SolBF(L,Var,M): inputs a list L of equations in the variables Var, and  apos. integer M, finds all solutions in integers from -M to M. Try`):
print(`SolBF([x^2-1,x+y-3,x+y+z-10],{x,y,z},5);`):

elif nargs=1 and args[1]=SolBF1 then
print(`SolBF1(P,Var,M): Give a polynomial P in the list of variables Var finds all solutions of P(Var)=0`):
print(`with values from -M to M by brute force. Try:`):
print(`SolBF1(x^2-2*y^2-1,[x,y],20);`):


elif nargs=1 and args[1]=SolDE2 then
print(`SolDE2(P,X,Y,t,M): Inputs a polynomial P of degree 2  w/o linear terms, with integer coefficients, in the variables X and Y`):
print(`and a variable t, and a positive integer M, finds, if it exists`):
print(`two rational functions of the form [(a0+a1*t)/(1-k*t+t^2),(b0+b1*t)/(1-k*t+t^2)] with`):
print(`k positive and larger than 2, and a0,a1,b0,b1 between -M and M, such that`):
print(`all the respective Taylor coefficients satisfy P(X,Y)=0. For example, to solve `):
print(`Pell's equation X^2-2*Y^2=1 type: `):
print(`SolDE2(X^2-2*Y^2-1,X,Y,t,6);`):


elif nargs=1 and args[1]=SolDE2m then
print(`SolDE2m(P,X,Y,t,M): Inputs a polynomial P , with integer coefficients, in the variables X and Y`):
print(`and a variable t, and a positive integer M, finds, if it exists`):
print(`two rational functions of the form [(a0+a1*t)/(1-k*t-t^2),(b0+b1*t)/(1-k*t-t^2)] with`):
print(`k positive and larger than 2, and a0,a1,b0,b1 between -M and M, such that`):
print(`all the respective Taylor coefficients a(n), b(n) satisfy P(a(n),b(n))=(-1)^n. For example, Try:`):
print(`SolDE2m(X^2-9*X*Y-Y^2,X,Y,t,11);`):


elif nargs=1 and args[1]=SolDE2e then
print(`SolDE2e(P,X,Y,t,M): Inputs a polynomial P , with integer coefficients, in the variables X and Y`):
print(`and a variable t, and a positive integer M, finds, if it exists`):
print(`two rational functions of the form [(a0+a1*t)/(1-k*t+t^2),(b0+b1*t)/(1-k*t+t^2)] with`):
print(`k positive and larger than 2, and a0,a1,b0,b1 between -M and M, such that`):
print(`all the respective EVEN Taylor coefficients satisfy P(X,Y)=0. For example, to solve `):
print(`Pell's equation X^2-2*Y^2=1 type: `):
print(`SolDE2e(X^2-2*Y^2-1,X,Y,t,6);`):

elif nargs=1 and args[1]=SortLE then
print(`SortLE(L,Vars): inputs a list of polynomials L in the set of variables Vars, outputs the sorted`):
print(`version in increasing number of variables. Try:`):
print(`SortLE([x+y+z,z,y+z,x+z],{x,y,z});`):


elif nargs=1 and args[1]=SR then
print(`SR(Zug,t,L,X,Y,K): Inputs a pair of rational functions Zug in the variable t, and a list of polynomials in the`):
print(`variables X and Y and a fairly large positive integer K (a parameter for guessing), outputs`):
print(`the generating functions for the members of L of the sequence obtained by replacing X and Y by the Taylor`):
print(`coefficients of Zug[1] and Zug[2] respectively. Try:`):
print(`SR([1/(1-9*t-t^2),t/(1-9*t-t^2)],t,[X^2+7*X*Y-9*Y^2, 2*X^2-4*X*Y+12*Y^2,2*X^2+10*Y^2],X,Y,10);`):
print(`SR([-(t-82)/(t^2-83*t+1), 9/(t^2-83*t+1)],t,[X^2+7*X*Y-9*Y^2, 2*X^2-4*X*Y+12*Y^2,2*X^2+10*Y^2],X,Y,10);`):
print(`SR([(8*t-1)/(t^2+9*t-1), -t/(t^2+9*t-1)],t,[X^2-9*X*Y-Y^2, 2*X^2-4*X*Y+12*Y^2,2*X^2+10*Y^2],X,Y,10);`):



elif nargs=1 and args[1]=SR1 then
print(`SR1(Zug,t,P,X,Y,K): Inputs a pair of rational functions Zug in the variable t, and a polynomial P in the`):
print(`variables X and Y and a fairly large positive integer K (a parameter for guessing), outputs`):
print(`the generating functions for  sequence obtained by replacing, in P,  X and Y by the Taylor`):
print(`coefficients of Zug[1] and Zug[2] respectively. Try:`):
print(`SR1([1/(1-9*t-t^2),t/(1-9*t-t^2)],t,X^2+7*X*Y-9*Y^2,X,Y,10);`):

else
 print(`There is no such thing as`, args):

fi:


end:



###start from Cfinite.txt

#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:

eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:




#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

for d from 1 to trunc(nops(L)/2)-2 do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:


#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:


#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):
L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
#print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(f):
fi:
end:


#RtoC(R,t): Given a rational function R(t) finds the
#C-finite description of its sequence of coefficients
#For example, try:
#RtoC(1/(1-t-t^2),t);
RtoC:=proc(R,t) local S1,S2,D1,R1,d,f1,i,a:
R1:=normal(R):
D1:=denom(R1):
d:=degree(D1,t):
if degree(numer(R1),t)>=d then
print(`The degree of the top must be less than the degree of the bottom`):
 RETURN(FAIL):
fi:

a:=coeff(D1,t,0):
S2:=[seq(-coeff(D1,t,i)/a,i=1..degree(D1,t))]:
f1:=taylor(R,t=0,d+1):
S1:=[seq(coeff(f1,t,i),i=0..d-1)]:
[S1,S2]:
end:

###end from Cfinite.txt





#DE(k,a0,a1,b0,b1,X,Y,t): Finds a polynomial P in the variables X and Y such that
#if a(n) b(n) are solutions of the recurrence
#x(n)=k*x(n-1)-x(n-2) with a(0)=a0, a(1)=a(1), b(0)=b0, b(1)=b1. 
#then P(a(n),b(n))=0
#It also returns the generating functions for a(n) and b(n)
#Try:
#DE(3,1,2,3,4,X,Y,t);

DE:=proc(k,a0,a1,b0,b1,X,Y,t) local gu1,gu2,P,c1,c2,c3,eq,var,i:

gu1:=SeqFromRec([[a0,a1],[k,-1]],15):

gu2:=SeqFromRec([[b0,b1],[k,-1]],15):

P:=c1*X^2+c2*X*Y+c3*Y^2+1:
var:={c1,c2,c3}:

eq:={seq(subs({X=gu1[i],Y=gu2[i]},P),i=1..6)}:

var:=solve(eq,var):

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


P:=numer(subs(var,P)):


if {seq(subs({X=gu1[i],Y=gu2[i]},P),i=1..15)}<>{0} then
 RETURN(FAIL):
fi:

[P,[CtoR([[a0,a1],[k,-1]],t),CtoR([[b0,b1],[k,-1]],t)]]:

end:





#DEm(k,a0,a1,b0,b1,X,Y,t): Finds a polynomial P in the variables X and Y such that
#if a(n) b(n) are solutions of the recurrence
#x(n)=k*x(n-1)+x(n-2) with a(0)=a0, a(1)=a(1), b(0)=b0, b(1)=b1. 
#then P(a(n),b(n))=(-1)^n
#It also returns the generating functions for a(n) and b(n)
#Try:
#DEm(9,0,1,1,9,X,Y,t);
DEm:=proc(k,a0,a1,b0,b1,X,Y,t) local gu1,gu2,P,c1,c2,c3,eq,var,i:

gu1:=SeqFromRec([[a0,a1],[k,1]],15):

gu2:=SeqFromRec([[b0,b1],[k,1]],15):

P:=c1*X^2+c2*X*Y+c3*Y^2:
var:={c1,c2,c3}:

eq:={seq(subs({X=gu1[i],Y=gu2[i]},P)-(-1)^(i-1),i=1..6)}:

var:=solve(eq,var):

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

P:=subs(var,P):


if {seq(subs({X=gu1[i],Y=gu2[i]},P)-(-1)^(i-1),i=1..15)}<>{0} then
 RETURN(FAIL):
fi:

[P,[CtoR([[a0,a1],[k,1]],t),CtoR([[b0,b1],[k,1]],t)]]:

end:




#FindEqs1(L,t,Vars,P,X,Y,K): inputs a pair of rational functions in t phrased in terms of parameters in Vars
#and a  polynomial P in the variables X and Y, finds a Groebner basis for the set of nops(Vars)+K equations 
#such that the  if a(n) and b(n) are the coefficients of t^n of L[1], L[2] respectively then
#P(a(n),b(n))=0. If there is no solution it returns FAIL. Try:
#FindEqs1([(a0+a1*t)/(1-6*t+t^2),(b0+b1*t)/(1-6*t+t^2)],t,[a0,a1,b0,b1],X^2-2*Y^2-1,X,Y,4);
FindEqs1:=proc(L,t,Vars,P,X,Y,K) local eq,gu1,gu2,i:

gu1:=taylor(L[1],t=0,nops(Vars)+K):
gu2:=taylor(L[2],t=0,nops(Vars)+K):

eq:=expand({seq(subs({X=coeff(gu1,t,i),Y=coeff(gu2,t,i)},P),i=1..nops(Vars)+3)}):

Groebner[Basis](eq,plex(op(Vars)) ):
end:





#FindVars(P,Vars): inputs a polynomial in a subset of the set of variables Vars,
#finds the subest that actually participates in it. Try:
#FindVars(x^3+y*z,{x,y,z,w});
FindVars:=proc(P,Vars) local gu,v:

gu:={}:

for v in Vars do
 if expand(diff(P,v))<>0 then
  gu:=gu union {v}:
 fi:
od:

gu:

end:


#SortLE(L,Vars): inputs a list of polynomials L in the set of variables Vars, outputs the sorted
#version in increasing number of variables. Try:
#SortLE([x+y+z,z,y+z,x+z],{x,y,z});
SortLE:=proc(L,Vars) local L1,i,gu,T,gu1:

gu:={seq(nops(FindVars(L[i],Vars)),i=1..nops(L))}:

for gu1 in gu do
 T[gu1]:={}:
od:

for  i from 1 to nops(L) do
 T[nops(FindVars(L[i],Vars))]:=T[nops(FindVars(L[i],Vars))] union {L[i]}:
od:


L1:= sort(convert({seq(nops(FindVars(L[i],Vars)),i=1..nops(L))},list)):


[seq(op(T[L1[i]]),i=1..nops(L1))]:

end:


#SolBF1a(P,Var,M): Give a polynomial P in the list of variables Var finds all solutions of P(Var)=0
#with values from -M to M by brute force. Try:
#SolBF1a(x^2-2*y^2-1,[x,y],M);
SolBF1a:=proc(P,Var,M) local x,Var1,lu,lu1,gu,P1,i:
if Var=[] then
RETURN(FAIL):
fi:

x:=Var[1]:
if nops(Var)=1 then

 gu:={}:

for  i from -M to M do
 if subs(x=i,P)=0 then
  gu:=gu union {[i]}:
 fi:
od:
RETURN(gu)
fi:

Var1:=[op(2..nops(Var),Var)]:

gu:={}:

for i from -M to M do
 P1:=subs(x=i,P):
 lu:=SolBF1a(P1,Var1,M):
 gu:=gu union {seq([i,op(lu1)],lu1 in lu)}:
od:
gu:
end:



#SolBF1(P,Var,M): Give a polynomial P in the list of variables Var finds all solutions of P(Var)=0
#with values from -M to M by brute force. Try:
#SolBF1(x^2-2*y^2-1,[x,y],M);
SolBF1:=proc(P,Var,M) local gu,gu1,i:

gu:=SolBF1a(P,Var,M):

{seq({seq(Var[i]=gu1[i],i=1..nops(Var))},gu1 in gu)}:
end:




#FindEqs(L,t,Vars,P,X,Y): inputs a pair of rational functions in t phrased in terms of parameters in Vars
#and a  polynomial P in the variables X and Y, finds a Groebner basis for the set of nops(Vars)+K equations 
#such that the  if a(n) and b(n) are the coefficients of t^n of L[1], L[2] respectively then
#P(a(n),b(n))=0. If there is no solution it returns FAIL. Try:
#FindEqs([(a0+a1*t)/(1-6*t+t^2),(b0+b1*t)/(1-6*t+t^2)],t,[a0,a1,b0,b1],X^2-2*Y^2-1,X,Y);
FindEqs:=proc(L,t,Vars,P,X,Y) local K,gu1,gu2,gu3:

gu1:=FindEqs1(L,t,Vars,P,X,Y,4):

gu2:=FindEqs1(L,t,Vars,P,X,Y,5):

for K from 6 while gu1<>gu2 do
 gu3:=FindEqs1(L,t,Vars,P,X,Y,K):
 gu1:=gu2:
 gu2:=gu3:
od:


SortLE(gu2,{op(Vars)}):


end:



#SolBF(L,Var,M): inputs a list L of equations in the set of variables Var of variables Var, 
#and  a pos. integer M, finds all solutions in integers from -M to M. Try
#SolBF([x^2-1,x+y-3,x+y+z-10],{x,y,z},5);
SolBF:=proc(L,Var,M) local P,var1,lu,gu,Var1,L1,lu1,ku1,ku:

if Var={} then
if convert(L,set)={0} then
 RETURN({{}}):
else
 RETURN({}):
fi:
fi:

if nops(L)=1 then
P:=L[1]:

if FindVars(P,Var)<>convert(Var,set) then
  RETURN(FAIL):
else

RETURN(SolBF1(P,convert(Var,list),M)):


fi:

fi:

P:=L[1]:

var1:=FindVars(P,Var):

lu:=SolBF1(P,convert(var1,list),M):

Var1:=Var minus var1:

L1:=[op(2..nops(L),L)]:

gu:={}:

for lu1 in lu do
 ku:=SolBF(subs(lu1,L1),Var1,M):
 gu:=gu union {seq(lu1 union ku1,ku1 in ku)}:
od:
gu:


end:




#SolDE2(P,X,Y,t,M): Inputs a polynomial P , with integer coefficients, in the variables X and Y
#and a variable t, and a positive integer M, finds, if it exists
#two rational functions of the form [(a0+a1*t)/(1-k*t+t^2),(b0+b1*t)/(1-k*t+t^2)] with
#k positive and larger than 2, and a0,a1,b0,b1 between -M and M, such that
#all the respective Taylor coefficients satisfy P(X,Y)=0. For example, to solve
#Pell's equation X^2-2*Y^2=1 type:
#SolDE2(X^2-2*Y^2-1,X,Y,t,6);
SolDE2:=proc(P,X,Y,t,M) local a0,a1,b0,b1,k,Zug, mu,Var,lu,lu1,gu,tov,i,gu1:

if degree(P,{X,Y})<>2 then
 print(P, `should be a quadartic polynomial in`, X,Y):
 RETURN(FAIL):
fi:

Zug:=[(a0+a1*t)/(1-k*t+t^2),(b0+b1*t)/(1-k*t+t^2)]:

mu:=FindEqs(Zug,t,[a0,a1,b0,b1,k],P,X,Y);

Var:={a0,a1,b0,b1,k}:

lu:=SolBF(mu,Var,M):


gu:={}:

for lu1 in lu do

if subs(lu1,k)>2 then
 gu:= gu union {normal(subs(lu1,Zug))}:
fi:

od:

gu:

tov:={}:

for gu1 in gu do
 if min(seq(coeff(taylor(gu1[1],t=0,11),t,i),i=1..10))>0  and
     min(seq(coeff(taylor(gu1[2],t=0,11),t,i),i=1..10))>0 then
  tov:=tov union {gu1}:
 fi:
od:

if tov<>{} then
 gu:=tov:
fi:

gu: 

end:





####start with (-1)^n




#FindEqsM1(L,t,Vars,P,X,Y,K): inputs a pair of rational functions in t phrased in terms of parameters in Vars
#and a  polynomial P in the variables X and Y, finds a Groebner basis for the set of nops(Vars)+K equations 
#such that the  if a(n) and b(n) are the coefficients of t^n of L[1], L[2] respectively then
#P(a(n),b(n))=(-1)^n. If there is no solution it returns FAIL. Try:
#FindEqsM1([(a0+a1*t)/(1-9*t-t^2),(b0+b1*t)/(1-9*t-t^2)],t,[a0,a1,b0,b1],X^2-9*X*Y-Y^2,X,Y,5);
FindEqsM1:=proc(L,t,Vars,P,X,Y,K) local eq,gu1,gu2,i:

gu1:=taylor(L[1],t=0,nops(Vars)+K):
gu2:=taylor(L[2],t=0,nops(Vars)+K):

eq:=expand({seq(subs({X=coeff(gu1,t,i),Y=coeff(gu2,t,i)},P)-(-1)^i,i=1..nops(Vars)+3)}):

Groebner[Basis](eq,plex(op(Vars)) ):
end:


#FindEqsM(L,t,Vars,P,X,Y): inputs a pair of rational functions in t phrased in terms of parameters in Vars
#and a  polynomial P in the variables X and Y, finds a Groebner basis for the set of nops(Vars)+K equations 
#such that the  if a(n) and b(n) are the coefficients of t^n of L[1], L[2] respectively then
#P(a(n),b(n))=(-1)^n. If there is no solution it returns FAIL. Try:
#FindEqsM([(a0+a1*t)/(1-9*t-t^2),(b0+b1*t)/(1-9*t-t^2)],t,[a0,a1,b0,b1],X^2-9*X*Y-Y^2,X,Y);
FindEqsM:=proc(L,t,Vars,P,X,Y) local K,gu1,gu2,gu3:

gu1:=FindEqsM1(L,t,Vars,P,X,Y,4):

gu2:=FindEqsM1(L,t,Vars,P,X,Y,5):

for K from 6 while gu1<>gu2 do
 gu3:=FindEqsM1(L,t,Vars,P,X,Y,K):
 gu1:=gu2:
 gu2:=gu3:
od:


SortLE(gu2,{op(Vars)}):


end:





#SolDE2m(P,X,Y,t,M): Inputs a polynomial P , with integer coefficients, in the variables X and Y
#and a variable t, and a positive integer M, finds, if it exists
#two rational functions of the form [(a0+a1*t)/(1-k*t+t^2),(b0+b1*t)/(1-k*t+t^2)] with
#k positive and larger than 2, and a0,a1,b0,b1 between -M and M, such that
#all the respective Taylor coefficients satisfy P(a(n),b(n))=(-1)^n. For example, try
#SolDE2m(X^2-9*X*Y-Y^2,X,Y,t,6);
SolDE2m:=proc(P,X,Y,t,M) local a0,a1,b0,b1,k,Zug, mu,Var,lu,lu1,gu,tov,i,gu1:


if degree(P,{X,Y})<>2 then
 print(P, `should be a quadartic polynomial in`, X,Y):
 RETURN(FAIL):
fi:

Zug:=[(a0+a1*t)/(1-k*t-t^2),(b0+b1*t)/(1-k*t-t^2)]:

mu:=FindEqsM(Zug,t,[a0,a1,b0,b1,k],P,X,Y);

Var:={a0,a1,b0,b1,k}:

lu:=SolBF(mu,Var,M):


gu:={}:

for lu1 in lu do

if subs(lu1,k)>2 then
 gu:= gu union {normal(subs(lu1,Zug))}:
fi:

od:

gu:

tov:={}:

for gu1 in gu do
 if min(seq(coeff(taylor(gu1[1],t=0,11),t,i),i=1..10))>0  and
     min(seq(coeff(taylor(gu1[2],t=0,11),t,i),i=1..10))>0 then
  tov:=tov union {gu1}:
 fi:
od:

if tov<>{} then
 gu:=tov:
fi:

gu: 

end:



#SR1(Zug,t,P,X,Y,K): Inputs a pair of rational functions Zug in the variable t, and a polynomial P in the
#variables X and Y and a fairly large positive integer K (a parameter for guessing), outputs
#the generating functions for  sequence obtained by replacing, in P,  X and Y by the Taylor
#coefficients of Zug[1] and Zug[2] respectively. Try:
#SR1([1/(1-9*t-t^2),t/(1-9*t-t^2)],t,X^2+7*X*Y-9*Y^2,X,Y,10);
SR1:=proc(Zug,t,P,X,Y,K) local gu1,gu2,gu,i:
gu1:=taylor(Zug[1],t=0,K+1):
gu2:=taylor(Zug[2],t=0,K+1):

gu1:=[seq(coeff(gu1,t,i),i=0..K)]:
gu2:=[seq(coeff(gu2,t,i),i=0..K)]:
gu:=[seq(subs({X=gu1[i],Y=gu2[i]},P),i=1..K+1)]:
gu:=GuessRec(gu) :

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

factor(CtoR(gu,t)):

end:




#SR(Zug,t,L,X,Y,K): Inputs a pair of rational functions Zug in the variable t, and a list of polynomials in the
#variables X and Y and a fairly large positive integer K (a parameter for guessing), outputs
#the generating functions for the members of L of the sequence obtained by replacing X and Y by the Taylor
#coefficients of Zug[1] and Zug[2] respectively. Try:
#SR([1/(1-9*t-t^2),t/(1-9*t-t^2)],t,[X^2+7*X*Y-9*Y^2, 2*X^2-4*X*Y+12*Y^2,2*X^2+10*Y^2],X,Y,10);
SR:=proc(Zug,t,L,X,Y,K) local i,gu,mu:

gu:=[]:

for i from 1 to nops(L) do
 mu:=SR1(Zug,t,L[i],X,Y,K):
  if mu=FAIL then
    RETURN(FAIL):
  else
  gu:=[op(gu),mu]:
 fi:
od:

gu:

end:

####start General









#GenPolHst(m,k,d,c,st): Inputs a symbol m, a pos. integer k, a pos. integer d, and a symbol c, outputs
#A generic homog. polynomial in m[1], ..,m[k] with indexed coefficients c[i] starting at i=st. 
#It also outputs the set of coefficients. Try:
#GenPolHst(m,2,2,c,1);
GenPolHst:=proc(m,k,d,c,st) local gu,kv,i,lu,st1:

if k=1 then
 RETURN(c[st]*m[1]^d, {c[st]}):
fi:

gu:=0:
kv:={}:

st1:=st:

for i from 0 to d do
 lu:=GenPolHst(m,k-1,d-i,c,st1):
 gu:=expand(gu+lu[1]*m[k]^i):
 kv:=kv union lu[2]:
 st1:=st1+nops(lu[2]):
od:
gu,kv:
end:

#GenPol(m,k,d,c): Inputs a symbol m, a pos. integer k, a pos. integer d, and a symbol c, outputs
#A generic polynomial in m[1], ..,m[k] with indexed coefficients c[i] and starting at i=st. 
#It also outputs the set of coefficients. Try:
#GenPol(m,2,2,c);
GenPol:=proc(m,k,d,c) local gu,kv,d1,mu,st:

gu:=1:
kv:={}:
st:=1:

for  d1 from 1 to d do

 mu:=GenPolHst(m,k,d1,c,st):

 st:=st+nops(mu[2]):

 gu:=gu+mu[1]:
 kv:=kv union mu[2]:

od:

gu,kv:

end:





#DEg1(L,t,X,d): inputs a list L of rational functions in t of length k say in the varible t,
#a letter X and a pos. integer d,outputs a polynomial in X[1], ..., X[k] of degree d 
#let's call it P(X[1], ..., X[k])
#such
#that if you plug-in the respective Taylor coefficients of L[i] are solutions of the
#Diphantine equation P(X[1], ..., X[k])=0.
#Try:
#DEg1([(1+t)/(1-3*t+t^2],4/(1-3*t+t^2],t,X,2);

DEg1:=proc(L,t,X,d) local k,P,K,var,eq,c,i,j:

k:=nops(L):

P:=GenPol(X,k,d,c);
var:=P[2]:
P:=P[1]:

K:=nops(var)+6:

eq:={seq(subs({seq(X[i]=coeff(taylor(L[i],t=0,K+1),t,j),i=1..k)},P),j=0..K)}:

var:=solve(eq,var):

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


P:=numer(subs(var,P)):


end:



#DEg(L,t,X,d): inputs a list L of rational functions in t of length k say in the varible t,
#a letter X and a pos. integer d,outputs a polynomial in X[1], ..., X[k] of degree <=d 
#let's call it P(X[1], ..., X[k])
#DEg([(1+t)/(1-3*t+t^2],4/(1-3*t+t^2],t,X,2);

DEg:=proc(L,t,X,d) local d1,gu:

for d1 from 1 to d do
 gu:=DEg1(L,t,X,d1):

 if gu<>FAIL then
   RETURN(gu):
 fi:
od:

FAIL:

end:






#DEgO(L,t,X,d): inputs a list L of rational functions in t of length k say in the varible t,
#a letter X and a pos. integer d,outputs a polynomial in X[1], ..., X[k] of homog. degree d  plust a constant
#let's call it P(X[1], ..., X[k])
#such
#that if you plug-in the respective Taylor coefficients of L[i] are solutions of the
#Diphantine equation P(X[1], ..., X[k])=1.
#Try:
#DEgO([(1+t)/(1-3*t+t^2],4/(1-3*t+t^2],t,X,2);

DEgO:=proc(L,t,X,d) local k,P,K,var,eq,c,i,j:

k:=nops(L):

P:=GenPolHst(X,k,d,c,1);
var:=P[2]:
P:=P[1]:

K:=nops(var)+6:

eq:={seq(subs({seq(X[i]=coeff(taylor(L[i],t=0,K+1),t,j),i=1..k)},P)-1,j=0..K)}:

var:=solve(eq,var):

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


P:=numer(subs(var,P)):


end:


#DEgM(L,t,X,d): inputs a list L of rational functions in t of length k say in the varible t,
#a letter X and a pos. integer d,outputs a polynomial in X[1], ..., X[k] of homog. degree d  
#let's call it P(X[1], ..., X[k])
#such
#that if you plug-in the respective Taylor coefficients of L[i] are solutions of the
#Diphantine equation P(X[1], ..., X[k])=(-1)^k.
#Try:
#DEgM([(1+t)/(1-3*t+t^2],4/(1-3*t+t^2],t,X,2);

DEgM:=proc(L,t,X,d) local k,P,K,var,eq,c,i,j:

k:=nops(L):

P:=GenPolHst(X,k,d,c,1);
var:=P[2]:
P:=P[1]:

K:=nops(var)+6:

eq:={seq(subs({seq(X[i]=coeff(taylor(L[i],t=0,K+1),t,j),i=1..k)},P)-(-1)^j,j=0..K)}:

var:=solve(eq,var):

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


P:=subs(var,P):
#P:=numer(subs(var,P)):


end:


#######################End  stuff from RatDio.txt


#FindDE(L,m,X): inputs a list L of polynomials in m[1], ..., m[k-1] (where k=nops(L),  outputs the minimal polynomial
#in X[1], ... X[k] such that
#P(X[1],..., X[k]]) such that P(L[1], ..., L[m]) is identically zero. Try:
#FindDE([m[1]^2+1,m[1]^3+2],m,X);
#FindDE([m[1]^2-m[2]^2,2*m[1]*m[2],m[1]^2+m[2]^2],m,X);
FindDE:=proc(L,m,X) local k,eq,var,i:
k:=nops(L):
var:=[seq(m[i],i=1..k-1),seq(X[i],i=1..k)]:
eq:={seq(X[i]-L[i],i=1..k)}:
Groebner[Basis](eq,plex(op(var)))[1]:
end:

#GenPol1(m,k,d,c,st): Inputs a symbol m, a pos. integer k, a pos. integer d, and a symbol c, outputs
#A generic homog. polynomial in m[1], ..,m[k] with indexed coefficients c[i] starting at i=st. 
#It also outputs the set of coefficients. Try:
#GenPol1(m,2,2,c,1);
GenPol1:=proc(m,k,d,c,st) local gu,kv,i,lu,st1:

if k=1 then
 RETURN(c[st]*m[1]^d, {c[st]}):
fi:

gu:=0:
kv:={}:

st1:=st:

for i from 0 to d do
 lu:=GenPol1(m,k-1,d-i,c,st1):
 gu:=expand(gu+lu[1]*m[k]^i):
 kv:=kv union lu[2]:
 st1:=st1+nops(lu[2]):
od:
gu,kv:
end:



#ExtractCoeffsOld(P,m,k): inputs a polynomial P in m[1], ..., m[k] outputs the set of coefficients of all monomials. Try:
#ExtractCoeffsOld(5*m[1]*m[2]+3*m[1]^2,m,2);
ExtractCoeffsOld:=proc(P,m,k) local gu,i:

if k=1 then
 RETURN({seq(coeff(P,m[1],i), i=0..degree(P,m[1]))}):
fi:

gu:={}:

for i from 0 to degree(P,m[k]) do
 gu:=gu union  ExtractCoeffsOld(coeff(P,m[k],i),m,k-1):
od:

gu minus {0}:

end:

#FindEqsPold(P,X,r,m,k,d,c): Inputs
#(i)a polynomial P in X[1], .., X[r], (where X is a symbol)
#(ii) a symbol m
#(iii) pos. integers k and d
#(iv) a symbol c. Outputs
#(i) a list of length r with tenative expressions of generic polynomials with indexed variables given by c[i]'s
#let's call it L
#(ii) The set of equations in the c[i]'s such that P(L[1], ..., L[r]) is identically zero
#(iii) the set of c[i]'2. Try:
#FindEqsPold(X[1]^2+X[2]^2-X[3]^2,X,3,m,2,2,c);
#FindEqsPold(X[1]^3+X[2]^3+X[3]^3+X[4]^3,X,4,m,2,2,c);
FindEqsPold:=proc(P,X,r,m,k,d,c) local L,st,i,lu,var,gu,eq:
st:=1:
L:=[]:
var:={}:
for i from 1 to r do
 lu:=GenPol1(m,k,d,c,st):
 L:=[op(L),lu[1]]:
 var:=var union lu[2]:
 st:=st+nops(lu[2]):
od:

gu:=expand(subs({seq(X[i]=L[i],i=1..r)},P));

eq:=ExtractCoeffsOld(gu,m,k):
eq:=convert(eq,list):

eq:=SortLE(eq,var):
eq,var,L:


end:




#SolDEPold(P,X,r,m,k,d,M): Inputs
#(i)a polynomial P in X[1], .., X[r], (where X is a symbol)
#(ii) a symbol m
#(iii) pos. integers k and d. Outputs a set of parametric solutions to
#the diophantine equation P(X[1], ..., X[r])=0 where each solution is a list
#of length r whose entries are  homog. polynomials in m[1],...,m[k] of total degree d, and the
#coefficients are between -M and M, excluding trivial solutions.
#a list of length r with proiv
#Try:
#SolDEPold(X[1]^2+X[2]^2-X[3]^2,X,3,m,2,2,1);
#SolDEPold(X[1]^3+X[2]^3+X[3]^3+X[4]^3,X,4,m,2,2,6);
SolDEPold:=proc(P,X,r,m,k,d,M) local gu,c,var,L,eq,lu,lu1,mu:
gu:=FindEqsPold(P,X,r,m,k,d,c):
eq:=gu[1]:
var:=gu[2]:
L:=gu[3]:

lu:=SolBF(eq,var,M):

gu:={}:

for lu1 in lu do

mu:=subs(lu1,L):

if not member(0,convert(mu,set)) then
 gu:=gu union {mu}:
fi:

od:

gu:

end:


#SolBF1a(P,Var,M): Give a polynomial P in the list of variables Var finds all solutions of P(Var)=0
#with values from -M to M by brute force. Try:
#SolBF1a(x^2-2*y^2-1,[x,y],M);
SolBF1a:=proc(P,Var,M) local x,Var1,lu,lu1,gu,P1,i:
if Var=[] then
#RETURN(FAIL):
RETURN({}):
fi:

x:=Var[1]:
if nops(Var)=1 then

 gu:={}:

for  i from -M to M do
 if subs(x=i,P)=0 then
  gu:=gu union {[i]}:
 fi:
od:
RETURN(gu)
fi:

Var1:=[op(2..nops(Var),Var)]:

gu:={}:

for i from -M to M do
 P1:=subs(x=i,P):
 lu:=SolBF1a(P1,Var1,M):
 gu:=gu union {seq([i,op(lu1)],lu1 in lu)}:
od:
gu:
end:



#SolBF1(P,Var,M): Give a polynomial P in the list of variables Var finds all solutions of P(Var)=0
#with values from -M to M by brute force. Try:
#SolBF1(x^2-2*y^2-1,[x,y],M);
SolBF1:=proc(P,Var,M) local gu,gu1,i:

gu:=SolBF1a(P,Var,M):

{seq({seq(Var[i]=gu1[i],i=1..nops(Var))},gu1 in gu)}:
end:




#FindEqs(L,t,Vars,P,X,Y): inputs a pair of rational functions in t phrased in terms of parameters in Vars
#and a  polynomial P in the variables X and Y, finds a Groebner basis for the set of nops(Vars)+K equations 
#such that the  if a(n) and b(n) are the coefficients of t^n of L[1], L[2] respectively then
#P(a(n),b(n))=0. If there is no solution it returns FAIL. Try:
#FindEqs([(a0+a1*t)/(1-6*t+t^2),(b0+b1*t)/(1-6*t+t^2)],t,[a0,a1,b0,b1],X^2-2*Y^2-1,X,Y);
FindEqs:=proc(L,t,Vars,P,X,Y) local K,gu1,gu2,gu3:

gu1:=FindEqs1(L,t,Vars,P,X,Y,4):

gu2:=FindEqs1(L,t,Vars,P,X,Y,5):

for K from 6 while gu1<>gu2 do
 gu3:=FindEqs1(L,t,Vars,P,X,Y,K):
 gu1:=gu2:
 gu2:=gu3:
od:


SortLE(gu2,{op(Vars)}):


end:



####start new stuff July 1, 2020



#ExtractCoeffs(P,v): inputs a polynomial P in the list of variables v, outputs the set of coefficients of all monomials. Try:
#ExtractCoeffs(5*m1*m2+3*m1^2,[m1,m2]);
ExtractCoeffs:=proc(P,v) local k,gu,i,v1:

k:=nops(v):
if k=1 then
 RETURN({seq(coeff(P,v[1],i), i=0..degree(P,v[1]))}):
fi:

gu:={}:

v1:=[op(1..k-1,v)]:

for i from 0 to degree(P,v[k]) do
 gu:=gu union  ExtractCoeffs(coeff(P,v[k],i),v1):
od:

gu minus {0}:

end:

#FindEqsP(P,x,L,v,Pars): Inputs
#(i)a polynomial P in the list of variables x
#(ii) a list L of the same length as x  of the variables v the coefficients may involve parameters.
#(iii) the set of parameters
#Outputs:
#The set of coefficients in the variables in v of all monomials. Try:
#FindEqsP(X^2+Y^2-Z^2,[X,Y,Z],[a20*m^2+a11*m*n+a02*n^2,b20*m^2+b11*m*n+b02*n^2,c20*m^2+c11*m*n+c02*n^2],[m,n],{a20,a11,a02,b20,b11,b02,c20,c11,c02});
FindEqsP:=proc(P,x,L,v,Pars) local gu,eq,i:

if not (type(x,list) and type(L,list) and nops(x)=nops(L) and type(v,list) and type(Pars,set)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=expand(subs({seq(x[i]=L[i],i=1..nops(x))},P));

eq:=ExtractCoeffs(gu,v):
eq:=convert(eq,list):

eq:=SortLE(eq,v):
eq:

end:



#IsVeryTrivial(L): decides whether a solution is trivial
IsVeryTrivial:=proc(L) local i,j:

 if member(0,convert(L,set)) then
   RETURN(true):
 fi:

if igcd(op(L))<>1 then
 RETURN(true):
fi:


false:

end:


#CleanUp1(S,var,i): Given a set of polynomial tuples in the list of variables var and i between 1 and nops(var)
#kicks out those obtained by var[i]->-var[i]. try:
#CleanUp1({m*n,-m*n},[m,n],1);
CleanUp1:=proc(S,var,i) local S1,S2,s,s1:
S1:={}:
S2:=S:

 while S2<>{} do
   s:=S2[1]:
   s1:=subs(var[i]=-var[i],s):
 if member(s1,S2) then
   S1:=S1 union {s}:
    S2:=S2 minus {s1,s}:
 fi:
od:

S1:

end:

#CleanUp(S,var): Given a list of tuples of polynomials in the list of variables var, kicks out redundanies, Try:
#CleanUp([m,-m],[m]); 
CleanUp:=proc(S,var) local i,S1:

S1:=S:
for i from 1 to nops(var) do
 S1:=CleanUp1(S,var,i):
od:
S1:
end:



#IsTrivial(L): decides whether a solution is trivial
IsTrivial:=proc(L) local i:

 if member(0,convert(L,set)) then
   RETURN(true):
 fi:

if {seq(type(normal(L[i]/L[1]),numeric), i=2..nops(L))}={true} then
   RETURN(true):
 fi:

false:

end:

#SolDEP(P,x,L,v,Pars,M): Inputs
#(i)a polynomial P in the list of variables x
#(ii) a list L of the same length as x  of the variables v the coefficients may involve parameters.
#(iii) the set of parameters
#(iv) a positive integer M
#Outputs:
#The set of specializations of L with integers between -M and M such that P(L[1], ..., L[nops(L)])=0. Try:
#SolDEP(X^2+Y^2-Z^2,[X,Y,Z],[m^2+a11*m*n+a02*n^2,b20*m^2+b11*m*n+b02*n^2,c20*m^2+c11*m*n+c02*n^2],[m,n],{a11,a02,b20,b11,b02,c20,c11,c02},2);
#SolDEP(X^3+Y^3+Z^3+W^3,[X,Y,Z,W],
#[m^2+a11*m*n+a02*n^2,2*m^2+b11*m*n+b02*n^2,-2*m^2+c11*m*n+c02*n^2,-m^2+d11*m*n+d02*n^2],[m,n].{a11,a02,b20,b11,b02,c11,c02,d11,d02},12);

SolDEP:=proc(P,x,L,v,Pars,M) local eq,lu,lu1,gu,gu1:
eq:=FindEqsP(P,x,L,v,Pars):

lu:=SolBF(eq,Pars,M):

gu:={}:

for lu1 in lu do
 gu1:=subs(lu1,L):
  if not IsTrivial(gu1) then
   gu:=gu union {gu1}:
  fi:
od:

CleanUp(gu,v):

end:



#RandPolH(m,k,d,K): A random polynomial in m[1], ..., m[k], of homeg. degree k with coefficients from -K to K. Try:
#RandPolH(m,2,2,3);
RandPolH:=proc(m,k,d,K) local ra,gu,i:

ra:=rand(-K..K):

if k=1 then
 RETURN(ra()*m[1]^d):
fi:

gu:=0:

for i from 0 to d do
 gu:=gu+RandPolH(m,k-1,d-i,K)*m[k]^i:
od:

gu:

end:


#Rootk(n,k): the cubic-root of n or FAIL. Try:
#Rootk(-8,3);
Rootk:=proc(n,k) local m:

if n=0 then
 RETURN(0):
fi:

if not type(n,integer) then
  RETURN(FAIL):
fi:

if k mod 2=1 then
m:=simplify((abs(n))^(1/k)):

if not type(m,integer) then
 RETURN(FAIL):
fi:


if n>0 then
 RETURN(m):
else
 RETURN(-m):
fi:

else

m:=simplify(n^(1/k)):

if not type(m,integer) then
 RETURN(FAIL):
else
 RETURN(m):
fi:

fi:
end:




#Box1(L): All the vectors x[1], .., x[nops(L)] such that  x[i] is in L[i]. Try:
#Box1([{-1,1},{-2,2}]);
Box1:=proc(L) local k,gu1,gu,L1,i:

k:=nops(L):

if k=1 then
 RETURN({seq([i], i in L[1])}):
fi:

L1:=[op(1..k-1,L)]:
gu:=Box1(L1):

{seq(seq([op(gu1),i],i in L[k]),gu1 in gu)}:

end:



#BFd1(A,d,K): inputs a list A of integers of length k say, and an odd positive ineger d
#finds all solutions of
#add(A[i]*X[i]^d,i=1..k)=0 with X[i], i from 1 to k-1 from 1 to K. Try:
#BFd1([1,1,1,1],3,5);
BFd1:=proc(A,d,K) local k,mu,mu1,gu,hal,i,S:

if not type(d,integer)  then
 print(d, `must have been an odd integer `):
 RETURN(FAIL):
fi:

k:=nops(A):

S:={seq(i,i=1..K)}:

mu:=Box1([S$(k-1)]):

gu:={}:

for mu1 in mu do
 hal:=-add(A[i]*mu1[i]^d,i=1..k-1)/A[k]:
 hal:=Rootk(hal,d):
  if hal<>FAIL then
    gu:=gu union {[op(mu1),hal]}:
  fi:
od:

gu:
 

end:

#Morph43(A,L1,L2,m,n,M): Given two specific solutions of add(A[i]*X[i]^3,i=1..4)=0, tries to find
#a doubly-infinite set of solutions where X[i] are quadratic polynomials in m and n,and it tries all the parameters from -M to M.  Try:
#Morph43([1,1,1,1],[3,4,5,-6],[-5,6,-3,-4],m,n,5):
Morph43:=proc(A,L1,L2,m,n,M) local c,P,i,L,X:

if  not (type(A,list) and nops(A)=4 and type(L1,list) and type(L2,list) and nops(L1)=4 and nops(L2)=4) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if add(A[i]*L1[i]^3,i=1..nops(L1))<>0 then
 print(L1 , `is not a solution of add(A[i]*X[i]^3, i=1..4)=0 `):
 RETURN(FAIL):
fi:


if add(A[i]*L2[i]^3,i=1..nops(L2))<>0 then
 print(L2 , `is not a solution of add(A[i]*X[i]^3, i=1..4)=0 `):
 RETURN(FAIL):
fi:

for i from 1 to 4 do
P[i]:=L1[i]*m^2+c[i]*m*n+L2[i]*n^2:
od:

L:=[seq(P[i],i=1..4)]:

SolDEP(add(A[i]*X[i]^3,i=1..4),[seq(X[i],i=1..4)],L,[m,n],{seq(c[i],i=1..4)},M):

end:






#BFdG(A,d,K): inputs a list A of integers of length k say, and an odd positive ineger d
#finds all solutions of
#add(A[i]*X[i]^d,i=1..k)=0 with X[i], i from 1 to k-1 from -K to K. Try:
#BFdG([1,1,1,1],3,5);
BFdG:=proc(A,d,K) local k,mu,mu1,gu,hal,i,S:

if not type(d,integer)  then
 print(d, `must have been an odd integer `):
 RETURN(FAIL):
fi:

k:=nops(A):
S:={seq(i,i=-K..K)} minus {0}:

mu:=Box1([S$(k-1)]):

gu:={}:

for mu1 in mu do
 hal:=-add(A[i]*mu1[i]^d,i=1..k-1)/A[k]:
 hal:=Rootk(hal,d):
  if hal<>FAIL then
  if not IsVeryTrivial([op(mu1),hal]) then
    gu:=gu union {[op(mu1),hal]}:
  fi:
  fi:
od:

gu:
 

end:

#EriJ(E,L1,L2): uses Eri Jabotinsky's method to get yet-anothe sultion of
#E[1]*X1^3+E[2]*X2^3+E[3]*X3^3+E[4]*X4^3=0 from two given solutions.
#Try:
#EriJ([1,1,1,1],[3,4,5,-6],[m,-m,n,-n]);
EriJ:=proc(E,L1,L2) local i,x,y,X,M:

if not (type(E,list) and type(L1,list) and type(L2,list) and nops(E)=4 and nops(L1)=4 and nops(L2)=4) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if add(E[i]*L1[i]^3,i=1..nops(E))<>0 then
 print(L1, `is not a solution of `, add(E[i]*X[i]^3,i=1..4)=0 ):
 RETURN(FAIL):
fi:


if add(E[i]*L2[i]^3,i=1..nops(E))<>0 then
 print(L2, `is not a solution of `, add(E[i]*X[i]^3,i=1..4)=0 ):
 RETURN(FAIL):
fi:

x:=-add(E[i]*L1[i]*L2[i]^2,i=1..nops(E)):
y:=add(E[i]*L1[i]^2*L2[i],i=1..nops(E)):

M:=normal([seq(x*L1[i]+y*L2[i],i=1..4)]):

if expand(add(E[i]*M[i]^3,i=1..nops(E)))<>0 then
 lprint(M, `did not work out`):
 RETURN(FAIL):
fi:

M:
end:



#EriJs(E,L1,L2,m,n): uses Eri Jabotinsky's method to get yet-anothe sultion of
#E[1]*X1^3+E[2]*X2^3+E[3]*X3^3+E[4]*X4^3=0 from two given solutions.
#here L1 and L2 are two known solutions possibly in terms of symbols m and n
#Try:
#EriJs([1,1,1,1],[3,4,5,-6],[m,-m,n,-n],m,n);
EriJs:=proc(E,L1,L2,m,n) local i,x,y,X,M,c:

if not (type(E,list) and type(L1,list) and type(L2,list) and nops(E)=4 and nops(L1)=4 and nops(L2)=4) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

if add(E[i]*L1[i]^3,i=1..nops(E))<>0 then
 print(L1, `is not a solution of `, add(E[i]*X[i]^3,i=1..4)=0 ):
 RETURN(FAIL):
fi:


if add(E[i]*L2[i]^3,i=1..nops(E))<>0 then
 print(L2, `is not a solution of `, add(E[i]*X[i]^3,i=1..4)=0 ):
 RETURN(FAIL):
fi:

x:=-add(E[i]*L1[i]*L2[i]^2,i=1..nops(E)):
y:=add(E[i]*L1[i]^2*L2[i],i=1..nops(E)):

M:=normal([seq(x*L1[i]+y*L2[i],i=1..4)]):

if expand(add(E[i]*M[i]^3,i=1..nops(E)))<>0 then
 lprint(M, `did not work out`):
 RETURN(FAIL):
fi:

c:=igcd(seq(op([coeff(M[i],m,2), coeff(coeff(M[i],m,1),n,1), coeff(M[i],n,2)]),i=1..nops(M))):

[seq(M[i]/c,i=1..nops(M))]:

end:


#AllSols(P,m,n,K): inputs a polynomial P in m and n and a positive integer K outputs all
#solutions in -K<=m,n<=K of P(m,n)=0. Try:
#AllSols(m^2-2*n^2-1,m,n,30);
AllSols:=proc(P,m,n,K) local gu, m1,n1:

gu:={}:

for m1 from -K to K do
 for n1 from -K to K do
  if subs({m=m1,n=n1},P)=0 then
   gu:=gu union {[m1,n1]}:
  fi:
 od:
od:

gu:

end:


#AllSolsP(P,m,n,K): inputs a polynomial P in m and n and a positive integer K outputs all
#solutions in 0<=m,n<=K of P(m,n)=0. Try:
#AllSolsP(m^2-2*n^2-1,m,n,30);
AllSolsP:=proc(P,m,n,K) local gu, m1,n1:

gu:={}:

for m1 from 0 to K do
 for n1 from 0 to K do
  if subs({m=m1,n=n1},P)=0 then
   gu:=gu union {[m1,n1]}:
  fi:
 od:
od:

gu:

end:



#FindGoodR(m,n,K): finds all the solutions to X1^3+X2^3+X3^3+X4^3 in terms of (m,n) that
#yields one solution in positive integers in m,n<=K to  equalling 1 or -1 for one of its componets.
#Try:
#FindGoodR(m,n,K);
FindGoodR:=proc(m,n,K) local gu,mu,lu,lu1,mu1,i,hal:

lu:=permute([3,4,5,-6]):
mu:=permute([m,-m,n,-n]):

gu:={}:

for lu1 in lu do
 for mu1 in mu do
  hal:=EriJs([1,1,1,1],lu1,mu1,m,n):
   for i from 1 to nops(hal) do
    if AllSolsP(hal[i]-1,m,n,K)<>{} or AllSolsP(hal[i]+1,m,n,K)<>{}  then
      gu:=gu union {[lu1,mu1,hal[i]]}:
    fi:
   od:
 od:
od:

gu:

end:




#Find43(A,L1,L2,m,n,M): Inputs a list of intgeres A of  length 4 and two numerical solutions L1, L2  to
#add(A[i]*X[i]^3,i=1..4)=0 and a pos. integer M, finds all the good Morph43 of L1 and a permutation of L2
#Morph43(A,L,M,m,n,M) (q.d.) where M is  a permutation of -L. Try:
#Find43([1,1,1,1],[3,4,5,-6],[-3,-4,-5,6],m,n,5);
Find43:=proc(A,L,m,n,M) local i,gu,mu,mu1,L1:

L1:=[seq(-L[i],i=1..nops(L))]:

mu:=permute(L1):

gu:={}:

for mu1 in mu do

gu:=gu union Morph43(A,L,mu1,m,n,M):

od:

gu:

end:


#Katan(P,m,n,K): inputs a polynomial P in m and n finds the minimum value (except 0),  in absolute value
#in [-K,K]^2 and where it is and whether it is posive or negative. Try:
#Katan(m^2-2*n^2,m,n,10);
Katan:=proc(P,m,n,K) local ka,m1,n1,gu:

ka:=min({seq(seq(abs(subs({m=m1,n=n1},P)),m1=-K..K),n1=-K..K)} minus {0}):

gu:=AllSols(P-ka,m,n,K):

if gu<>{} then
 RETURN(ka,gu):
fi:

gu:=AllSols(P+ka,m,n,K):

 RETURN(-ka,gu):

end:










#BFd(A,d,K): inputs a list A of integers of length k say, and an odd positive ineger d
#finds all good solutions of
#add(A[i]*X[i]^d,i=1..k)=0 with X[i], i from 1 to k-1 from 1 to K. Try:
#BFd([1,1,1,1],3,5);
BFd:=proc(A,d,K) local gu,gu1,mu,mu1,i:
gu:=BFd1(A,d,K):

gu:=gu union {seq(-gu1,gu1 in gu)}:

mu:={seq(op(permute(gu1)),gu1 in gu)}:

mu:
gu:={}:

for mu1 in mu do

if expand(add(A[i]*mu1[i]^3,i=1..nops(mu1)))=0 then
 mu1:=mu1/igcd(op(mu1)):
 gu:=gu union {mu1}:
fi:
od:

gu:


end:


#Cands(A,m,n,K,r): inputs a list of integers A of length 4 and symbols m,n,and a pos. integer K, and r
#outputs a set of size <=r of quadruples of quadratic polynomials
#such that for each member L such that add(A[i]*L[i]^3,i=1..4)=0 
#the empty set. K is the trial parameter. Try:
#Cands([1,1,1,1],m,n,5,1);
Cands:=proc(A,m,n,K,r) local gu,i,j,hal,S:

gu:=BFd(A,3,K):

S:={}:

for i from 1 to nops(gu) while nops(S)<=r do

 for j from 1 to nops(gu) while nops(S)<=r do

  hal:=Morph43(A,gu[i],gu[j],m,n,K):
  if hal<>{} then
   S:=S union hal:
  fi:

 od:

od:

S:
end:


#ParAB3(A,B,m,n,K): inputs positive integers A and B and outputs a parametric equation, in terms of quadratic polynomials in m and n
#to the Diphantine equation
#A*X^3+A*Y^3+B*Z^3+B*W^3=0 using Eri Jabotinsky's "trick" by first guessing a particular solution of size <=K. If
#none found it returns FAIL. Try:
#ParAB3(2,3,m,n,8);
ParAB3:=proc(A,B,m,n,K) local gu1,gu:
gu1:=BFd([A,A,B,B],3,K):

if gu1={} then
  RETURN(FAIL):
fi:


gu1:=gu1[1]:

gu:=EriJs([A,A,B,B],gu1,[m,-m,n,-n],m,n);

if expand(A*(gu[1]^3+gu[2]^3)+B*(gu[3]^3+gu[4]^3))<>0 then
 print(`Something went wrong`):
 RETURN(FAIL):
fi:

[gu1,gu]:

end:


#Reject(P,var,N): Given a  polynomial with integer coefficients in the set of variables var and a positive integer N
#decides whether there exists a prime power less than N such that if you solve it mod this power there is no solution of
#P=0 and hence there is no solution. try:
#Reject(19*X^2-Y^2-2,[X,Y],30);
Reject:=proc(P,var,N) local gu,gu1,n,p,e,i,j:

p:=2:

while p<=N do
 for e from 1 while p^e<=N do
  n:=p^e:
  gu:=Box1([{seq(i,i=0..n-1)}$nops(var)]):
  if not member(0,  {seq(subs({seq(var[j]=gu1[j],j=1..nops(gu1))},P) mod n,gu1 in gu)}) then
    RETURN(n):
  fi:
 od:
p:=nextprime(p):
od:
FAIL:
end:





#MakeRama1(L,m,n,t,K): inputs a list L of homog. quadratic polynomials in m and n
#finds the best entry i if it exists together with the constant. The output is [i,c,PairOfRationalFunctions of L[i]=c. Try:
#MakeRama1([m^2-7*m*n-9*n^2, 2*m^2+4*m*n+12*n^2, -2*m^2-10*n^2, -m^2-9*m*n+n^2],m,n,t,100);
MakeRama1:=proc(L,m,n,t,K) local i,gu:

for i from 1 to nops(L) do

 gu:=SolQuad(L[i],m,n,t,K):
 
  if gu<>FAIL then
   RETURN([i,gu[3],gu[1],gu[2]]):
 fi:

od:

FAIL:

end:




#MakeRama(L,m,n,t,K): inputs a list L of homog. quadratic polynomials in m and n, a variable t and a large pos. integer K
#outputs what MakeRama1(L,m,n,t,K) (q.v.)  let's call it [i,c,Pair]
#together with a tuple of the same length as L where
#the i-th entry is replaced by c and the other ones, way L[j], is the generating functions of the expressions for L[j]
#if m and n are replaced by Pair[1],Pair[2]. Try:
#MakeRama([m^2-7*m*n-9*n^2, 2*m^2+4*m*n+12*n^2, -2*m^2-10*n^2, -m^2-9*m*n+n^2],m,n,t,100);
MakeRama:=proc(L,m,n,t,K) local  gu,i,c,Pair,L1,mu,i1:
gu:=MakeRama1(L,m,n,t,K):
if gu=FAIL then
 RETURN(FAIL):
fi:

i:=gu[1]:

c:=gu[2]:

Pair:=[gu[3],gu[4]]:

L1:=[seq(L[i1]/c,i1=1..i-1),seq(L[i1]/c,i1=i+1..nops(L))]:
mu:=SR(Pair,t,L1,m,n,K):

mu:=[op(1..i-1,mu),1,op(i..nops(mu),mu)]:

c:=lcm(seq(coeff(denom(mu[i]),t,0),i=1..nops(mu))):

[gu,c*mu]:

end:







#MakeRamaSol1(L,m,n,t,K): inputs a list L of homog. quadratic polynomials in m and n
#finds the best entry i if it exists together with the constant. The output is [i,c,PairOfRationalFunctions of L[i]=c. Try:
#MakeRamaSol1([m^2-7*m*n-9*n^2, 2*m^2+4*m*n+12*n^2, -2*m^2-10*n^2, -m^2-9*m*n+n^2],m,n,t,100);
MakeRamaSol1:=proc(L,m,n,t,K) local i,gu:

for i from 1 to nops(L) do
 gu:=SolDE2(L[i]-1,m,n,t,K):

  if gu<>{} then
    gu:=gu[1]:
   RETURN([i,1,gu[1],gu[2]]):
 fi:

od:

FAIL:

end:




#MakeRamaSol(L,m,n,t,K): inputs a list L of homog. quadratic polynomials in m and n, a variable t and a large pos. integer K
#outputs what MakeRama1(L,m,n,t,K) (q.v.)  let's call it [i,c,Pair]
#together with a tuple of the same length as L where
#the i-th entry is replaced by c and the other ones, way L[j], is the generating functions of the expressions for L[j]
#if m and n are replaced by Pair[1],Pair[2]. Try:
#MakeRamaSol([m^2-7*m*n-9*n^2, 2*m^2+4*m*n+12*n^2, -2*m^2-10*n^2, -m^2-9*m*n+n^2],m,n,t,100);
MakeRamaSol:=proc(L,m,n,t,K) local  gu,i,c,Pair,L1,mu,i1:
gu:=MakeRamaSol1(L,m,n,t,K):
if gu=FAIL then
 RETURN(FAIL):
fi:

i:=gu[1]:

c:=gu[2]:

Pair:=[gu[3],gu[4]]:

L1:=[seq(L[i1]/c,i1=1..i-1),seq(L[i1]/c,i1=i+1..nops(L))]:
mu:=SR(Pair,t,L1,m,n,K):

mu:=[op(1..i-1,mu),1,op(i..nops(mu),mu)]:

c:=lcm(seq(coeff(denom(mu[i]),t,0),i=1..nops(mu))):

[gu,c*mu]:

end:


#TheoremAB(A,B,m,n,t,X,Y,Z,K1,K2): Inputs positive integers A and B, letters m,n, t and X and guessing parameters K1 and K2
#(fairly large positive integers, start K1 around 10 and K2 around 100), outputs
#(i) A list of length 4, let's call it L, consisting of parametric equation for the solution of A*X[1]^3+A*X[2]^3+B*X[3]^3+B*X[4]^3=0 in terms of
# quadratic polynomials in m, n
#(ii) A list of length 4 of the form [Place,constant,f1,f2] where 1<=Place<=4, constant is a rational number and f1 and f2 are rational functions of
# the variable t such that if a(i) is the coefficient of t^i  in the power-series expansion of f1
#and b(i) is the coefficient of t^i  in the power-series expansion of f2 then L[Place](a(i),b(i))=consant for all i. In other wors
#it is the solution of the Diphantine equation  L[Place](m,n)=constant in the C-finite ansatz
#(iii) A list of length 4 there of whose entries are rational functions in t and one of them is a constant such that
#they form solutions of A*X^3+B*(Y^3+Z^3)=CONSTANT
#(iv) A Diophantine equation of the from C1*X[1]^3+C2*X[2]^3+C2*X[3]^4=CONSTANT where two of the {C1,C2,C3} are equal (they are either A or B)
#(v): A triple of rational functions in t whose MacLaurin coefficients form solutions of that Diophantine equation. Try:
#TheoremAB(1,2,m,n,t,X,K1,K2);
TheoremAB:=proc(A,B,m,n,t,X,Y,Z,K1,K2) local gu,lu,ku,Eq,f1,f2,f3,i:
option remember:

if not (type(A,integer) and type(B,integer) and type(m,symbol) and type(n,symbol) and type(X,symbol) and type(Y,symbol) and type(Z,symbol) 
and type(K1,integer) and type(K2,integer)) then
 print(`Bad input`):
 RETURN(type):
fi:

ku:=[]:

gu:=ParAB3(A,B,m,n,K1):
if gu=FAIL then
 RETURN(FAIL):
fi:

ku:=[gu]:

lu:=MakeRama(gu[2],m,n,t,K2):

if lu=FAIL then
 RETURN(ku):
fi:

ku:=[op(ku),op(lu)]:

if ku[2][1]<>1 then
# print(`Not yet implemented`):
 RETURN(ku):
fi:

if normal(diff(ku[3][1],t))<>0 then
print(`Somethins is wrong`):
 RETURN(FAIL):
fi:

Eq:=A*X^3+B*Y^3+B*Z^3=-A*lu[2][1]^3:

ku:=[gu,[op(2..4,lu[1])],Eq,[op(2..4,lu[2])]]:

f1:=taylor(lu[2][2],t=0,101):
f2:=taylor(lu[2][3],t=0,101):
f3:=taylor(lu[2][4],t=0,101):

if {seq(evalb(subs({X=coeff(f1,t,i),Y=coeff(f2,t,i),Z=coeff(f3,t,i)},Eq)),i=1..100)}<>{true} then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

ku:
end:



#TheoremABs(A,B,t,X,Y,Z,K1,K2): States a theorem about infinitely many  solutions  to the
#diophantine equation
#A*X^3+B*Y^3+B*Z^3=Constant, for some constant.
#Try:
#TheoremABs(1,2,t,X,Y,Z,10,100):
TheoremABs:=proc(A,B,t,X,Y,Z,K1,K2) local m,n,gu,a,b,c,i:
gu:=TheoremAB(A,B,m,n,t,X,Y,Z,K1,K2):

if gu=FAIL or nops(gu)=1 then
 RETURN(FAIL):
fi:

print(`Infinitely many solutions to the Diophantine Equation`, gu[3]):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let a(i), b(i), c(i) be defined in terms of the following generating functions`):
print(``):
print(Sum(a(i)*t^i,i=0..infinity)=gu[4][1]):
print(``):
print(Sum(b(i)*t^i,i=0..infinity)=gu[4][2]):
print(``):
print(Sum(c(i)*t^i,i=0..infinity)=gu[4][3]):
print(``):

print(`In Maple notation, these generating functions are`):

lprint(gu[4][1]):
print(``):
lprint(gu[4][2]):
print(``):
lprint(gu[4][3]):
print(``):

print(`Then for all i>=0 we have`):
print(``):
print(subs({X=a(i),Y=b(i),Z=c(i)},gu[3])):
print(``):
print(`Proof: Since everything is in the C-finite ansatz, checking it for i from 0 to 10 suffices, which we did and you are welcome to check.`):
print(``):
end:




#TheoremABs1(A,B,t,X,Y,Z,K1,K2,co): States  theorem number co  about infinitely many  solutions  to the
#diophantine equation
#A*X^3+B*Y^3+C*Z^3=Constant, for some constant.
#Try:
#TheoremABs1(1,2,t,X,Y,Z,20,200,1):
#TheoremABs1(1,2,m,n,t,X,10,100,1);
TheoremABs1:=proc(A,B,t,X,Y,Z,K1,K2,co) local m,n,gu,a,b,c,i:
gu:=TheoremAB(A,B,m,n,t,X,Y,Z,K1,K2):
if gu=FAIL or nops(gu)=1 then
 RETURN(FAIL):
fi:
print(``):
print(`------------------------------------------`):
print(``):
print(`Theorem Number `, co, `: Let a(i), b(i), c(i) be defined in terms of the following generating functions`):
print(``):
print(Sum(a(i)*t^i,i=0..infinity)=gu[4][1]):
print(``):
print(Sum(b(i)*t^i,i=0..infinity)=gu[4][2]):
print(``):
print(Sum(c(i)*t^i,i=0..infinity)=gu[4][3]):
print(``):

print(`In Maple notation, these generating functions are`):

lprint(gu[4][1]):
print(``):
lprint(gu[4][2]):
print(``):
lprint(gu[4][3]):
print(``):

print(`Then for all i>=0 we have`):
print(``):
print(subs({X=a(i),Y=b(i),Z=c(i)},gu[3])):
print(``):
print(`Proof: Since everything is in the C-finite ansatz, checking it for i from 0 to 10 suffices, which we did and you are welcome to check.`):
print(``):
end:





#Paper(C,t,X,Y,Z,K1,K2): Inputs a positive integer B, variables t, X,Y,Z, and positive integers K1 and K2 (guessing parameters)
#outputs an article about Diophantine equations of the form
#A*X^3+B*Y^3+C*Z^3=Constant, for some constant for A<B<=C where A and B are relatively prime.
#Try:
#Paper(3,t,X,Y,Z,20,200):
Paper:=proc(C,t,X,Y,Z,K1,K2) local m,n,A,B,co,gu,t0,CONSTANT:
t0:=time():
print(``):
print(`Inifintely many solutions to Diophantine Equations of the form`, A*X^3 + B*Y^3+ B*Z^3=CONSTANT,  `for A and B At Most`, C):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

co:=0:

for A from 1 to C do
 for B from A+1 to C do
  if igcd(A,B)=1 then
   gu:=TheoremAB(A,B,m,n,t,X,Y,Z,K1,K2):
    if (gu<>FAIL and nops(gu)=4) then
     co:=co+1:
     TheoremABs1(A,B,t,X,Y,Z,K1,K2,co):
    fi:
  fi:
 od:
od:

print(`-----------------------------------------------------`):
print(``):
print(`This ends this fascinating article that stated`, co , `theorems and   took`, time()-t0, `seconds to generate`):
print(``):
end:







#TheoremABsV(A,B,m,n,t,X,Y,Z,K1,K2): States and EXPLAINS a  theorem about infinitely many  solutions  to the
#diophantine equation
#A*X^3+B*Y^3+B*Z^3=Constant, for some constant.
#Try:
#TheoremABsV(1,2,m,n,t,X,Y,Z,10,100):
TheoremABsV:=proc(A,B,m,n,t,X,Y,Z,K1,K2) local gu,a,b,c,i,x,y,W,ka:
gu:=TheoremAB(A,B,m,n,t,X,Y,Z,K1,K2):

if gu=FAIL or nops(gu)=1 then
 RETURN(FAIL):
fi:

print(`Infinitely many solutions to the Diophantine Equation`, gu[3]):
print(`And Explaining How They Were Discovered`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`We are interested in finding infinitely many solutions of the cubic diophantine equation`):
print(``):
print(gu[3]):
print(``):
print(`We first consider the diophantine equation in FOUR variables, W,X,Y,Z`):
print(``):
print(A*W^3+A*X^3+B*Y^3+B*Z^3=0):
print(``):
prin(`We first seach, by brute force solutions where W,X,Y,Z are less than`, K1):
print(``):
print(`and get the following special case`):
print(``):
print(W=gu[1][1][1], X=gu[1][1][2], Y=gu[1][1][3], Z=gu[1][1][4]):
print(``):
print(`Obviously, the following doubly-infinite family is also a solution`):
print(``):
print(W=m, X=-m, Y=n, Z=-n):
print(``):
print(`For any integers m and n`):
print(``):
print(`Using the Eri Jabotinsky trick writing`):
print(``):
print(W=gu[1][1][1]*x+m*y, X=gu[1][1][2]*x-m*y,  Y=gu[1][1][3]*x+n*y, Z=gu[1][1][4]*x-n*y):
print(``):
print(`Plugging into the euqation`):
print(``):
print(A*W^3+A*X^3+B*Y^3+B*Z^3=0):
print(``):
ka:=factor(A*(gu[1][1][1]*x+m*y)^3+A*(gu[1][1][2]*x-m*y)^3+B*(gu[1][1][3]*x+n*y)^3+B*(gu[1][1][4]*x-n*y)^3):
print(`we get `):
print(ka=0):
print(``):
print(`dividing by `, x*y, ` we have `):
print(``):
ka:=normal(ka/(x*y)):
ka:=factor(coeff(coeff(ka,x,1),y,0))*x+factor(coeff(coeff(ka,y,1),x,0))*y:
print(``):
print(ka=0):
print(``):
print(`Obviously the following values of x and y will make it zero`):
print(``):
print(x=-coeff(ka,y,1)	, y=coeff(ka,x,1)):
print(``):
print(`Plugging them back, we get the following parametric solution to `, A*W^3+A*X^3+B*Y^3+B*Z^3=0):
print(``):
print(W=gu[1][2][1], X=gu[1][2][2], Y=gu[1][2][3], Z=gu[1][2][4]):
print(``):

if expand(A*gu[1][2][1]^3+A*gu[1][2][2]^3+B*gu[1][2][3]^3+B*gu[1][2][4]^3)<>0 then
  print(`Something went wrong `):
  RETURN(FAIL):
fi:
print(``):
print(`Check!`):
print(``):
print(`We now need a Lemma `):
print(``):
print(`Lemma: Let d(i), e(i) be defined in terms of the following generating functions `):
print(``):
print(Sum(d(i)*t^i,i=0..infinity)=	gu[2][2]):
print(``):
print(Sum(e(i)*t^i,i=0..infinity)=	gu[2][3]):
print(``):
print(`Then m=d(i), n=e(i) satisfy the quadratic diophantine equation`):
print(``):
print(gu[1][2][1]=gu[2][1]):
print(``):
print(`The proof is routine, since everything is in the C-finite ansatz, and in fact checing it for the first three cases suffices,`):
print(``):
print(`but it was discovered using the algorithm for solving qudaratic diophanine equations`):
print(``):
print(`Now that we have C-finite representations for m and n, we can plug them into`):
print(``):
print( X=gu[1][2][2], Y=gu[1][2][3], Z=gu[1][2][4]):
print(``):
print(`and get C-finite representations for X,Y,Z, leading to the following theorem `):
print(``):

print(`Theorem: Let a(i), b(i), c(i) be defined in terms of the following generating functions `):
print(``):
print(Sum(a(i)*t^i,i=0..infinity)=gu[4][1]):
print(``):
print(Sum(b(i)*t^i,i=0..infinity)=gu[4][2]):
print(``):
print(Sum(c(i)*t^i,i=0..infinity)=gu[4][3]):
print(``):

print(`In Maple notation, these generating functions are`):

lprint(gu[4][1]):
print(``):
lprint(gu[4][2]):
print(``):
lprint(gu[4][3]):
print(``):

print(`Then for all i>=0 we have`):
print(``):
print(subs({X=a(i),Y=b(i),Z=c(i)},gu[3])):
print(``):
print(`Proof: Since everything is in the C-finite ansatz, checking it for i from 0 to 10 suffices, which we did and you are welcome to check.`):
print(``):
print(`--------------------------------`):
print(``):
end:






#PaperV(C,t,X,Y,Z,K1,K2): Inputs a positive integer B, variables t, X,Y,Z, and positive integers K1 and K2 (guessing parameters)
#outputs a MOTIVATED article about Diophantine equations of the form
#A*X^3+B*Y^3+C*Z^3=Constant, for some constant for A<B<=C where A and B are relatively prime.
#Try:
#PaperV(3,t,X,Y,Z,20,200):
PaperV:=proc(C,m,n,t,X,Y,Z,K1,K2) local A,B,co,gu,t0,CONSTANT:
t0:=time():
print(``):
print(`Finding Inifintely many solutions to Diophantine Equations of the form`, A*X^3 + B*Y^3+ B*Z^3=CONSTANT,  `for A and B At Most`, C):
print(`With Detailed Explanation and Motovation`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

co:=0:

for A from 1 to C do
 for B from A+1 to C do
  if igcd(A,B)=1 then
   gu:=TheoremAB(A,B,m,n,t,X,Y,Z,K1,K2):
    if (gu<>FAIL and nops(gu)=4) then
     co:=co+1:
 print(`Chapter `, co):


     TheoremABsV(A,B,m,n,t,X,Y,Z,K1,K2,co):

    fi:
  fi:
 od:
od:

print(`-----------------------------------------------------`):
print(``):
print(`This ends this fascinating article that contained `, co , `chapters and   took`, time()-t0, `seconds to generate`):
print(``):
end: