src/float/transcendental/cl_LF_zeta3.cc

   1 // zeta3().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "float/transcendental/cl_F_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/lfloat.h"
  13 #include "float/transcendental/cl_LF_tran.h"
  14 #include "float/lfloat/cl_LF.h"
  15 #include "cln/integer.h"
  16 #include "base/cl_alloca.h"
  17
  18 namespace cln {
  19
  20 const cl_LF zeta3 (uintC len)
  21 {
  22         struct rational_series_stream : cl_pqa_series_stream {
  23                 uintC n;
  24                 static cl_pqa_series_term computenext (cl_pqa_series_stream& thisss)
  25                 {
  26                         var rational_series_stream& thiss = (rational_series_stream&)thisss;
  27                         var uintC n = thiss.n;
  28                         var cl_pqa_series_term result;
  29                         if (n==0) {
  30                                 result.p = 1;
  31                         } else {
  32                                 result.p = -expt_pos(n,5);
  33                         }
  34                         result.q = expt_pos(2*n+1,5)<<5;
  35                         result.a = 205*square((cl_I)n) + 250*(cl_I)n + 77;
  36                         thiss.n = n+1;
  37                         return result;
  38                 }
  39                 rational_series_stream ()
  40                         : cl_pqa_series_stream (rational_series_stream::computenext),
  41                           n (0) {}
  42         } series;
  43         // Method:
  44         //            /infinity                                  \
  45         //            | -----       (n + 1)       2              |
  46         //        1   |  \      (-1)        (205 n  - 160 n + 32)|
  47         //        -   |   )     ---------------------------------|
  48         //        2   |  /              5                 5      |
  49         //            | -----          n  binomial(2 n, n)       |
  50         //            \ n = 1                                    /
  51         //
  52         // The formula used to compute Zeta(3) has reference in the paper
  53         // "Hypergeometric Series Acceleration via the WZ method" by
  54         // T. Amdeberhan and Doron Zeilberger,
  55         // Electronic J. Combin. 4 (1997), R3.
  56         //
  57         // Computation of the sum:
  58         // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
  59         // with appropriate N, and
  60         //   a(n) = 205*n^2+250*n+77, b(n) = 1,
  61         //   p(0) = 1, p(n) = -n^5 for n>0, q(n) = 32*(2n+1)^5.
  62         var uintC actuallen = len+2; // 2 guard digits
  63         var uintC N = ceiling(actuallen*intDsize,10);
  64         // 1024^-N <= 2^(-intDsize*actuallen).
  65         var cl_LF sum = eval_rational_series<false>(N,series,actuallen,actuallen);
  66         return scale_float(shorten(sum,len),-1);
  67 }
  68 // Bit complexity (N := len): O(log(N)^2*M(N)).
  69
  70 // Timings of the above algorithm, on an i486 33 MHz, running Linux.
  71 //    N   sum_exp sum_cvz1 sum_cvz2 hypgeom
  72 //    10     1.17    0.081   0.125   0.013
  73 //    25     5.1     0.23    0.50    0.045
  74 //    50    15.7     0.66    1.62    0.14
  75 //   100    45.5     1.93    5.4     0.44
  76 //   250   169      13.1    25.1     2.03
  77 //   500   436      56.5    70.6     6.44
  78 //  1000           236     192      18.2
  79 //  2500                            78.3
  80 //  5000                           202
  81 // 10000                           522
  82 // 25000                          1512
  83 // 50000                          3723
  84 // asymp.    FAST     N^2    FAST    FAST
  85 // (FAST means O(log(N)^2*M(N)))
  86
  87 }  // namespace cln