#####################################################################
## Pell.txt Save this file as  Pell.txt  to use it,                  #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read Pell.txt <Enter>                             #
## 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 Pell.txt, it is one of the Maple packages `):
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/Pell.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(`For a list of the supporting functions type: ezra1();`):
print():
 



ezra1:=proc()
if args=NULL then

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

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

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

print(`The MAIN procedures are:`):
print(` AllPellGF, AllPellGFv, DErat, PellFt, PellFtV, PellGF, PellLike, PellPaper, 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,400);`):

elif nargs=1 and args[1]=AllPellGFv then
print(`AllPellGFv(N,x,K,X,Y): Verbose form of AllPellGF(N,x,K) for Pell's Equation X^2-n*Y^2=1 for all non-square n from 2 to N^2-1 . Try:`):
print(`AllPellGFv(10,t,400,X,Y);`):

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]=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]=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]=PellFtV then
print(`PellFtV(N1,N2,M,t,X,Y): Like PellFt(N1,N2,M,t,X,Y) but in the format [N1*X^2-N2*Y^2=CONST,f1(t),f2(t)] where`):
print(`f1(t), f2(t) are the generating functions of the solutions (X,Y)=(a(i),b(i)) to N1*X^2-N2*Y^2=CONST. Try:`):
print(`PellFtV(1,19,100,t,X,Y);`):


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]=PellPaper then
print(`PellPaper(M,K,t,X,Y):  outputs generating function solutions to Pell-Like equations`):
print(`N1*X^2-N2*Y^2=CONSTANT for various constants in terms of t for all N2 from 2 to M that is not a perfect square`):
print(`and N1 from 1 to N2-1 not relatively prime to N2 with guessing parameter K. Try:`):
print(`PellPaper(10,500,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]=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]=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:



#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:

#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:



##############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:




#PellFtV(N1,N2,M,t,X,Y): Like PellFt(N1,N2,M,t,X,Y) but in the format [N1*X^2-N2*Y^2=CONST,f1(t),f2(t)] where
#f1(t), f2(t) are the generating functions of the solutions (X,Y)=(a(i),b(i)) to N1*X^2-N2*Y^2=CONST. Try:
#PellFtV(1,19,100,t,X,Y);
PellFtV:=proc(N1,N2,M,t,X,Y) local  gu,i:
gu:=PellFt(N1,N2,M,t):

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

[seq([N1*X^2-N2*Y^2=gu[i][1],gu[i][2],gu[i][3]],i=1..nops(gu))]:

end:


#PellPaper(M,K,t,X,Y):  outputs generating function solutions to Pell-Like equations
#N1*X^2-N2*Y^2=CONSTANT for various constants in terms of t for all N2 from 2 to M that is not a perfect square
#and N1 from 1 to N2-1 not relatively prime to N2 with guessing parameter K. Try:
#PellPaper(10,500,t,X,Y);
PellPaper:=proc(M,K,t,X,Y) local N1,N2,CONSTANT,gu,j,t0:

t0:=time():

print(`Generating Functions For solutions of Pell-Like Equations of the from `, N1*X^2-N2*Y^2=CONSTANT, ` for various small CONSTANT`):
print(`for Non-Square N2 less than`, M, `and N1 less than N2 `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`In the table below, given in computer-format, so that it can be transported, the convention is `):
print(``):
print([P(X,Y), f1(t),f2(t)]):
print(``):
print(`meaning that if a(i) and b(i) are the coefficients of t^i in the power-series expansion of f1(t), f2(t) respetively`):
print(`then P(a(i),b(i)) is correct`):

print(``):

for N2 from 1 to M do
 if not type(sqrt(N2),integer) then
print(``):
 print(`doing N2=`, N2):
print(``):
 for N1 from 1 to N2-1 do
   if igcd(N1,N2)=1 then
    gu:=PellFtV(N1,N2,K,t,X,Y):
     if gu<>FAIL then
       for j from 1 to nops(gu) do
        lprint(gu[j]):
        print(``):
       od:
      fi:
   fi:
 od:
 fi:
od:

print(``):
print(`This ends this book that took`, time()-t0, `seconds to generate. `):

    
end:



#AllPellGFv(N,x,K,X,Y): Verbose form of AllPellGF(N,x,K) for Pell's Equation X^2-n*Y^2=1 for all non-square n from 2 to N^2-1 . Try:
#AllPellGFv(10,t,400,X,Y);
AllPellGFv:=proc(N,x,K,X,Y) local gu,i, t0,n,k,b:

t0:=time():

gu:=AllPellGF(N,x,K):

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

print(`Generating Functions for soltions of Pell's Equation `, X^2-n*Y^2=1 , `for all non-perfect-square n from 2 to `, N^2-1):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):


print(`In the table below, given in computer-format, so that it can be transported, the convention is `):
print(``):
print([P(X,Y), [k,b], f1(t),f2(t)]):
print(``):
print(`meaning that if a(i) and b(i) are the coefficients of t^i in the power-series expansion of f1(t), f2(t) respetively`):
print(`then P(a(i),b(i)) is correct. [k,b] is the smallest non-trivial solution (of course [1,0] is also a solution)`):

for i from 1 to nops(gu) do
print(``):
k:=-coeff(numer(gu[i][2][1]),t,1):
b:=coeff(numer(gu[i][2][2]),t,1):
lprint(X^2-gu[i][1]*Y^2=1,[k,b],gu[i][2][1],gu[i][2][2]):
print(``):
od:

print(``):
print(`------------------------------------------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `seconds to generate. `):
print(``):
end: