12 #include "cln/lfloat.h"
13 #include "cl_LF_tran.h"
15 #include "cln/integer.h"
16 #include "cl_alloca.h"
20 const cl_LF zeta3 (uintC len)
24 // | ----- (n + 1) 2 |
25 // 1 | \ (-1) (205 n - 160 n + 32)|
26 // - | ) ---------------------------------|
28 // | ----- n binomial(2 n, n) |
31 // The formula used to compute Zeta(3) has reference in the paper
32 // "Hypergeometric Series Acceleration via the WZ method" by
33 // T. Amdeberhan and Doron Zeilberger,
34 // Electronic J. Combin. 4 (1997), R3.
36 // Computation of the sum:
37 // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
38 // with appropriate N, and
39 // a(n) = 205*n^2+250*n+77, b(n) = 1,
40 // p(0) = 1, p(n) = -n^5 for n>0, q(n) = 32*(2n+1)^5.
41 var uintC actuallen = len+2; // 2 Schutz-Digits
42 var uintC N = ceiling(actuallen*intDsize,10);
43 // 1024^-N <= 2^(-intDsize*actuallen).
45 var cl_I* av = (cl_I*) cl_alloca(N*sizeof(cl_I));
46 var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
47 var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
49 for (n = 0; n < N; n++) {
50 init1(cl_I, av[n]) (205*square((cl_I)n) + 250*(cl_I)n + 77);
52 init1(cl_I, pv[n]) (1);
54 init1(cl_I, pv[n]) (-expt_pos(n,5));
55 init1(cl_I, qv[n]) (expt_pos(2*n+1,5)<<5);
57 var cl_pqa_series series;
59 series.pv = pv; series.qv = qv; series.qsv = NULL;
60 var cl_LF sum = eval_rational_series(N,series,actuallen);
61 for (n = 0; n < N; n++) {
66 return scale_float(shorten(sum,len),-1);
68 // Bit complexity (N := len): O(log(N)^2*M(N)).
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
84 // asymp. FAST N^2 FAST FAST
85 // (FAST means O(log(N)^2*M(N)))