12 #include "cln/lfloat.h"
13 #include "cl_LF_tran.h"
15 #include "cln/integer.h"
16 #include "cl_alloca.h"
20 #define floor cln_floor
24 const cl_LF compute_exp1 (uintC len)
26 // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
27 // with appropriate N, and
28 // a(n) = 1, b(n) = 1, p(n) = 1, q(n) = n for n>0.
29 var uintC actuallen = len+1; // 1 Schutz-Digit
30 // How many terms to we need for M bits of precision? N terms suffice,
33 // <== N*(log(N)-1) > M*log(2)
34 // First approximation:
35 // N0 = M will suffice, so put N<=N0.
36 // Second approximation:
37 // N1 = floor(M*log(2)/(log(N0)-1)), slightly too small, so put N>=N1.
38 // Third approximation:
39 // N2 = ceiling(M*log(2)/(log(N1)-1)), slightly too large.
40 // N = N2+2, two more terms for safety.
41 var uintL N0 = intDsize*actuallen;
42 var uintL N1 = (uintL)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0));
43 var uintL N2 = (uintL)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0))+1;
46 var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
47 var uintL* qsv = (uintL*) cl_alloca(N*sizeof(uintL));
49 for (n = 0; n < N; n++) {
50 init1(cl_I, qv[n]) (n==0 ? 1 : n);
52 var cl_q_series series;
54 series.qsv = (len >= 1000 && len <= 10000 ? qsv : 0); // tiny speedup
55 var cl_LF fsum = eval_rational_series(N,series,actuallen);
56 for (n = 0; n < N; n++) {
59 return shorten(fsum,len); // verkürzen und fertig
61 // Bit complexity (N = len): O(log(N)*M(N)).
63 // Timings of the above algorithm, on an i486 33 MHz, running Linux.
78 const cl_LF exp1 (uintC len)
80 var uintC oldlen = TheLfloat(cl_LF_exp1)->len; // vorhandene Länge
82 return shorten(cl_LF_exp1,len);
86 // TheLfloat(cl_LF_exp1)->len um mindestens einen konstanten Faktor
87 // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
88 var uintC newlen = len;
89 oldlen += floor(oldlen,2); // oldlen * 3/2
93 // gewünschte > vorhandene Länge -> muß nachberechnen:
94 cl_LF_exp1 = compute_exp1(newlen); // (exp 1)
95 return (len < newlen ? shorten(cl_LF_exp1,len) : cl_LF_exp1);