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 uintL actuallen = len + 2; // 2 Schutz-Digits
30 var sintL 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 var uintL actuallen = len + 2; // 2 Schutz-Digits
54 // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
55 // with appropriate N, and
56 // a(n) = 1, b(n) = 2*n+1,
57 // p(n) = n for n>0, q(n) = 2*(2*n+1) for n>0.
58 var uintL N = (intDsize/2)*actuallen;
59 // 4^-N <= 2^(-intDsize*actuallen).
61 var cl_I* bv = (cl_I*) cl_alloca(N*sizeof(cl_I));
62 var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
63 var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
65 init1(cl_I, bv[0]) (1);
66 init1(cl_I, pv[0]) (1);
67 init1(cl_I, qv[0]) (1);
68 for (n = 1; n < N; n++) {
69 init1(cl_I, bv[n]) (2*n+1);
70 init1(cl_I, pv[n]) (n);
71 init1(cl_I, qv[n]) (2*(2*n+1));
73 var cl_pqb_series series;
75 series.pv = pv; series.qv = qv; series.qsv = NULL;
76 var cl_LF fsum = eval_rational_series(N,series,actuallen);
77 for (n = 0; n < N; n++) {
83 scale_float(The(cl_LF)(3*fsum)
84 + The(cl_LF)(pi(actuallen))
85 * ln(cl_I_to_LF(2,actuallen)+sqrt(cl_I_to_LF(3,actuallen))),
87 return shorten(g,len); // verkürzen und fertig
89 // Bit complexity (N := len): O(log(N)^2*M(N)).
91 const cl_LF compute_catalanconst_expintegral1 (uintC len)
93 // We compute f(x) classically and g(x) using the partial sums of f(x).
94 var uintC actuallen = len+2; // 2 Schutz-Digits
95 var uintL x = (uintL)(0.693148*intDsize*actuallen)+1;
96 var uintL N = (uintL)(2.718281828*x);
97 var cl_LF fterm = cl_I_to_LF(1,actuallen);
98 var cl_LF fsum = fterm;
99 var cl_LF gterm = fterm;
100 var cl_LF gsum = gterm;
103 // fterm = x^n/n!, fsum = 1 + x/1! + ... + x^n/n!,
104 // gterm = S_n*x^n/n!, gsum = S_0*x^0/0! + ... + S_n*x^n/n!.
105 for (n = 1; n < N; n++) {
106 fterm = The(cl_LF)(fterm*x)/n;
108 gterm = The(cl_LF)(gterm*x)/n;
110 gterm = gterm + fterm/square((cl_I)(2*n+1));
112 gterm = gterm - fterm/square((cl_I)(2*n+1));
115 var cl_LF result = gsum/fsum;
116 return shorten(result,len); // verkürzen und fertig
118 // Bit complexity (N = len): O(N^2).
120 // Same algorithm as expintegral1, but using binary splitting to evaluate
122 const cl_LF compute_catalanconst_expintegral2 (uintC len)
124 var uintC actuallen = len+2; // 2 Schutz-Digits
125 var uintL x = (uintL)(0.693148*intDsize*actuallen)+1;
126 var uintL N = (uintL)(2.718281828*x);
128 var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
130 for (n = 0; n < N; n++) {
132 init1(cl_I, args[n].p) (1);
133 init1(cl_I, args[n].q) (1);
135 init1(cl_I, args[n].p) (x);
136 init1(cl_I, args[n].q) (n);
138 init1(cl_I, args[n].d) (evenp(n)
139 ? square((cl_I)(2*n+1))
140 : -square((cl_I)(2*n+1)));
142 var cl_LF result = eval_pqd_series(N,args,actuallen);
143 for (n = 0; n < N; n++) {
148 return shorten(result,len); // verkürzen und fertig
150 // Bit complexity (N = len): O(log(N)^2*M(N)).
152 // Using Cohen-Villegas-Zagier acceleration, but without binary splitting.
153 const cl_LF compute_catalanconst_cvz1 (uintC len)
155 var uintC actuallen = len+2; // 2 Schutz-Digits
156 var uintL N = (uintL)(0.39321985*intDsize*actuallen)+1;
158 var cl_LF fterm = cl_I_to_LF(2*(cl_I)N*(cl_I)N,actuallen);
159 var cl_LF fsum = fterm;
160 var cl_LF gterm = fterm;
161 var cl_LF gsum = gterm;
164 // fterm = (N+n)!N/(2n+2)!(N-n-1)!*2^(2n+2), fsum = ... + fterm,
165 // gterm = S_n*fterm, gsum = ... + gterm.
166 for (n = 1; n < N; n++) {
167 fterm = The(cl_LF)(fterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
169 gterm = The(cl_LF)(gterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
171 gterm = gterm + fterm/square((cl_I)(2*n+1));
173 gterm = gterm - fterm/square((cl_I)(2*n+1));
176 var cl_LF result = gsum/(cl_I_to_LF(1,actuallen)+fsum);
178 // Take advantage of the fact that fterm and fsum are integers.
179 var cl_I fterm = 2*(cl_I)N*(cl_I)N;
180 var cl_I fsum = fterm;
181 var cl_LF gterm = cl_I_to_LF(fterm,actuallen);
182 var cl_LF gsum = gterm;
185 // fterm = (N+n)!N/(2n+2)!(N-n-1)!*2^(2n+2), fsum = ... + fterm,
186 // gterm = S_n*fterm, gsum = ... + gterm.
187 for (n = 1; n < N; n++) {
188 fterm = exquopos(fterm*(2*(cl_I)(N-n)*(cl_I)(N+n)),(cl_I)(2*n+1)*(cl_I)(n+1));
190 gterm = The(cl_LF)(gterm*(2*(cl_I)(N-n)*(cl_I)(N+n)))/((cl_I)(2*n+1)*(cl_I)(n+1));
192 gterm = gterm + cl_I_to_LF(fterm,actuallen)/square((cl_I)(2*n+1));
194 gterm = gterm - cl_I_to_LF(fterm,actuallen)/square((cl_I)(2*n+1));
197 var cl_LF result = gsum/cl_I_to_LF(1+fsum,actuallen);
199 return shorten(result,len); // verkürzen und fertig
201 // Bit complexity (N = len): O(N^2).
203 // Using Cohen-Villegas-Zagier acceleration, with binary splitting.
204 const cl_LF compute_catalanconst_cvz2 (uintC len)
206 var uintC actuallen = len+2; // 2 Schutz-Digits
207 var uintL N = (uintL)(0.39321985*intDsize*actuallen)+1;
209 var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
211 for (n = 0; n < N; n++) {
212 init1(cl_I, args[n].p) (2*(cl_I)(N-n)*(cl_I)(N+n));
213 init1(cl_I, args[n].q) ((cl_I)(2*n+1)*(cl_I)(n+1));
214 init1(cl_I, args[n].d) (evenp(n)
215 ? square((cl_I)(2*n+1))
216 : -square((cl_I)(2*n+1)));
218 var cl_pqd_series_result sums;
219 eval_pqd_series_aux(N,args,sums);
220 // Here we need U/(1+S) = V/D(Q+T).
222 cl_I_to_LF(sums.V,actuallen) / The(cl_LF)(sums.D * cl_I_to_LF(sums.Q+sums.T,actuallen));
223 for (n = 0; n < N; n++) {
228 return shorten(result,len); // verkürzen und fertig
230 // Bit complexity (N = len): O(log(N)^2*M(N)).
232 // Timings of the above algorithms, on an i486 33 MHz, running Linux.
233 // N ram ramfast exp1 exp2 cvz1 cvz2
234 // 10 0.055 0.068 0.32 0.91 0.076 0.11
235 // 25 0.17 0.26 0.95 3.78 0.23 0.43
236 // 50 0.43 0.73 2.81 11.5 0.63 1.36
237 // 100 1.32 2.24 8.82 34.1 1.90 4.48
238 // 250 6.60 10.4 48.7 127.5 10.3 20.8
239 // 500 24.0 30.9 186 329 38.4 58.6
240 // 1000 83.0 89.0 944 860 149 163
241 // 2500 446 352 6964 3096 1032 545
243 // asymp. N^2 FAST N^2 FAST N^2 FAST
244 // (FAST means O(log(N)^2*M(N)))
246 // The "exp1" and "exp2" algorithms are always about 10 to 15 times slower
247 // than the "ram" and "ramfast" algorithms.
248 // The break-even point between "ram" and "ramfast" is at about N = 1410.
250 const cl_LF compute_catalanconst (uintC len)
253 return compute_catalanconst_ramanujan_fast(len);
255 return compute_catalanconst_ramanujan(len);
257 // Bit complexity (N := len): O(log(N)^2*M(N)).
259 const cl_LF catalanconst (uintC len)
261 var uintC oldlen = TheLfloat(cl_LF_catalanconst)->len; // vorhandene Länge
263 return shorten(cl_LF_catalanconst,len);
265 return cl_LF_catalanconst;
267 // TheLfloat(cl_LF_catalanconst)->len um mindestens einen konstanten Faktor
268 // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
269 var uintC newlen = len;
270 oldlen += floor(oldlen,2); // oldlen * 3/2
274 // gewünschte > vorhandene Länge -> muß nachberechnen:
275 cl_LF_catalanconst = compute_catalanconst(newlen);
276 return (len < newlen ? shorten(cl_LF_catalanconst,len) : cl_LF_catalanconst);