4 #include "base/cl_sysdep.h"
7 #include "float/transcendental/cl_F_tran.h"
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"
20 const cl_LF zeta3 (uintC len)
22 struct rational_series_stream : cl_pqa_series_stream {
24 static cl_pqa_series_term computenext (cl_pqa_series_stream& thisss)
26 var rational_series_stream& thiss = (rational_series_stream&)thisss;
27 var uintC n = thiss.n;
28 var cl_pqa_series_term result;
32 result.p = -expt_pos(n,5);
34 result.q = expt_pos(2*n+1,5)<<5;
35 result.a = 205*square((cl_I)n) + 250*(cl_I)n + 77;
39 rational_series_stream ()
40 : cl_pqa_series_stream (rational_series_stream::computenext),
45 // | ----- (n + 1) 2 |
46 // 1 | \ (-1) (205 n - 160 n + 32)|
47 // - | ) ---------------------------------|
49 // | ----- n binomial(2 n, n) |
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.
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);
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)))