###################################################################### ## RatDio.txt Save this file as RatDio.txt to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read RatDio.txt # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## DoronZeil at gmail dot com # ###################################################################### print(`First Written: June 2020: tested for Maple 2018 `): print(`Version : June 2020: `): print(): print(`This is RatDio.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/XXX .`): print(`Please report all bugs to: DoronZeil at gmail dot com .`): print(): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): print(``): print(`For a list of the supporting functions type: ezra1();`): print(): ezra1:=proc() if args=NULL then print(`The SUPPORTING procedures are`): print(` DE,DEm, DEgO, FindEqs, FindEqs1 , FindEqsM, FindVars, GenPol, GenPolHst, SolBF, SolBF1, SortLE, SR1 `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The MAIN procedures are: DEg, DEgM, GT, Mamar, RT, RTm, RThm, 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]=GT then print(`GT(d,t,K,c,k): inputs a positive integer d, a variable t, and a pos. integer K `): print(` and symbols c and k (for indexed variables) outputs `): print(` a generic tuple of length d of rational functions whose denominator are the same and have the from `): print(` 1+c1*t+...+c[d-1]*t^(d-1) +t^d (if d is even) or `): print(` 1+c1*t+...+c[d-1]*t^(d-1) -t^d (if d is odd) and numerators are different `): print(` and are polynomials of degree <=d-1 and the coefficients are all from -K to K. Try:`): print(`GT(4,t,10,c,k);`): elif nargs=1 and args[1]=Mamar then print(`Mamar(d,t,K,X,CO): produceds an article with CO theorems using`): print(`RThm1(d,t,K,X) (q.v.). Try:`): print(`Mamar(2,t,20,X,10);`): elif nargs=1 and args[1]=RT then print(`RT(d,t,K): inputs a positive integer d, a variable t, and a pos. integer K outputs`): print(`a random tuple of length d of rational functions whose denominator are the same and have the from`): print(`1+c1*t+...+c[d-1]*t^(d-1) +t^d (if d is even) or`): print(`1+c1*t+...+c[d-1]*t^(d-1) -t^d (if d is odd) and numerators are different`): print(`and are polynomials of degree <=d-1 and the coefficients are all from -K to K. Try:`): print(`RT(4,t,10);`): elif nargs=1 and args[1]=RThm then print(`RThm(d,t,K,X): inputs a positive integer d, a variable t, and a pos. integer K outputs`): print(`a random theorem about a certain diophantine equation having many solutions. Try:`): print(`RThm(4,t,10,X);`): elif nargs=1 and args[1]=RTm then print(`RTm(d,t,K): inputs a positive integer d, a variable t, and a pos. integer K outputs`): print(`a random tuple of length d of rational functions whose denominator are the same and have the from`): print(`1+c1*t+...+c[d-1]*t^(d-1) -t^d (if d is even) or`): print(`1+c1*t+...+c[d-1]*t^(d-1) +t^d (if d is odd) and numerators are different`): print(`and are polynomials of degree <=d-1 and the coefficients are all from -K to K. Try:`): print(`RTm(4,t,10);`): 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 Nexpand(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: #RT(d,t,K): inputs a positive integer d, a variable t, and a pos. integer K outputs #a random tuple of length d of rational functions whose denominator are the same and have the from #1+c1*t+...+c[d-1]*t^(d-1) +t^d (if d is even) or #1+c1*t+...+c[d-1]*t^(d-1) -t^d (if d is odd) and numerators are different #and are polynomials of degree <=d-1 and the coefficients are all from -K to K. Try: #RT(4,t,10); RT:=proc(d,t,K) local ra,Q,P,i,j,gu: if not (type(d,integer) and d>1 and type(K,integer) and K>0 and type(t,symbol) ) then print(`Bad input`): RETURN(FAIL): fi: ra:=rand(-K..K): if d mod 2=0 then Q:=1+add(ra()*t^i,i=1..d-1)+t^d: else Q:=1+add(ra()*t^i,i=1..d-1)-t^d: fi: gu:=[]: for i from 1 to d do P:=add(ra()*t^j,j=0..d-1): gu:=[op(gu),P/Q]: od: gu: end: #RTm(d,t,K): inputs a positive integer d, a variable t, and a pos. integer K outputs #a random tuple of length d of rational functions whose denominator are the same and have the from #1+c1*t+...+c[d-1]*t^(d-1) -t^d (if d is even) or #1+c1*t+...+c[d-1]*t^(d-1) +t^d (if d is odd) and numerators are different #and are polynomials of degree <=d-1 and the coefficients are all from -K to K. Try: #RTm(4,t,10); RTm:=proc(d,t,K) local ra,Q,P,i,j,gu: if not (type(d,integer) and d>1 and type(K,integer) and K>0 and type(t,symbol) ) then print(`Bad input`): RETURN(FAIL): fi: ra:=rand(-K..K): if d mod 2=0 then Q:=1+add(ra()*t^i,i=1..d-1)-t^d: else Q:=1+add(ra()*t^i,i=1..d-1)+t^d: fi: gu:=[]: for i from 1 to d do P:=add(ra()*t^j,j=0..d-1): gu:=[op(gu),P/Q]: od: gu: end: #GT(d,t,K,c,k): inputs a positive integer d, a variable t, and a pos. integer K #and symbols c and k (for indexed variables) outputs #a generic tuple of length d of rational functions whose denominator are the same and have the from #1+c1*t+...+c[d-1]*t^(d-1) +t^d (if d is even) or #1+c1*t+...+c[d-1]*t^(d-1) -t^d (if d is odd) and numerators are different #and are polynomials of degree <=d-1 and the coefficients are all from -K to K. Try: #GT(4,t,10,c,k); GT:=proc(d,t,K,c,k) local Q,P,i,j,gu: if not (type(d,integer) and d>1 and type(K,integer) and K>0 and type(t,symbol) ) then print(`Bad input`): RETURN(FAIL): fi: if d mod 2=0 then Q:=1-add(k[i]*t^i,i=1..d-1)+t^d: else Q:=1-add(k[i]*t^i,i=1..d-1)-t^d: fi: gu:=[]: for i from 1 to d do P:=add(c[i,j]*t^j,j=0..d-1): gu:=[op(gu),P/Q]: od: gu: end: #RThm(d,t,K,X): inputs a positive integer d, a variable t, and a pos. integer K outputs #a random theorem about a certain diophantine equation having many solutions. Try: #RThm(4,t,10,X); RThm:=proc(d,t,K,X) local gu,P,a,i,j,P1: gu:=RT(d,t,K,X): P1:=DEg(gu,t,X,d): print(`Infinitely many solutions to the Diophantine Equation`): print(``): print(P1=0): print(``): print(`By Shalosh B. Ekhad`): print(``): print(`Theorem: `): print(`Let `): print(P(seq(X[i],i=1..d))=P1): print(``): print(`Define `, d, `sequences a[i,j] i goes from 1 to`, d, `and j from 0 to infinity in terms of the generating functions`): print(``): for i from 1 to d do print(``): print(Sum(a[i,j]*t^j,j=0..infinity)=gu[i]): print(``): od: print(``): print(`then for each j from 0 to infinity we have`): print(``): print(P(seq(a[i,j],i=1..d))=0): print(``): print(`-----------------------------`): print(``): end: #RThm1(d,t,K,X,co): inputs a positive integer d, a variable t, and a pos. integer K outputs #a random theorem about a certain diophantine equation having many solutions. Try: #RThm1(4,t,10,X,co); RThm1:=proc(d,t,K,X,co) local gu,P,a,i,j,P1: gu:=RT(d,t,K,X): P1:=DEg(gu,t,X,d): print(`Theorem Number`, co): print(`Let `): print(P(seq(X[i],i=1..d))=P1): print(``): print(`Define `, d, `sequences a[i,j] i goes from 1 to`, d, `and j from 0 to infinity in terms of the generating functions`): print(``): for i from 1 to d do print(``): print(Sum(a[i,j]*t^j,j=0..infinity)=gu[i]): print(``): od: print(``): print(`then for each j from 0 to infinity we have`): print(``): print(P(seq(a[i,j],i=1..d))=0): print(``): print(`-----------------------------`): print(``): end: #Mamar(d,t,K,X,CO): produceds an article with CO theorems using #RThm1(d,t,K,X) (q.v.). Try: #Mamar(2,t,20,X,10); Mamar:=proc(d,t,K,X,CO) local co,t0: t0:=time(): print(``): print(`Solutions of a selection of `, CO, `diophantine equations with`, d, `variables of degree`, d): print(``): print(`By Shalosh B. Ekhad `): print(``): for co from 1 to CO do RThm1(d,t,K,X,co): od: print(``): print(`----------------------------------------`): print(``): print(`This concludes this article that took`, time()-t0, `seconds to generate`): print(``): end: