// Program to calculate e using binary splitting.  There is now a new library
// e.frink which allows resumable calculations and caching of e.  This program
// contains a simpler version of the algorithm for benchmarking and testing.
//
// See:
//    http://numbers.computation.free.fr/Constants/Algorithms/splitting.html
//    http://www.ginac.de/CLN/binsplit.pdf

digits = 10000

if length[ARGS] >= 1
   digits = eval[ARGS@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;
}

setPrecision[digits+3]

//println["k=$k"]
s1 = now[]
start = now[]
p = P[0,k]
q = Q[0,k]
end = now[]
println["Time spent in binary splitting: " + format[end-start, "s", 3]]

start = now[]
e = 1 + (1. * p)/q
end = now[]
println["Time spent in combining operations: " + format[end-start, "s", 3]]

//println[e]                      // Rational number
setPrecision[digits]

start = now[]
e = 1. * e
end = now[]
println["Time spent in floating-point conversion: " + format[end-start, "s", 3]]

start = now[]
es = toString[e]
end = now[]
println["Time spent in radix conversion: " + format[end-start, "s", 3]]

println[e]

e1 = now[]
println["Total time spent:  " + format[e1-s1, "s", 3]]
println[e]

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]
}