src/float/transcendental/cl_LF_zeta_int.cc

   1 // zeta().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_F_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/lfloat.h"
  13 #include "cl_LF_tran.h"
  14 #include "cl_LF.h"
  15 #include "cln/integer.h"
  16 #include "cln/abort.h"
  17 #include "cl_alloca.h"
  18
  19 namespace cln {
  20
  21 const cl_LF compute_zeta_exp (int s, uintC len)
  22 {
  23         // Method:
  24         // zeta(s) = 1/(1-2^(1-s)) sum(n=0..infty, (-1)^n/(n+1)^s),
  25         // with convergence acceleration through exp(x), and evaluated
  26         // using the binary-splitting algorithm.
  27         var uintC actuallen = len+2; // 2 Schutz-Digits
  28         var uintL x = (uintL)(0.693148*intDsize*actuallen)+1;
  29         var uintL N = (uintL)(2.718281828*x);
  30         CL_ALLOCA_STACK;
  31         var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
  32         var uintL n;
  33         for (n = 0; n < N; n++) {
  34                 if (n==0) {
  35                         init1(cl_I, args[n].p) (1);
  36                         init1(cl_I, args[n].q) (1);
  37                 } else {
  38                         init1(cl_I, args[n].p) (x);
  39                         init1(cl_I, args[n].q) (n);
  40                 }
  41                 init1(cl_I, args[n].d) (evenp(n)
  42                                         ? expt_pos(n+1,s)
  43                                         : -expt_pos(n+1,s));
  44         }
  45         var cl_LF result = eval_pqd_series(N,args,actuallen);
  46         for (n = 0; n < N; n++) {
  47                 args[n].p.~cl_I();
  48                 args[n].q.~cl_I();
  49                 args[n].d.~cl_I();
  50         }
  51         result = shorten(result,len); // verkürzen und fertig
  52         // Zum Schluss mit 2^(s-1)/(2^(s-1)-1) multiplizieren:
  53         return scale_float(result,s-1) / (ash(1,s-1)-1);
  54 }
  55 // Bit complexity (N = len): O(log(N)^2*M(N)).
  56
  57 const cl_LF compute_zeta_cvz1 (int s, uintC len)
  58 {
  59         // Method:
  60         // zeta(s) = 1/(1-2^(1-s)) sum(n=0..infty, (-1)^n/(n+1)^s),
  61         // with Cohen-Villegas-Zagier convergence acceleration.
  62         var uintC actuallen = len+2; // 2 Schutz-Digits
  63         var uintL N = (uintL)(0.39321985*intDsize*actuallen)+1;
  64         var cl_I fterm = 2*(cl_I)N*(cl_I)N;
  65         var cl_I fsum = fterm;
  66         var cl_LF gterm = cl_I_to_LF(fterm,actuallen);
  67         var cl_LF gsum = gterm;
  68         var uintL n;
  69         // After n loops
  70         //   fterm = (N+n)!N/(2n+2)!(N-n-1)!*2^(2n+2), fsum = ... + fterm,
  71         //   gterm = S_n*fterm, gsum = ... + gterm.
  72         for (n = 1; n < N; n++) {
  73                 fterm = exquopos(fterm*(2*(cl_I)(N-n)*(cl_I)(N+n)),(cl_I)(2*n+1)*(cl_I)(n+1));
  74                 fsum = fsum + fterm;
  75                 gterm = The(cl_LF)(gterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
  76                 if (evenp(n))
  77                         gterm = gterm + cl_I_to_LF(fterm,actuallen)/expt_pos(n+1,s);
  78                 else
  79                         gterm = gterm - cl_I_to_LF(fterm,actuallen)/expt_pos(n+1,s);
  80                 gsum = gsum + gterm;
  81         }
  82         var cl_LF result = gsum/cl_I_to_LF(1+fsum,actuallen);
  83         result = shorten(result,len); // verkürzen und fertig
  84         // Zum Schluss mit 2^(s-1)/(2^(s-1)-1) multiplizieren:
  85         return scale_float(result,s-1) / (ash(1,s-1)-1);
  86 }
  87 // Bit complexity (N = len): O(N^2).
  88
  89 const cl_LF compute_zeta_cvz2 (int s, uintC len)
  90 {
  91         // Method:
  92         // zeta(s) = 1/(1-2^(1-s)) sum(n=0..infty, (-1)^n/(n+1)^s),
  93         // with Cohen-Villegas-Zagier convergence acceleration, and
  94         // evaluated using the binary splitting algorithm.
  95         var uintC actuallen = len+2; // 2 Schutz-Digits
  96         var uintL N = (uintL)(0.39321985*intDsize*actuallen)+1;
  97         CL_ALLOCA_STACK;
  98         var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
  99         var uintL n;
 100         for (n = 0; n < N; n++) {
 101                 init1(cl_I, args[n].p) (2*(cl_I)(N-n)*(cl_I)(N+n));
 102                 init1(cl_I, args[n].q) ((cl_I)(2*n+1)*(cl_I)(n+1));
 103                 init1(cl_I, args[n].d) (evenp(n)
 104                                         ? expt_pos(n+1,s)
 105                                         : -expt_pos(n+1,s));
 106         }
 107         var cl_pqd_series_result sums;
 108         eval_pqd_series_aux(N,args,sums);
 109         // Here we need U/(1+S) = V/D(Q+T).
 110         var cl_LF result =
 111           cl_I_to_LF(sums.V,actuallen) / The(cl_LF)(sums.D * cl_I_to_LF(sums.Q+sums.T,actuallen));
 112         for (n = 0; n < N; n++) {
 113                 args[n].p.~cl_I();
 114                 args[n].q.~cl_I();
 115                 args[n].d.~cl_I();
 116         }
 117         result = shorten(result,len); // verkürzen und fertig
 118         // Zum Schluss mit 2^(s-1)/(2^(s-1)-1) multiplizieren:
 119         return scale_float(result,s-1) / (ash(1,s-1)-1);
 120 }
 121 // Bit complexity (N = len): O(log(N)^2*M(N)).
 122
 123 // Timings of the above algorithm, on an i486 33 MHz, running Linux.
 124 //    s                5                           15
 125 //    N     sum_exp sum_cvz1 sum_cvz2   sum_exp sum_cvz1 sum_cvz2
 126 //    10       2.04    0.09    0.17        8.0     0.11    0.49
 127 //    25       8.6     0.30    0.76       30.6     0.37    2.36
 128 //    50      25.1     0.92    2.49       91.1     1.15    7.9
 129 //   100               2.97    8.46                3.75   24.5
 130 //   250              16.7    36.5                21.7   108
 131 //   500              64.2   106                  85.3   295
 132 //  1000             263     285                 342     788
 133 // asymp.      FAST    N^2     FAST        FAST    N^2     FAST
 134 //
 135 // The break-even point between cvz1 and cvz2 seems to grow linearly with s.
 136
 137 const cl_LF zeta (int s, uintC len)
 138 {
 139         if (!(s > 1))
 140                 cl_abort();
 141         if (len < 280*(uintL)s)
 142                 return compute_zeta_cvz1(s,len);
 143         else
 144                 return compute_zeta_cvz2(s,len);
 145 }
 146 // Bit complexity (N = len): O(log(N)^2*M(N)).
 147
 148 }  // namespace cln