/**
   Formula for calculating the Riemann Zeta function.
  
   This follows the very efficient algorithm set out by P. Borwein in 
   "An Efficient Algorithm for the Riemann Zeta Function", Jan. 20, 1995.
  
   http://eprints.cecm.sfu.ca/archive/00000107/
  
   or
   https://vdoc.pub/download/an-efficient-algorithm-for-the-riemann-zeta-function-6p64dg2ligs0
  
   This implements Algorithm 2 of the paper.
  
   As noted by Borwein, "These algorithms do not compete with the
   Riemann-Siegel formula for computations concerning zeros on the critical
   line (Im[s] = 1/2) where multiple low precision evaluations are required."

   Also see:
   "Numerical evaluation of the Riemann Zeta-function", Xavier Gourdon and
   Pascal Sebah
  http://numbers.computation.free.fr/Constants/Miscellaneous/zetaevaluations.pdf
  
   Also see: https://www.hpmuseum.org/cgi-bin/archv021.cgi?read=235199 which
   states: "It turns out that for 12 digits precision one needs about
   n=12-15. This is in accordance with n=1.3d+0.9t suggested by Gourdon and
   Sebah ("d" is the number of digits, "t" is the absolute value of the
   imaginary part of the argument)."  He, however, calculates an extra value
   in the loops of both RiemannD and RiemannZeta.
  
   This means that it'll work around the critical line, but there are known
   faster algorithms if you just need low precision and only work around the
   critical line.
  
   This is the prototype of the (not-yet-implemented) Riemann Zeta function
   in Frink.
*/
class RiemannZeta
{
   // The last value of n we calculated for.
   class var lastN = -1

   // Contains d_k terms.  This is recalculated on demand by recalcRiemannD.
   class var d = new array
   
   // Calculate the value of the Riemann Zeta function at s
   // digits is the approximate number of digits of accuracy.
   class zeta[s, digits=getPrecision[]] :=
   {
      n = ceil[1.3 * digits + 0.9 abs[Im[s]]]
      
      sum = 0
      recalcRiemannD[n]
      
      dn = d@n
      //println["dn is $dn"]
      
      for k = 0 to n-1
      {
         dk = d@k
         sum = sum + (-1)^k (dk - dn)/(k+1)^s
      }

      return sum * (-1/(dn * (1 - 2^(1-s))))
   }

   // n is the approximate number of digits of precision in the result.
   class recalcRiemannD[n] :=
   {
      if n == lastN
         return
      
      sum = 0
      
      for i = 0 to n
      {
         sum = sum + ((n+i-1)! 4^i)/((n-i)! (2i)!)
         d@i = n sum
      }

      lastN = n
   }
}