share/numeric/bffac.mac

   1 /* Author: Bill Gosper */
   2
   3 /*
   4   Initial version generated from bffac.lisp with the commands
   5
   6   load("bffac.lisp");
   7   bffac_functions:[OBFAC,AZETB,VONSCHTOONK,DIVRLST,BURN,OBZETA,BZETA,
   8   BFPSI,BFFAC,BFHZETA,CBFFAC,BFPSI0,BFZETA];
   9   apply(stringout,cons("bffac.mac",map('fundef,bffac_functions))),grind:true;
  10 */
  11
  12
  13 prodhack_product (expr, i%, i%0, i%1) :=
  14     (expr: subst ('i%, i%, expr),
  15     if i%0 > i%1 then
  16         subst (i%, 'i%, 1 / product (expr, i%, i%1 + 1, i%0 - 1))
  17     else
  18         subst (i%, 'i%, product (expr, i%, i%0, i%1)));
  19
  20 obfac(z,fpprec):=block([y:z-entier(z),x:entier(7*fpprec/3),m:0],
  21       for k from entier(26*fpprec/3) step -1 thru 1 do m:x*(1/(y+k)-m/k),
  22       bfloat(m)*x^y*prodhack_product(z-i+1,i,1,entier(z)));
  23
  24 azetb(s,fpp):=block([fpprec:15],
  25       if s > 0
  26           then block([m:0,
  27                       p:print(entier(log(
  28                                       bfloat(2*%pi)^s
  29                                        *bfloat(10)^(fpp:entier(fpp))
  30                                        /bfloat(%pi*(s-1)!))
  31                                       /bfloat(2*%pi)))],fpprec:fpp,
  32                      for k from entier(ev(%pi*p,numer)) step -1 thru 0 do
  33                          m:(s+2*k)*(s+2*k-1)*m/((2*k+1)*(2*k+2)*p^2)
  34                            +p*bern(2*k)/(s-1),(1/2+m)/p^s+sum(k^-s,k,1,p-1),
  35                      bfloat(%\%))
  36           else (if s = 0 then -1/2
  37                     else ((2*%pi)^s*sin(%pi*s/2)*bffac(-s,fpprec)
  38                                    *bzeta(1-s,fpprec)
  39                            /%pi,bfloat(%\%))));
  40
  41 /* -----------------------------------------------------------------------------
  42  * Function burn:
  43  *
  44  * Calculate an approximation of Bernoulli numbers for even integers using
  45  * the formula:
  46  *
  47  *
  48  *                          n - 1  1 - 2 n
  49  *                     (- 1)      2        zeta(2 n) (2 n)!
  50  *    bernoulli(2 n) = ------------------------------------
  51  *                                       2 n
  52  *                                    %pi
  53  *
  54  * burn(n) calculates a rational number, which is a very good approximation of
  55  * the Bernoulli number for big integers n. With bern the time to calculate a
  56  * bernoulli number grows very fast for the first call.
  57  * The function burn is much faster. The following table shows the computation
  58  * time in seconds for bern (first call) and burn. Furthermore, the relative
  59  * error (bern(n)-burn(n))/bern(n) is shown:
  60  *
  61  *     n     bern(n)     burn(n)     relative Error
  62  *    256     0,1200 s    0,0040 s      3.41b-302
  63  *    512     0,7800 s    0,0120 s      1.94b-758
  64  *   1024     6,2320 s    0,0320 s      6.38b-1822
  65  *   2048    59,5040 s    0,1360 s      2.50b-4258
  66  *   4096     -           0,9720 s
  67  *   8192     -           7,3680 s
  68  *
  69  * Burn invokes the approximation for even integers n > 255. For odd integers
  70  * and n <= 255 the function bern is called.
  71  *
  72  * The functions vonschtoonk and divrlst are helper functions for burn.
  73  * -------------------------------------------------------------------------- */
  74
  75 vonschtoonk(p):=apply("*",
  76             subst(["*" = lambda([[l]],1),"^" = lambda([[l]],1)],
  77                   factor(1+divrlst(p))));
  78
  79 divrlst(p):=(block([temp:0],
  80         temp: expand(ratsimp(subst("^" = lambda([a,b],
  81                                                (1-part(a)^(b/2+1))
  82                                                 /(1-part(a))),factor(p^2)))),
  83     if atom(temp) then cons(temp,[]) else args(temp),
  84         ev(%\%,part)));
  85
  86 burn(p):=if evenp(p) and p>255
  87       then block([d:vonschtoonk(p)],
  88                  block([fpprec:entier(ev(
  89                                        (log(8*%pi*p)/2+p*(log(p/(2*%pi))-1)
  90                                                       +2*log(d))
  91                                         /log(10),numer))
  92                                +1],
  93                        block([pi:bfloat(%pi)],
  94                              entier(2*d*p!*bfzeta(p,fpprec)/(2^p*pi^p)+1/2)
  95                               /d)))
  96       *(-1)^(p/2-1) else bern(p);
  97
  98 obzeta(s,fpprec):=if s > 0
  99         then block([m:0,
 100                     p:entier(ev((log(2*%pi)*s+log(10)*fpprec-log(%pi*(s-1)!))
 101                                  /(2*%pi),numer))],
 102                    for k from entier(ev(%pi*p,numer)) step -1 thru 0 do
 103                        m:(s+2*k)*(s+2*k-1)*m/((2*k+1)*(2*k+2)*p^2)
 104                          +p*bern(2*k)/(s-1),(1/2+m)/p^s+sum(k^-s,k,1,p-1),
 105                    bfloat(%\%))
 106         else (if s = 0 then -1/2
 107                   else ((2*%pi)^s*sin(%pi*s/2)*bffac(-s,fpprec)
 108                                         *bzeta(1-s,fpprec)
 109                                 /%pi,bfloat(%\%)));
 110
 111 bzeta(s,fpp):=if s > 0
 112        then block([fpprec:7,m:0,t:0.0b0,p,
 113                    k:2*max(1,
 114                            entier(ev(3/2-(s-1/2)*log(s/(2*%pi))+log(10)*fpp
 115                                         -log(%pi),numer)
 116                                    /2))],
 117                   p:1+entier(bfloat((2*10^fpp*product(s+k-i-2,i,0,k-2)
 118                                      /(2*3.14159265b0)^k)
 119                                      ^(1.0b0/(s+k-1)))),fpprec:fpp,
 120                   for k from k/2 step -1 thru 0 do
 121                       m:(s+2*k)*(s+2*k-1)*m/((2*k+1)*(2*k+2)*p^2)
 122                         +p*bern(2*k)/(s-1),
 123                   if s < 25 then t:sum(k^-s,k,1,p-1)
 124                       else for k from p-1 step -1 thru 1 do
 125                                (fpprec:fpp-entier(ev(s*log(k)/log(10),numer)),
 126                                 t:t+k^-s),bfloat((1/2+m)/p^s+t))
 127        else (if s = 0 then -1/2
 128                  else ((2*%pi)^s*sin(%pi*s/2)*bffac(-s,fpprec)
 129                                 *bzeta(1-s,fpprec)
 130                         /%pi,bfloat(%\%)));
 131
 132 bfpsi(n,z,fpprec):=if equal(n,0) then bfpsi0(z,fpprec)
 133        else bfhzeta(n+1,z,fpprec)*(-1)^(n-1)*n!;
 134
 135 bffac(z,fpprec):=if z < 0 then bfloat(%pi*z/sin(%pi*z))/bffac(-z,fpprec)
 136        else block([k:2*(entier(0.41*fpprec)+1)],
 137                   block([m:1,y:(z+k)^2,x:0.0b0],
 138                         for i thru k/2 do
 139                             (m:m*(z+2*i-1)*(z+2*i),
 140                              x:(x+bern(k-2*i+2)/((k-2*i+1)*(k-2*i+2)))/y),
 141                         bfloat(sqrt(2)*sqrt(%pi)*sqrt(z+k)
 142                                       *%e^((z+k)*(log(z+k)+x-1)))
 143                          /m));
 144
 145 /* Hurwitz Zeta function
 146  *
 147  * bfhzeta(s,h) = sum(1/(h+k)^s, k, 0, inf)
 148  */
 149 bfhzeta(s,h,fpprec):=if s > 0
 150          then block([m:0,p,
 151                      q:entier(ev((-log(%pi*(s-1)!)+log(2*%pi)*s
 152                                                   +log(10)*fpprec)
 153                                   /(2*%pi)
 154                                   -h+1,numer))
 155                        +h],p:max(q,h+1),
 156                     for k from entier(ev(%pi*q,numer)) step -1 thru 0 do
 157                         m:m*(s+2*k-1)*(s+2*k)/((2*k+1)*(2*k+2)*(p-1)^2)
 158                           +bern(2*k)*(p-1)/(s-1),
 159                     (m+1/2)/(p-1)^s+sum(1/(k+h)^s,k,0,entier(p-h-1.9b0)),
 160                     bfloat(%\%))
 161          else (if s = 0 then 1/2-h else funmake('bfhzeta,[s,h,fpprec]));
 162
 163 cbffac(z,fpprec):=rectform(
 164        if realpart(z:rectform(z)) < 0
 165            then bfloat(%pi*z/sin(%pi*z))/cbffac(-z,fpprec)
 166            else block([k:2*(entier(0.41*fpprec)+1)],
 167                       block([m:1,y:rectform(1/(z+k)^2),x:0.0b0],
 168                             for i thru k/2 do
 169                                 (m:expand(m*(z+2*i-1)*(z+2*i)),
 170                                  x:expand(
 171                                    (x+bern(k-2*i+2)/((k-2*i+1)*(k-2*i+2)))
 172                                     *y)),
 173                             bfloat(sqrt(2)*sqrt(%pi)*sqrt(z+k)
 174                                           *%e^((z+k)*(log(z+k)+x-1)))
 175                              /m)));
 176
 177 /* psi[0](z) - digamma function
 178  *
 179  * A&S 6.3.5:  Recurrence formula
 180  *
 181  *   psi[0](1+z) = 1/z+psi[0](z)
 182  *
 183  * A&S 6.3.6:  Recurrence formula
 184  *
 185  *   psi[0](n+z) = 1/(z+n-1)+1/(z+n-2)+...+1/(z+2)+1/(1+z)+psi[0](1+z)
 186  *
 187  * A&S 6.3.7:  Reflection formula
 188  *
 189  *   psi[0](1-z) = psi[0](z) + %pi*cot(%pi*z)
 190  *
 191  * A&S 6.3.18:  Asymptotic formula
 192  *
 193  *   psi[0](z) ~ log(z) - sum(bern(2*n)/(2*n*z^(2*n)), n, 1, inf)
 194  *
 195  * So use reflection formula if z < 0.  For z > 0, use the recurrence
 196  * formula to increase the argument and then apply the asymptotic formula.
 197  *
 198  * For small z with absolute value < 1E-6 use formula 6.3.5
 199  *
 200  */
 201 bfpsi0(z,fpprec):=if z < 0 then bfloat(%pi*cot(%pi*(-z)))+bfpsi0(1-z,fpprec)
 202         elseif 0 < realpart(z) and cabs(z) < 1E-6 then bfloat(rectform(bfpsi0(1+z,fpprec) - 1/z))
 203         else block([k:2*(entier(0.41*fpprec)+1)],
 204                    block([m:0,y:(z+k)^2,x:0.0b0],
 205                          for i thru k/2 do
 206                              (m:rectform(1/(z+2*i-1)+1/(z+2*i-2)+m),
 207                               x:rectform((x+bern(k-2*i+2)/(k-2*i+2))/y)),
 208                          bfloat(rectform(log(z+k)-1/(2*(z+k))-x-m))));
 209
 210 bfzeta (s, fpp) :=
 211   block ([sigma: realpart (s), t: imagpart (s)],
 212     if is (not (numberp (sigma)) or
 213            not (numberp (t)) or
 214            not (integerp (fpp) and fpp > 0)) then 'bfzeta(s,fpp)
 215     else (
 216       if is (s = 0) then bfloat (-1/2)
 217       elseif is (s = 1) then inf
 218       elseif is (t = 0 and sigma > 0) then bfzeta_positive (s, fpp)
 219       elseif is (t # 0 and sigma >= 1/2) then bfzeta_borwein (s, fpp)
 220       else
 221         /*
 222           bfzeta_negative calls bfzeta with 1-s. However, we know that
 223           sigma < 1/2, so 1-sigma > 1/2 and we won't recurse infinitely.
 224         */
 225         bfzeta_negative (s, fpp)))$
 226
 227 /*
 228   A&S 23.2.6
 229 */
 230 bfzeta_negative (s,fpp):=
 231   block([fpprec: fpp],
 232     rectform (bfloat (bfzeta(1-s,fpprec) * %pi^(s-1) * 2^s *
 233                       bffac(-s,fpprec) * sin(%pi*s/2))))$
 234
 235 /*
 236   Calculates zeta(s) for positive real s using an Euler-Maclaurin
 237   expansion (The fastest of the algorithms, but only works for real s)
 238 */
 239 bfzeta_positive (s, fpp) :=
 240   block ([fpprec:7,m:0,t:0.0b0,p,
 241           k:2*max(1,
 242                   entier(ev(-(s-1/2)*log(s/(2*%pi))
 243                             +log(10)*fpp-log(%pi)+3/2,numer)/2))],
 244           p:entier(bfloat(2^(1/(s+k-1))*10^(fpp/(s+k-1))
 245                    *product(s+k-i-2,i,0,k-2)
 246                    ^(1/(s+k-1))
 247                    /6.2831853b0^(k/(s+k-1))))+1,
 248           fpprec:fpp,
 249           for k from k/2 step -1 thru 0 do
 250             m:m*(s+2*k-1)*(s+2*k)/((2*k+1)*(2*k+2)*p^2)
 251               +bern(2*k)*p/(s-1),
 252           if s < 25 then t:sum(1/k^s,k,1,p-1)
 253           else for k from p-1 step -1 thru 1 do
 254                  (fpprec:fpp-entier(ev(log(k)*s/log(10),numer)),
 255                   t:t+1/k^s),
 256           bfloat(t+(m+1/2)/p^s))$
 257
 258 /*
 259   The bfzeta_borwein function calculates zeta(S) to within 10^(-FPP)
 260   using the methods described in Borwein, "An Efficient Algorithm for
 261   the Riemann Zeta Function" (MR1777614), described slightly
 262   differently in Gourdon and Sebah, "Numerical evaluation of the
 263   Riemann Zeta-function" (we use their notation).
 264
 265   Write s = sigma + %i*t. The expansion is valid for sigma >= 1/2.
 266
 267   The bound given in Algorithm 2 of the papers shows that, in order to
 268   get an error below 10^(-N) when evaluating zeta(s) = zeta(sigma+it)
 269   with sigma >= 1/2, we need n terms, where
 270
 271     log(3+sqrt(8))*n
 272      >=
 273     log(3) + log(1+2*abs(t)) + %pi/2*abs(t) + N*log(10) - log(abs(1-2^(1-s)))
 274
 275   To calculate the sum, we need binomial coefficients of the form
 276
 277     b[i] = n*(n+i-1)!4^i/((n-i)!(2i)!)
 278
 279   with recurrence given by
 280
 281     b[i] = b[i+1] * (2*i+1)*(i+1)/(2*(n+i)*(n-i))
 282
 283   where b[n] = 4^n / 2. Then let d[k] = sum (b[i], i, k, n) (with the
 284   very easy recurrence that d[k] = n*b[k] + d[k+1]) and
 285
 286     zeta(s) = (d[0]*(1-2^(1-s)))^(-1) *
 287               sum ((-1)^(k-1) * d[k] / k^s) + err
 288
 289 */
 290 bfzeta_borwein_terms_needed (s, t, N) :=
 291   ceiling (bfloat (
 292       (log (3) + log (1+2*abs (t)) + %pi/2*abs(t)
 293                + N*log(10) - log(abs(1-2^(1-s))))
 294       / 1.76b0))$
 295
 296 bfzeta_borwein (s, fpp) := block (
 297   [t: bfloat (imagpart (s)), sigma: bfloat (realpart (s)), n],
 298   n: bfzeta_borwein_terms_needed (s, t, fpp+1),
 299   /* b is the binomial coefficient; d is dk */
 300   block ([fpprec: ceiling (log (n) + fpp), acc: [0,0],
 301           b: 4^n/2, d: 0, lbk],
 302     for k:n step -1 thru 1 do (
 303       tlbk: -t * log (bfloat (k)),
 304       d: d + b,
 305       b: b * (2*k-1)*k / (2*(n+k-1)*(n-k+1)),
 306       acc: acc + (-1)^(k-1)*d*k^(-sigma)*[cos(tlbk), sin(tlbk)]),
 307     d: d + b,
 308     rectform (bfloat ((d * (1-(2^(1-s))))^-1) *
 309               (first (acc) + %i*second (acc)))))$