src/float/transcendental/cl_LF_zeta3.cc

   1 // cl_zeta3().
   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 "cl_lfloat.h"
  13 #include "cl_LF_tran.h"
  14 #include "cl_LF.h"
  15 #include "cl_integer.h"
  16 #include "cl_alloca.h"
  17
  18 const cl_LF cl_zeta3 (uintC len)
  19 {
  20         // Method:
  21         //            /infinity                                  \
  22         //            | -----       (n + 1)       2              |
  23         //        1   |  \      (-1)        (205 n  - 160 n + 32)|
  24         //        -   |   )     ---------------------------------|
  25         //        2   |  /              5                 5      |
  26         //            | -----          n  binomial(2 n, n)       |
  27         //            \ n = 1                                    /
  28         //
  29         // The formula used to compute Zeta(3) has reference in the paper
  30         // "Acceleration of Hypergeometric Series via the WZ method" by
  31         // T. Amdeberhan and Doron Zeilberger, to appear in the Electronic
  32         // Journal of Combinatorics [Wilf Festschrift Volume].
  33         //
  34         // Computation of the sum:
  35         // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
  36         // with appropriate N, and
  37         //   a(n) = 205*n^2+250*n+77, b(n) = 1,
  38         //   p(0) = 1, p(n) = -n^5 for n>0, q(n) = 32*(2n+1)^5.
  39         var uintC actuallen = len+2; // 2 Schutz-Digits
  40         var uintL N = ceiling(actuallen*intDsize,10);
  41         // 1024^-N <= 2^(-intDsize*actuallen).
  42         CL_ALLOCA_STACK;
  43         var cl_I* av = (cl_I*) cl_alloca(N*sizeof(cl_I));
  44         var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  45         var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  46         var uintL n;
  47         for (n = 0; n < N; n++) {
  48                 init1(cl_I, av[n]) (205*square((cl_I)n) + 250*(cl_I)n + 77);
  49                 if (n==0)
  50                         init1(cl_I, pv[n]) (1);
  51                 else
  52                         init1(cl_I, pv[n]) (-expt_pos(n,5));
  53                 init1(cl_I, qv[n]) (expt_pos(2*n+1,5)<<5);
  54         }
  55         var cl_pqa_series series;
  56         series.av = av;
  57         series.pv = pv; series.qv = qv; series.qsv = NULL;
  58         var cl_LF sum = eval_rational_series(N,series,actuallen);
  59         for (n = 0; n < N; n++) {
  60                 av[n].~cl_I();
  61                 pv[n].~cl_I();
  62                 qv[n].~cl_I();
  63         }
  64         return scale_float(shorten(sum,len),-1);
  65 }
  66 // Bit complexity (N := len): O(log(N)^2*M(N)).
  67
  68 // Timings of the above algorithm, on an i486 33 MHz, running Linux.
  69 //    N   sum_exp sum_cvz1 sum_cvz2 hypgeom
  70 //    10     1.17    0.081   0.125   0.013
  71 //    25     5.1     0.23    0.50    0.045
  72 //    50    15.7     0.66    1.62    0.14
  73 //   100    45.5     1.93    5.4     0.44
  74 //   250   169      13.1    25.1     2.03
  75 //   500   436      56.5    70.6     6.44
  76 //  1000           236     192      18.2
  77 //  2500                            78.3
  78 //  5000                           202
  79 // 10000                           522
  80 // 25000                          1512
  81 // 50000                          3723
  82 // asymp.    FAST     N^2    FAST    FAST
  83 // (FAST means O(log(N)^2*M(N)))