// Program to calculate e using binary splitting.  This is very much
// faster than the simple method used in e.frink.
//
// See:
//    http://numbers.computation.free.fr/Constants/Algorithms/splitting.html
//    http://www.ginac.de/CLN/binsplit.pdf

digits = 10000
setPrecision[digits+3]

sum = 0

// Find number of terms to calculate.  ln[x!] = ln[1] + ln[2] + ... + ln[x]
k = 1;
logFactorial = 0.;
logMax = digits * ln[10]
while (logFactorial < logMax)
{
   logFactorial = logFactorial + ln[k];
   k = k + 1;
}

println["k=$k"]

start = now[];
e = 1 + P[0,k]/Q[0,k]
end = now[];

println["Calculation complete."]

println[e]                      // Rational number
println[1. * e]                 // Force to floating-point

println["Time is " + (end-start)]

P[a,b] :=
{
   if (b-a) == 1
      return 1
   m = (a+b) div 2
   r = P[a,m] Q[m,b] + P[m,b]
   return r
}

Q[a,b] :=
{
   if (b-a) == 1
      return b
   m = (a+b) div 2
   return Q[a,m] Q[m,b]
}