a gen. base expansion program in Maple

[Simon Plouffe]

##  General base expansion program written by Simon Plouffe, April 1996.
##  This program expand a real number with several algorithms and tries
##  the gfun engine to find an expression for that number.
##  There are 13 algorithms in all and 2 calling procedures.
##  Gfun is a program in Maple that was originally written by F. Bergeron
##  and Simon Plouffe and later greatly improved by Salvy and Zimmermann (INRIA)
##  It is now part of the share package in Maple, to get it
##  type with(share);readshare(gfun,calculus);
##
##  Trydev and Tryhard
##  They are both called with 1 parameter : the number.
##  That number is evaluated in floating point within the procedure, so
##  a call with Trydev(cosh(1)+sin(1)) will work.
##  Both procedures accept numeric and symbolic answers, complex numbers
##  are converted to absolute values.
##  Trydev is using the convert(",ratpoly) of Maple to guess the rational
##  fraction of the elements (in power series) of the expansion.
##  Tryhard will use guessgf(",x) of the GFUN package. (5 minutes with 24 
##  digits precision). 
##  
##  To test a number X , type Trydev(X);  for a quick answer,
##                            Tryhard(X); for a deeper answer (longer)
###############################################################################
#                               Continued Fraction 
#
#                                 y(n) = [1/x(n)]
#                                x(n+1) = {1/x(n)}
#
#
r2cf:=proc(s)
local max,liste,S,liste2,nn;
    S :=  evalf(abs(s));
    max :=  10^(Digits-2);
    liste2:=convert(S,confrac);
    nn:=nops(liste2)-3;
    liste:=[op(1..nn,liste2)];
     RETURN(liste);
end:
######################################################################
#                               Egyptian fraction 
#
#                               y(n) = 1+[1/x(n)]
#                              x(n+1) = x(n)-1/y(n)
#
#
r2ef:=proc(s)
local max,liste,prod,S,T;
    S :=  evalf(frac(abs(s)));
    max :=  10^(Digits-5);
    prod :=  1;
    liste :=  [trunc(s)];
    while prod <=  max do
        T :=  1+trunc(1/S);
        S :=  S-1/T;
        liste :=  [op(liste),T];
        prod :=  prod*T
    od;
     RETURN(liste);
end:
###############################################################################
#                               Egyptian product 
#
#                               y(n) = 1+[1/x(n)]
#                              x(n+1) = {x(n)*y(n)}
#
#
r2ep:=proc(s)
local max,liste,prod,S,T;
    S :=  evalf(frac(abs(s)));
    max :=  10^(Digits-5);
    prod :=  1;
    liste :=  [trunc(s)];
    while prod <=  max do
        T :=  1+trunc(1/S);
        S :=  frac(S*T);
        liste :=  [op(liste),T];
        prod :=  prod*T
    od;
     RETURN(liste);
end:
###############################################################################
#                         Alternated Egyptian fraction
#
#                                y(n) = [1/x(n)]
#                            x(n+1) = 1/x(n)-y(n)
#
#
r2aef:=proc(s)
local max,liste,prod,S,T;
    S :=  evalf(frac(abs(s)));
    max :=  10^(Digits-5);
    prod :=  1;
    liste :=  [trunc(s)];
    while prod <=  max do
        T :=  trunc(1/S);
        S :=  1/T-S;
        liste :=  [op(liste),T];
        prod :=  prod*T
    od;
     RETURN(liste);
end:
###############################################################################
#                         Alternated Egyptian product
#
#                               y(n) = [1/x(n)]
#                            x(n+1) = 1-{(x(n)*y(n)}
#
r2aep:=proc(s)
local max,liste,prod,S,T;
    S :=  evalf(frac(abs(s)));
    max :=  10^(Digits-5);
    prod :=  1;
    liste :=  [trunc(s)];
    while prod <=  max do
        T :=  trunc(1/S);
        S :=  1-frac(S*T);
        liste :=  [op(liste),T];
        prod :=  prod*T
    od;
     RETURN(liste);
end:
###############################################################################
#                               Base k algorithm
#
#                               y(n) = [k*x(n)]
#                              x(n+1) = {k*x(n)}
#
#
r2bk:=proc(s,b) local i,j,v,premier,fin;
v:=abs(evalf(s));
fin:=min(trunc(evalf(Digits/log10(b)))-5,32);
[trunc(v),seq(trunc(v*b**(i+1))-b*trunc(v*b**(i)),i=0..fin)];
end:

###############################################################################
#                               Factorial base
#
#                               y(n) = [n!*x(n)]
#                              x(n+1) = x(n) - y(n)/n!
#
#
r2fb:=proc(s)
local max,liste,prod,S,T,k;
    S :=  evalf(frac(abs(s)));
    max :=  10^(Digits-2);
    prod:=1;k:=1;
    liste :=  [trunc(s)];
    while prod <=  max do
        T :=  trunc(k!*S);
        S :=  S-T/k!;
        liste :=  [op(liste),T];
        k:=k+1;
        prod :=  prod*k;
    od;
     RETURN(liste);
end:
###############################################################################
r2p:=proc(s) local max,liste,prod,S,T;
S:=evalf(abs(s));
max:=10**(Digits-2);
if S>1 then
    prod:=1:
    liste:=[trunc(abs(s))];
    while prod<=max do
          T:=1+floor(1/(S-1));
          S:=S/(T+1)*T;
          liste:=[op(liste),T];
          prod:=prod*T;
    od;
elif S<1 then
    prod:=1:liste:=[trunc(abs(s))];
    while prod<=max do
          T:=1+floor(1/abs(S-1));
          S:=S/(T-1)*T;
          liste:=[op(liste),T];
          prod:=prod*T;
    od;
fi;
RETURN (liste);
end:

Trydev:=proc(s)
local v,w,nn,aa,i,wildguess;
v:=evalf(s);
w:=1/v;
##################
r2cf(v);aa:=";
print(`guessing a gen. fun. for continued fraction : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2ef(v);aa:=";
print(`guessing a gen. gun. for egyptian fraction of x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2ep(v);aa:=";
print(`guessing a gen. gun. for egyptian product of x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2ep(w);aa:=";
print(`guessing a gen. gun. for egyptian product of 1/x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2ef(w);aa:=";
print(`guessing a gen. gun. for egyptian fraction of 1/x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2p(v);aa:=";
print(`guessing a gen. gun. for infinite product of x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2aep(v);aa:=";
print(`guessing a gen. gun. for alt. egyptian product of x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2aef(v);aa:=";
print(`guessing a gen. gun. for alt. egyptian fraction of x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else print(`FAIL`) fi; 
##################
r2aep(w);aa:=";
print(`guessing a gen. gun. for alt. egyptian product of 1/x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2aef(w);aa:=";
print(`guessing a gen. gun. for alt. egyptian fraction of 1/x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2fb(v);aa:=";
print(`guessing a gen. gun. for factorial base of x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else
print(`FAIL`) fi;
##################
r2fb(w);aa:=";
print(`guessing a gen. gun. for factorial base of 1/x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else
print(`FAIL`) fi;
##################
r2bk(v,2);aa:=";
print(`guessing a gen. gun. for binary expansion of x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else
print(`FAIL`) fi;
##################
r2bk(w,2);aa:=";
print(`guessing a gen. gun. for binary expansion of 1/x : `);
wildguess:=listtoratpoly(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else
print(`FAIL`) fi;
##################
end:



Tryhard:=proc(s)
local v,w,nn,aa,i,wildguess;
v:=evalf(s);
w:=1/v;
##################
r2cf(v);aa:=";
print(`guessing a gen. fun. for continued fraction : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2ef(v);aa:=";
print(`guessing a gen. gun. for egyptian fraction of x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2ep(v);aa:=";
print(`guessing a gen. gun. for egyptian product of x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2ep(w);aa:=";
print(`guessing a gen. gun. for egyptian product of 1/x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2ef(w);aa:=";
print(`guessing a gen. gun. for egyptian fraction of 1/x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2p(v);aa:=";
print(`guessing a gen. gun. for infinite product of x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2aep(v);aa:=";
print(`guessing a gen. gun. for alt. egyptian product of x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2aef(v);aa:=";
print(`guessing a gen. gun. for alt. egyptian fraction of x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else print(`FAIL`) fi; 
##################
r2aep(w);aa:=";
print(`guessing a gen. gun. for alt. egyptian product of 1/x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2aef(w);aa:=";
print(`guessing a gen. gun. for alt. egyptian fraction of 1/x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else 
print(`FAIL`) fi;
##################
r2fb(v);aa:=";
print(`guessing a gen. gun. for factorial base of x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else
print(`FAIL`) fi;
##################
r2fb(w);aa:=";
print(`guessing a gen. gun. for factorial base of 1/x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else
print(`FAIL`) fi;
##################
r2bk(v,2);aa:=";
print(`guessing a gen. gun. for binary expansion of x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else
print(`FAIL`) fi;
##################
r2bk(w,2);aa:=";
print(`guessing a gen. gun. for binary expansion of 1/x : `);
wildguess:=guessgf(aa,x);
if length(wildguess)<3*length(aa) then print(wildguess) else
print(`FAIL`) fi;
##################
end: