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 compute_catalanconst_ramanujan (uintC len)
22 // [Jonathan M. Borwein, Peter B. Borwein: Pi and the AGM.
23 // Wiley 1987. Section 11.3, exercise 16 g, p. 386]
24 // G = 3/8 * sum(n=0..infty, n!^2 / (2n+1)!*(2n+1))
25 // + pi/8 * log(2+sqrt(3)).
26 // Every summand gives 0.6 new decimal digits in precision.
27 // The sum is best evaluated using fixed-point arithmetic,
28 // so that the precision is reduced for the later summands.
29 var uintC actuallen = len + 2; // 2 guard digits
30 var sintC scale = intDsize*actuallen;
33 var cl_I factor = ash(1,scale);
34 while (!zerop(factor)) {
35 sum = sum + truncate1(factor,2*n+1);
37 factor = truncate1(factor*n,2*(2*n+1));
39 var cl_LF fsum = scale_float(cl_I_to_LF(sum,actuallen),-scale);
41 scale_float(The(cl_LF)(3*fsum)
42 + The(cl_LF)(pi(actuallen))
43 * ln(cl_I_to_LF(2,actuallen)+sqrt(cl_I_to_LF(3,actuallen))),
45 return shorten(g,len); // verkürzen und fertig
47 // Bit complexity (N := len): O(N^2).
49 const cl_LF compute_catalanconst_ramanujan_fast (uintC len)
51 // Same formula as above, using a binary splitting evaluation.
52 // See [Borwein, Borwein, section 10.2.3].
53 struct rational_series_stream : cl_pqb_series_stream {
55 static cl_pqb_series_term computenext (cl_pqb_series_stream& thisss)
57 var rational_series_stream& thiss = (rational_series_stream&)thisss;
59 var cl_pqb_series_term result;
67 result.q = result.b << 1; // 2*(2*n+1)
72 rational_series_stream ()
73 : cl_pqb_series_stream (rational_series_stream::computenext),
76 var uintC actuallen = len + 2; // 2 guard digits
77 // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
78 // with appropriate N, and
79 // a(n) = 1, b(n) = 2*n+1,
80 // p(n) = n for n>0, q(n) = 2*(2*n+1) for n>0.
81 var uintC N = (intDsize/2)*actuallen;
82 // 4^-N <= 2^(-intDsize*actuallen).
83 var cl_LF fsum = eval_rational_series<false>(N,series,actuallen,actuallen);
85 scale_float(The(cl_LF)(3*fsum)
86 + The(cl_LF)(pi(actuallen))
87 * ln(cl_I_to_LF(2,actuallen)+sqrt(cl_I_to_LF(3,actuallen))),
89 return shorten(g,len); // verkürzen und fertig
91 // Bit complexity (N := len): O(log(N)^2*M(N)).
93 const cl_LF compute_catalanconst_expintegral1 (uintC len)
95 // We compute f(x) classically and g(x) using the partial sums of f(x).
96 var uintC actuallen = len+2; // 2 guard digits
97 var uintC x = (uintC)(0.693148*intDsize*actuallen)+1;
98 var uintC N = (uintC)(2.718281828*x);
99 var cl_LF fterm = cl_I_to_LF(1,actuallen);
100 var cl_LF fsum = fterm;
101 var cl_LF gterm = fterm;
102 var cl_LF gsum = gterm;
105 // fterm = x^n/n!, fsum = 1 + x/1! + ... + x^n/n!,
106 // gterm = S_n*x^n/n!, gsum = S_0*x^0/0! + ... + S_n*x^n/n!.
107 for (n = 1; n < N; n++) {
108 fterm = The(cl_LF)(fterm*x)/n;
110 gterm = The(cl_LF)(gterm*x)/n;
112 gterm = gterm + fterm/square((cl_I)(2*n+1));
114 gterm = gterm - fterm/square((cl_I)(2*n+1));
117 var cl_LF result = gsum/fsum;
118 return shorten(result,len); // verkürzen und fertig
120 // Bit complexity (N = len): O(N^2).
122 // Same algorithm as expintegral1, but using binary splitting to evaluate
124 const cl_LF compute_catalanconst_expintegral2 (uintC len)
126 var uintC actuallen = len+2; // 2 guard digits
127 var uintC x = (uintC)(0.693148*intDsize*actuallen)+1;
128 var uintC N = (uintC)(2.718281828*x);
130 var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
132 for (n = 0; n < N; n++) {
134 init1(cl_I, args[n].p) (1);
135 init1(cl_I, args[n].q) (1);
137 init1(cl_I, args[n].p) (x);
138 init1(cl_I, args[n].q) (n);
140 init1(cl_I, args[n].d) (evenp(n)
141 ? square((cl_I)(2*n+1))
142 : -square((cl_I)(2*n+1)));
144 var cl_LF result = eval_pqd_series(N,args,actuallen);
145 for (n = 0; n < N; n++) {
150 return shorten(result,len); // verkürzen und fertig
152 // Bit complexity (N = len): O(log(N)^2*M(N)).
154 // Using Cohen-Villegas-Zagier acceleration, but without binary splitting.
155 const cl_LF compute_catalanconst_cvz1 (uintC len)
157 var uintC actuallen = len+2; // 2 guard digits
158 var uintC N = (uintC)(0.39321985*intDsize*actuallen)+1;
160 var cl_LF fterm = cl_I_to_LF(2*(cl_I)N*(cl_I)N,actuallen);
161 var cl_LF fsum = fterm;
162 var cl_LF gterm = fterm;
163 var cl_LF gsum = gterm;
166 // fterm = (N+n)!N/(2n+2)!(N-n-1)!*2^(2n+2), fsum = ... + fterm,
167 // gterm = S_n*fterm, gsum = ... + gterm.
168 for (n = 1; n < N; n++) {
169 fterm = The(cl_LF)(fterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
171 gterm = The(cl_LF)(gterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
173 gterm = gterm + fterm/square((cl_I)(2*n+1));
175 gterm = gterm - fterm/square((cl_I)(2*n+1));
178 var cl_LF result = gsum/(cl_I_to_LF(1,actuallen)+fsum);
180 // Take advantage of the fact that fterm and fsum are integers.
181 var cl_I fterm = 2*(cl_I)N*(cl_I)N;
182 var cl_I fsum = fterm;
183 var cl_LF gterm = cl_I_to_LF(fterm,actuallen);
184 var cl_LF gsum = gterm;
187 // fterm = (N+n)!N/(2n+2)!(N-n-1)!*2^(2n+2), fsum = ... + fterm,
188 // gterm = S_n*fterm, gsum = ... + gterm.
189 for (n = 1; n < N; n++) {
190 fterm = exquopos(fterm*(2*(cl_I)(N-n)*(cl_I)(N+n)),(cl_I)(2*n+1)*(cl_I)(n+1));
192 gterm = The(cl_LF)(gterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
194 gterm = gterm + cl_I_to_LF(fterm,actuallen)/square((cl_I)(2*n+1));
196 gterm = gterm - cl_I_to_LF(fterm,actuallen)/square((cl_I)(2*n+1));
199 var cl_LF result = gsum/cl_I_to_LF(1+fsum,actuallen);
201 return shorten(result,len); // verkürzen und fertig
203 // Bit complexity (N = len): O(N^2).
205 // Using Cohen-Villegas-Zagier acceleration, with binary splitting.
206 const cl_LF compute_catalanconst_cvz2 (uintC len)
208 var uintC actuallen = len+2; // 2 guard digits
209 var uintC N = (uintC)(0.39321985*intDsize*actuallen)+1;
211 var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
213 for (n = 0; n < N; n++) {
214 init1(cl_I, args[n].p) (2*(cl_I)(N-n)*(cl_I)(N+n));
215 init1(cl_I, args[n].q) ((cl_I)(2*n+1)*(cl_I)(n+1));
216 init1(cl_I, args[n].d) (evenp(n)
217 ? square((cl_I)(2*n+1))
218 : -square((cl_I)(2*n+1)));
220 var cl_pqd_series_result<cl_I> sums;
221 eval_pqd_series_aux(N,args,sums);
222 // Here we need U/(1+S) = V/D(Q+T).
224 cl_I_to_LF(sums.V,actuallen) / The(cl_LF)(sums.D * cl_I_to_LF(sums.Q+sums.T,actuallen));
225 for (n = 0; n < N; n++) {
230 return shorten(result,len); // verkürzen und fertig
232 // Bit complexity (N = len): O(log(N)^2*M(N)).
235 const cl_LF compute_catalanconst_lupas (uintC len)
237 // [Alexandru Lupas. Formulae for Some Classical Constants.
238 // Proceedings of ROGER-2000.
239 // http://www.lacim.uqam.ca/~plouffe/articles/alupas1.pdf]
240 // [http://mathworld.wolfram.com/CatalansConstant.html]
241 // G = 19/18 * sum(n=0..infty,
242 // mul(m=1..n, -32*((80*m^3+72*m^2-18*m-19)*m^3)/
243 // (10240*m^6+14336*m^5+2560*m^4-3072*m^3-888*m^2+72*m+27))).
244 struct rational_series_stream : cl_pq_series_stream {
246 static cl_pq_series_term computenext (cl_pq_series_stream& thisss)
248 var rational_series_stream& thiss = (rational_series_stream&)thisss;
249 var cl_I n = thiss.n;
250 var cl_pq_series_term result;
255 // Compute -32*((80*n^3+72*n^2-18*n-19)*n^3) (using Horner scheme):
256 result.p = (19+(18+(-72-80*n)*n)*n)*n*n*n << 5;
257 // Compute 10240*n^6+14336*n^5+2560*n^4-3072*n^3-888*n^2+72*n+27:
258 result.q = 27+(72+(-888+(-3072+(2560+(14336+10240*n)*n)*n)*n)*n)*n;
263 rational_series_stream ()
264 : cl_pq_series_stream (rational_series_stream::computenext),
267 var uintC actuallen = len + 2; // 2 guard digits
268 var uintC N = (intDsize/2)*actuallen;
269 var cl_LF fsum = eval_rational_series<false>(N,series,actuallen,actuallen);
270 var cl_LF g = fsum*cl_I_to_LF(19,actuallen)/cl_I_to_LF(18,actuallen);
271 return shorten(g,len);
273 // Bit complexity (N = len): O(log(N)^2*M(N)).
275 // Timings of the above algorithms, on an AMD Opteron 1.7 GHz, running Linux/x86.
276 // N ram ramfast exp1 exp2 cvz1 cvz2 lupas
277 // 25 0.0011 0.0010 0.0094 0.0107 0.0021 0.0016 0.0034
278 // 50 0.0030 0.0025 0.0280 0.0317 0.0058 0.0045 0.0095
279 // 100 0.0087 0.0067 0.0854 0.0941 0.0176 0.0121 0.0260
280 // 250 0.043 0.029 0.462 0.393 0.088 0.055 0.109
281 // 500 0.15 0.086 1.7 1.156 0.323 0.162 0.315
282 // 1000 0.57 0.25 7.5 3.23 1.27 0.487 0.864
283 // 2500 3.24 1.10 52.2 12.4 8.04 2.02 3.33
284 // 5000 13.1 3.06 218 32.7 42.1 5.59 8.80
285 // 10000 52.7 8.2 910 85.3 216.7 15.3 22.7
286 // 25000 342 29.7 6403 295 1643 54.5 77.3
287 // 50000 89.9 139 195
289 // asymp. N^2 FAST N^2 FAST N^2 FAST FAST
291 // (FAST means O(log(N)^2*M(N)))
293 // The "ramfast" algorithm is always fastest.
295 const cl_LF compute_catalanconst (uintC len)
297 return compute_catalanconst_ramanujan_fast(len);
299 // Bit complexity (N := len): O(log(N)^2*M(N)).
301 const cl_LF catalanconst (uintC len)
303 var uintC oldlen = TheLfloat(cl_LF_catalanconst)->len; // vorhandene Länge
305 return shorten(cl_LF_catalanconst,len);
307 return cl_LF_catalanconst;
309 // TheLfloat(cl_LF_catalanconst)->len um mindestens einen konstanten Faktor
310 // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
311 var uintC newlen = len;
312 oldlen += floor(oldlen,2); // oldlen * 3/2
316 // gewünschte > vorhandene Länge -> muß nachberechnen:
317 cl_LF_catalanconst = compute_catalanconst(newlen);
318 return (len < newlen ? shorten(cl_LF_catalanconst,len) : cl_LF_catalanconst);