12 #include "cln/lfloat.h"
13 #include "cl_LF_tran.h"
15 #include "cln/integer.h"
16 #include "cln/abort.h"
17 #include "cl_alloca.h"
21 const cl_LF compute_zeta_exp (int s, uintC len)
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 uintC x = (uintC)(0.693148*intDsize*actuallen)+1;
29 var uintC N = (uintC)(2.718281828*x);
31 var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
33 for (n = 0; n < N; n++) {
35 init1(cl_I, args[n].p) (1);
36 init1(cl_I, args[n].q) (1);
38 init1(cl_I, args[n].p) (x);
39 init1(cl_I, args[n].q) (n);
41 init1(cl_I, args[n].d) (evenp(n)
45 var cl_LF result = eval_pqd_series(N,args,actuallen);
46 for (n = 0; n < N; n++) {
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);
55 // Bit complexity (N = len): O(log(N)^2*M(N)).
57 const cl_LF compute_zeta_cvz1 (int s, uintC len)
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 uintC N = (uintC)(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;
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));
75 gterm = The(cl_LF)(gterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
77 gterm = gterm + cl_I_to_LF(fterm,actuallen)/expt_pos(n+1,s);
79 gterm = gterm - cl_I_to_LF(fterm,actuallen)/expt_pos(n+1,s);
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);
87 // Bit complexity (N = len): O(N^2).
89 const cl_LF compute_zeta_cvz2 (int s, uintC len)
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 uintC N = (uintC)(0.39321985*intDsize*actuallen)+1;
98 var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
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)
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).
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++) {
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);
121 // Bit complexity (N = len): O(log(N)^2*M(N)).
123 // Timings of the above algorithm, on an i486 33 MHz, running Linux.
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
135 // The break-even point between cvz1 and cvz2 seems to grow linearly with s.
137 const cl_LF zeta (int s, uintC len)
141 if (len < 280*(uintL)s)
142 return compute_zeta_cvz1(s,len);
144 return compute_zeta_cvz2(s,len);
146 // Bit complexity (N = len): O(log(N)^2*M(N)).