#
x=0.001
DIGITS   500  1000  2000  5000  10000
x        
exp         0.9    1.6      5.4    43.9    248
exp1       3.1    6.0    20.4    154     828       [without ln(2.0)]
exp2       1.3    2.4      7.9     80      461
exp3       1.3    2.4      8.2    68.7     300
ln(2.0)    1.3     2.4     8.0    76.9     418
PI           0.8     1.0     1.4     4.7      15.9
fE           5.4    10.6    19.5   47.1    97.8
ln(3.0)     1.0     2.3     7.8    77.6
ln1(3.0)   2.2     3.6     9.5    54       218
ln2(3.0)   1.4    2.4      7.4    58.5
#

if type(lntwo)<>DOM_PROC then
lntwo:=proc() local d;
begin 
   if DIGITS>lntwo(0) then
      d:=DIGITS; DIGITS:=max(d,round(1.1*lntwo(0))); print("lntwo",DIGITS); 
      lntwo(0):=DIGITS; lntwo(1):=ln(2.0); DIGITS:=d end_if;
   float(lntwo(1))
end_proc:
lntwo(0):=0: # current precision #
end_if:

if type(fPI)<>DOM_PROC then
fPI:=proc() local d;
begin
   if DIGITS>fPI(0) then
      d:=DIGITS; DIGITS:=max(d,round(1.1*fPI(0))); print("PI",DIGITS);
      fPI(0):=DIGITS; fPI(PI):=float(PI); DIGITS:=d end_if;
   float(fPI(PI))
end_proc:
fPI(0):=0: # current precision #
end_if:

fE:=proc() local u,s,t,p,k,st;
begin
   p:=round(3.3219281*DIGITS+3);
   u:=2^p; s:=u; k:=1; t:=u; 
st:=time();
   repeat
      t:=t div k; s:=s+t; k:=k+1; t:=t div k; s:=s+t; k:=k+1; t:=t div k; s:=s+t; k:=k+1;
   until t=0 end_repeat;
print("loop=",time()-st);
   float(s)/float(u)
end_proc:

# using exp(x)=(1+r+r^2/2!+r^3/3!+...)^(2^K)*2^n
   where x=n*ln(2)+(2^K)*r
   number of operations = O(K+DIGITS*ln(10)/K/ln(2))
#
exp1:=proc(x,K) local d,n,r,s,t,k,e,dd,DIGITS;
begin
   x:=float(x);
   d:=DIGITS; DIGITS:=9; n:=round(x/0.6931471805);
   if args(0)=1 then K:=round(sqrt(5*d)); print(K) end_if;
   dd:=d+round(K*0.6931471805); DIGITS:=dd;
   r:=(if n<>0 then n*lntwo() else 0 end_if);
   r:=x-r; r:=r/2^K;
   # r^k/k!+r^(k+1)/(k+1)!+... <= r^k/k!*(1+r/(k+1)+...) <= r^k/k!/(1-r/(k+1)) #
   k:=0; t:=1; s:=1.0; DIGITS:=9; e:=10.0^(-d-3-K*0.3);
   repeat
      DIGITS:=dd;
      k:=k+1;
      t:=t*r/k;
      s:=s+t;
   until (DIGITS:=9;abs(t)<=s*e*(1-r/k)) end_repeat;
   DIGITS:=dd; 
   for k from 1 to K do s:=s*s end_for;
   if n<>0 then s:=s*2^n end_if; DIGITS:=d; float(s)
end_proc:

# using Newton's method with ln #
exp2:=proc(x) local d;
begin
   d:=DIGITS; DIGITS:=9; exp(x); DIGITS:=d;
   newton(subs(fun(args(1)*(x+1-ln(args(1)))),hold(x)=x),%2)
end_proc:

# using Newton's method with best ln #
exp3:=proc(x) local d;
begin
   d:=DIGITS; DIGITS:=9; exp(x); DIGITS:=d;
   newton(subs(fun(args(1)*(x+1-Ln(args(1)))),hold(x)=x),%2)
end_proc:

# iterates f from x #
newton:=proc(f,x) local oldx,l,k;
begin
   l:=[DIGITS]; while l[1]>9 do l:=[(l[1]+1) div 2].l end_while;
   for k from 2 to nops(l) do DIGITS:=l[k]; oldx:=x; x:=f(x) end_for; x
end_proc:

Ln:=fun(
   (if DIGITS<2457 then ln elif DIGITS<4100 then ln2 else ln1 end_if)(args(1))
):

ln2:=proc(x) local d; # using Newton's method #
begin
   d:=DIGITS; DIGITS:=9; ln(x); DIGITS:=d;
   newton(subs(fun(args(1)-1+x/exp(args(1))),hold(x)=x),%2)
end_proc:

AG:=proc(a,b) 
begin
   repeat
      (a+b)/2.0; b:=(a*b)^(1/2); a:=%2
   until b-a=0.0 end_repeat;
   b
end_proc:

# search m such that x*2^m>10^(DIGITS/2), i.e. ln(x)+m*ln(2)>DIGITS/2*ln(10) #
ln1:=proc(x) local d,m,s,k;
begin
   d:=DIGITS; x:=float(x);
   DIGITS:=9; m:=(1.1513*d-ln(x))/0.693; m:=ceil(m+ln(m));
   DIGITS:=(d:=d+3); x:=float(x); x:=x*2^m;
   x:=fPI()/2/AG(1.0,4/x)-m*lntwo();
   DIGITS:=d-3; float(x)
end_proc: