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


print(` Written: March 2020 `):
print():
print(`This is GCF.txt, A Maple package for conjecturing AND proving explicit values of some infinite generalized continued fractions.`):
print(` It accompanies    Robert Dougherty-Bliss's and Doron Zeilberger's  article `):
print(`Automatic Conjecturing and Proving Of Exact Values of Some Infinite Families of Infinite Continued Fractions`):
print(`avaliable from the authors' websites and soon from arxiv.org `):
print(``):
print(`This project was inspired by the innovative article "The Ramanujan machine: Automatically Generated Conjectures on Fundeamental Conjectures"`):
print(``):
print(`By  Gal Raayoni, George Pisha, Yahel Manor, Uri Mendlovic, Doron Haviv, Yaron Hadas, and Ido Kaminer`):
print(``):
print(`arXiv: 1907.00205v3 `):
print(` https://arxiv.org/abs/1907.00205`):


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


print():
print(`---------------------------------------------`):
print():
print(`For a list of the SUPPORTING procedures`):
print(` type "ezra1();". For specific help type "ezra(procedure_name);" `):
print():
print():
print(`---------------------------------------------`):
print():
 

print():
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():
print():
print(`---------------------------------------------`):
print():
 


ezra1:=proc()
if args=NULL then
 print(`The supporting procedures are Doron `):
 print(`Gk `):
else
 ezra(args):
fi:

end:

ezra:=proc()
if args=NULL then

print(` GCF.txt: A Maple package for  conjecturing and proving explicit expressions for some infinite generalized continued fractions.`):
print(``):
print(`The MAIN procedures are: Eq2, GCF, GCFab, RDB , RDBsym, Yaron, YaronV, Euler, Eulersym`):
print(``):


elif nargs=1 and args[1]=Eq2 then
print(`Eq2(k,a): The right side of Eq(2) in the paper. Try:`):
print(`Eq2(4,2);`):

elif nargs=1 and args[1]=Doron then
print(`Doron(n,a,b): Does RDBsym(a*n^2+b*n+1, -a*n^2-b*n,n); a and b are symbolic or numeric, but with enhanced limit taking. Try:`):
print(`Doron(n,a,b);`):

elif nargs=1 and args[1]=GCF then
print(`GCF(L): Inputs a finite list of pairs [[a1,b1],[a2,b2], ..., [ak,bk]] and outputs the value of the generalized continued fraction`):
print(`a1+b1/(a2+b2/(a3+..)... . Try:`):
print(`GCF([[1,1]$10]);`):
print(`GCF([ seq([4*i-2,1],i=1..20) ]);`):

elif nargs=1 and args[1]=GCFab then
print(`GCFab(a,b,n,K): The generalized continued fraction K(a(n),b(n)) up to K places.`):
print(`Inputs expressions a and b in the symbol n, and a positive integer K. Outputs the truncated continued fractions`):
print(`a(1)+b(1)/(a(2)+ b(2)/.... `):
print(`Try: `):
print(`evalf(GCFab(1,1,n,50)-(1+sqrt(5))/2 );`):
print(`evalf(GCFab(4*n-2,1,n,40)-(exp(1)+1)/(exp(1)-1));`):
print(`evalf(GCFab(n+2,-n,n,50)-exp(1.));`):
print(`evalf(GCFab(2*n-1,n^2,n,50)-4/Pi);`):
print(`evalf(GCFab((n-1)^3+n^3,-n^6,n,100)-1/Zeta(3));`):



elif nargs=1 and args[1]=Gk then
print(`Gk(n,a,G,k)  Inputs a symbol n, an integer (or symbol) a, and a symbol G, and a non-negative integer k, outputs G[n+k] in terms `):
print(`G[n] where G[n] stands for the incomplete Gamma function GAMMA(n+2,a):=int(exp(-x)*x^(n+1),x=a..infinity).`):
print(`Try `):
print(`Gk(n,a,G,2); `):

elif nargs=1 and args[1]=RDB then
print(`RDB(a,b,n): Inputs expressions a and b in n and outputs (if successful) the explicit expressions for the numerator and denominator`):
print(`of the n-th convergent of the infinite continued fraction, as well as its limit, and the error by using`):
print(`the 50-th. If it fails, it returns FAIL. Try:`):
print(`RDB(n+3,-n,n);`):
print(`RDB(n+11,-3*n,n);`):
print(`To prove the conjecture mentioned in our paper on p. 8. Type:`):
print(`RDB(4*n^2+6*n+1,-4*n^2-6*n,n);`):
print(``):
print(`To prove the impressive  "The Ramanujan Machine" infinite continuted fraction for 1/Zeta(3), type: `):
print(` RDB((n-1)^3+n^3,-n^6,n);`):

elif nargs=1 and args[1]=RDBsym then
print(`RDBsym(a,b,n): Like RDB(a,b,n) but now the expressions a and b in terms of n may depend on symbolic parameter`):
print(`so there is no  convergence-checking. Try:`):
print(`RDBsym(2*n-2,1,n); `):
print(`RDBsym(n+3,-n,n);`):
print(`RDBsym(c*n+1,-c*n,n);`):


elif nargs=1 and args[1]=Yaron then
print(`Yaron(k1,a1,n,G): Inputs a pos. integer k1 and a non-zero integer a1, and a symbol n, and a symbol G`):
print(`outputs the output of RDB(n+k1,a1*n,n) where GAMMA(n+2,-a1) is denoted by G[n]. It also proves it rigorously. Try:`):
print(`Yaron(4,1,n,G);`):


elif nargs=1 and args[1]=YaronV then
print(`YaronV(k1,a1,n,G): verbose form of Yaron(k1,a1,n,G) (q.v. ). Try:`):
print(`YaronV(4,1,n,G);`):

elif nargs=1 and args[1]=Euler then
print(`Euler(r, n): Derives Euler's continued fraction for the sequence r(n). Try:`):
print(`Euler(1 / n, n)`);
print(`Euler(1, n)`);


elif nargs=1 and args[1]=Eulersym then
print(`Eulersym(r, n): Derives Euler's continued fraction for the sequence r(n) without numeric checking. Try:`):
print(`Euler(k / n, n)`);

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

fi:


end:


Digits:=100:

#GCF(L): Inputs a finite list of pairs [[a1,b1],[a2,b2], ..., [ak,bk]] and outputs the value of the generalized continued fraction
#a1+b1/(a2+b2/(a3+..)... . Try:
#GCF([[1,1]$10]);
#GCF([ seq([1,4*i-2],i=1..20) ]);

GCF:=proc(L) local gu:
if nops(L)=1 then
 RETURN(L[1][1]):
fi:

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

normal(L[1][1]+L[1][2]/gu):

end:

#GCFab(a,b,n,K): The generalized continued fraction K(a(n),b(n)) up to K
#Try
#GCFab(1,1,n,10);
#GCFab(2*n-2,1,n,10);
#GCFab(n+3,-n,n,10);
GCFab:=proc(a,b,n,K) local L,n1:
L:=[seq([eval(a, n=n1),eval(b, n=n1)],n1=1..K)]:
GCF(L):
end:


#p_soln := rsolve({p(n + 2) - (n + 6) * p(n + 1) + (n + 2) * p(n), p(0) = 4, p(1) = 19}, p(n)):
#q_soln := rsolve({q(n + 2) - (n + 6) * q(n + 1) + (n + 2) * q(n), q(0) = 1, q(1) = 5}, q(n)):

#RDB(a,b,n): Inputs expressions a and b in n and outputs the explicit expressions for the numerator and denominator
#of the n-th convergent of the infinite continued fraction, as well as its limit, and the error by using
#the 50-th. Try:
#RDB(2*n-2,1,n); 
#RDB(n+3,-n,n);
RDB:=proc(a,b,n) local p,ps,qs,c,er:

ps := simplify(rsolve({p(n + 2) - eval(a, n=n+3)*p(n + 1) - eval(b, n=n + 2) * p(n), 
p(0) = eval(a, n=1), p(1) =eval(a, n = 1)*eval(a, n=2)+eval(b, n=1)}, p(n))):

if type(ps, function) then
 RETURN(FAIL):
fi:

qs := simplify(rsolve({p(n + 2) - subs(n=n+3,a)*p(n + 1) - subs(n=n + 2,b) * p(n), p(0) = 1, p(1) =subs(n=2,a)}, p(n))):


c:= limit(simplify(ps/ qs), n=infinity):

er:=evalf(GCFab(a,b,n,300)-c):
if abs(er)>1/10^5 then
 print(` WARNING: the convergence may be very slow. OR IT IS WRONG` ):
fi:

[ps, qs,c,er];


end:



#Gk(n,a,G,k)  Inputs a symbol n, an integer a, a symbol G, and a non-negative integer l, outputs G[n+k] in terms 
#G[n] where G[n] stands for the incomplete Gamma function GAMMA(n+2,a)=int(exp(-x)*x^(n+1),x=a..infinity).
#Try
#Gk(n,a,G,2);
Gk:=proc(n,a,G,k)  local gu:
option remember:
if k=0 then
 RETURN(G[n]):
fi:

if k=1 then 
RETURN(a^(n+2)*exp(-a)+(n+2)*G[n]):
fi:


gu:=Gk(n,a,G,k-1):
gu:=subs(n=n+1,gu):
gu:=collect(expand(subs(G[n+1]=Gk(n,a,G,1),gu)),G[n]):

RETURN(gu):


end:

#Yaron(k1,a1,n,G): Inputs a pos. integer k1 and a non-zero integer a1, and a symbol n, and a symbol G
#outputs the output of RDB(n+k1,a1*n,n) where GAMMA(n+2,-a1) is denoted by G[n]. Try:
#Yaron(4,1,n,G);
Yaron:=proc(k1,a1,n,G) local gu,f1,f2,a,b,hal1,hal2:
a:=n+k1:
b:=a1*n:
gu:=RDB(n+k1,a1*n,n):


gu:=convert(gu,factorial):
gu:=subs(GAMMA(n+2,-a1)=G[n],gu):
f1:=gu[1]:
f2:=gu[2]:
hal1:=expand(subs(n=n+2,f1) - subs(n=n+3,a)*subs(n=n+1,f1) - subs(n=n + 2,b) * f1):
hal1:=simplify(expand(subs({G[n+2]=Gk(n,-a1,G,2),G[n+1]=Gk(n,-a1,G,1)},hal1))):
hal2:=expand(subs(n=n+2,f2) - subs(n=n+3,a)*subs(n=n+1,f2) - subs(n=n + 2,b) * f2):
hal2:=simplify(expand(subs({G[n+2]=Gk(n,-a1,G,2),G[n+1]=Gk(n,-a1,G,1)},hal2))):

if [hal1,hal2]<>[0,0] then
 print(`Something terrible happened`):
 print(gu):
 print(`did not work out`):
 RETURN(FAIL):
fi:

gu:

end:


#YaronV(k1,a1,n,G): a verbose version of Yaron(ka,a1,nG). Outputs an article. Try
##YaronV(4,1,n,G);
YaronV:=proc(k1,a1,n,G) local gu,a,b,x,lu1,lu2,P,Q,t0:

t0:=time():
lu1:=GCFab(n+k1,a1*n,n,1):
lu2:=GCFab(n+k1,a1*n,n,2):

gu:=Yaron(k1,a1,n,G):

if gu=FAIL then
 RETURN(FAIL):
fi:
print(``):
print(`The exact value of the infinite generalized continued fraction a(1)+b(1)/(a(2)+b(2)/... where `, a(n)=n+k1, `and `, b(n)= a1*n):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Consider the following Infinite continued fraction, let's call it C:`):
print(``):
print(`C=a(1)+ b(1)/(a(2)+ b(2)/(a(3)+b(3)/(a(4)+.. `):
print(`where `):
print(a(n)=n+k1, b(n)=a1*n ):
print(``):
print(`Then the exact value of C equals`):
print(``):
print(gu[3]):
print(``):
print(`Its floating-point approximation is`):
print(``):
print(evalf(gu[3])):
print(``):
print(`Proof: `):
print(``):
print(`Let's denote  the incomplete Gamma function`):
print(``):
print( G[n]= GAMMA(n+2,-a1) ):
print(``):
print(`in other words `):
print(``):
print(G[n]=Int(x^(n+1)*exp(-x),x=-a1..infinity)):
print(``):
print(`By integration by parts it is very easy to deduce that G[n] satisfies the first-order (inhomog.) linear recurrence `):
print(``):
print(G[n+1]=Gk(n,a1,G,1)):
print(``):
print(`and hence, iterating, we also have `):
print(``):
print(G[n+2]=Gk(n,a1,G,2)):
print(``):
print(`Let P(n) and Q(n) be the numerator and denominator of the truncated version of the above infinite continued fraction`):
print(``):
print(`truncated at a(n),b(n) `):
print(``):
print(`Then we have the following explicit expressions for P(n) and Q(n) `):
print(``):
print(P(n)=gu[1]):
print(``):
print(Q(n)=gu[2]):
print(``):
print(`These follow from the fact that by the general theory of continued fractions (and is also easily seen) P(n) and Q(n) both satisfy`):
print(``):
print(`the following second-order linear recurrence `):
print(``):
print(P(n+2) - subs(n=n+3,n+k1)*P(n+1) - subs(n=n + 2,a1*n)*P(n)=0):
print(``):
print(``):
print(Q(n+2) - subs(n=n+3,n+k1)*Q(n+1) - subs(n=n + 2,a1*n)*Q(n)=0):
print(``):
print(`subject to the initial conditions `):
print(``):
print(P(0)=numer(lu1), P(1)=numer(lu2)):
print(``):
print(Q(0)=denom(lu1), Q(1)=denom(lu2)):
print(``):
print(`But the above proposed expressions for P(n) and Q(n) in terms of G[n] and n satisfy the correct initial conditions and satisfy the`):
print(`the recurrence, as is readily checked by using the above expressions of G[n+1] and G[n+2] in terms of G[n] and simplifying`):
print(``):
print(`Finally, taking the limit as n goes to infinity of P(n)/Q(n)`):
print(``):
print(`we get that C indeed equals `):
print(``):
print(gu[3]):
print(``):
print(`QED. `):
print(``):
print(`Comment: Just as a sanity check, the difference between C and P(100)/Q(100) equals`, gu[4]):
print(``):
print(`This ends this article, that took`, time()-t0, `seconds to generate. `):
print(``):
end:






#RDBsym(a,b,n): Like RDB(a,b,n) but now the expressions a and b in terms of n may depend on symbolic paratmer
#so there is no error-checking.
#RDBsym(2*n-2,1,n); 
#RDBsym(n+3,-n,n);
#RDBsym(c*n+1,-c*n,n);
RDBsym:=proc(a,b,n) local lu1,lu2,p,ps,qs,c:
lu1:=GCFab(a,b,n,1):
lu2:=GCFab(a,b,n,2):

ps := simplify(rsolve({p(n + 2) - subs(n=n+3,a)*p(n + 1) - subs(n=n + 2,b) * p(n), p(0) = numer(lu1), p(1) =numer(lu2)}, p(n))):

if type(ps, function) then
 RETURN(FAIL):
fi:

qs := simplify(rsolve({p(n + 2) - subs(n=n+3,a)*p(n + 1) - subs(n=n + 2,b) * p(n), p(0) = denom(lu1), p(1) =denom(lu2)}, p(n))):




c:= limit(simplify(ps/ qs), n=infinity):

[ps, qs,c];


end:




#Doron(n,a,b): Does RDBsym(a*n^2+b*n+1, -a*n^2-b*n,n); a and b are symbolic or numeric, but with enhanced limit taking. Try:
#Doron(n,a,b);
Doron:=proc(n,a,b) local ka2,lu,lu1,lu2:
#ka1:=hypergeom([1],[3, (3*a+b)/a],1/a):
ka2:=hypergeom([1],[n+3, ((n+3)*a+b)/a],1/a):
lu:=RDBsym(a*n^2+b*n+1, -a*n^2-b*n,n):
lu1:=subs(ka2=1,lu[1]):
lu2:=subs(ka2=1,lu[2]):
[lu[1],lu[2],limit(lu1/lu2, n=infinity)]:

end:


#Eq2(k,a): The right side of Eq(2) in the paper
Eq2:=proc(k,a) local s:
normal(a^(a+k+1)/(a+k)!/(exp(a)-add(a^s/s!,s=0..a+k))):
end:


# Euler(r, n): Derives Euler's continued fraction for the sequence r(n).
# Returns [p, q, limit, approx, error], where

#       `p` and `q` are the numerator and denominator sequences of the
#       convergents after the initial terms;
#
#       `limit` is the full limit of the continued fraction;
#
#       `approx` is the numerical approximation of the 1000th convergent; and
#
#       `error` is the numerical error of the 1000th convergent.

# Try:
# Euler(1 / n, n)
# Euler(1, n)
Euler:=proc(r, n)
    local shift_r, a, b, res, lim, appx, e:
    shift_r := subs(n = n - 1, r):
    a := piecewise(n = 1, 1, n >= 2, shift_r + 1):
    b := -r:
    res := RDB(a, b, n):
    lim := simplify(1 / res[3]):
    appx := evalf(1 / GCFab(a, b, n, 1000));
    e := evalf(lim - appx);

    [res[1], res[2], lim, appx, e];
end:

# EulerGCF(r, n, K): Computes the Kth convergent of Euler's continued fraction.
# Try:
# EulerGCF(1 / n, n, 1000)
# EulerGCF(1, n, 1000)
EulerGCF:=proc(r, n, K)
    local shift_r, a, b:
    shift_r := subs(n = n - 1, r):
    a := piecewise(n = 1, 1, n >= 2, shift_r + 1):
    b := -r:
    1 / GCFab(a, b, n, K);
end:

# Eulersym(r, n): The same as Euler(r, n), but without numeric approximations,
# so symbolic arguments can be used.
# Returns [p, q, limit], where

#       `p` and `q` are the numerator and denominator sequences of the
#       convergents after the initial terms. Specifically, 1 / (1 - p / q)
#
#       `limit` is the full limit of the continued fraction;

# Try:
# Euler(k / n, n)
Eulersym:=proc(r, n)
    local shift_r, a, b, res, lim:
    shift_r := subs(n = n - 1, r):
    a := piecewise(n = 1, 1, n >= 2, shift_r + 1):
    b := -r:
    res := RDBsym(a, b, n):
    if res[3] = 0 then
        lim := infinity:
    else
        lim := 1 / res[3]:
    fi:
    [res[1], res[2], lim];
end: