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"
19 #define floor cln_floor
23 const cl_LF compute_exp1 (uintC len)
25 // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
26 // with appropriate N, and
27 // a(n) = 1, b(n) = 1, p(n) = 1, q(n) = n for n>0.
28 var uintC actuallen = len+1; // 1 guard digit
29 // How many terms do we need for M bits of precision? N terms suffice,
32 // <== N*(log(N)-1) > M*log(2)
33 // First approximation:
34 // N0 = M will suffice, so put N<=N0.
35 // Second approximation:
36 // N1 = floor(M*log(2)/(log(N0)-1)), slightly too small, so put N>=N1.
37 // Third approximation:
38 // N2 = ceiling(M*log(2)/(log(N1)-1)), slightly too large.
39 // N = N2+2, two more terms for safety.
40 var uintC N0 = intDsize*actuallen;
41 var uintC N1 = (uintC)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0));
42 var uintC N2 = (uintC)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0))+1;
44 struct rational_series_stream : cl_q_series_stream {
46 static cl_q_series_term computenext (cl_q_series_stream& thisss)
48 var rational_series_stream& thiss = (rational_series_stream&)thisss;
49 var uintC n = thiss.n;
50 var cl_q_series_term result;
51 result.q = (n==0 ? 1 : n);
55 rational_series_stream()
56 : cl_q_series_stream (rational_series_stream::computenext),
59 var cl_LF fsum = eval_rational_series<false>(N,series,actuallen);
60 return shorten(fsum,len); // verkürzen und fertig
62 // Bit complexity (N = len): O(log(N)*M(N)).
64 // Timings of the above algorithm, on an i486 33 MHz, running Linux.
79 const cl_LF exp1 (uintC len)
81 var uintC oldlen = TheLfloat(cl_LF_exp1())->len; // vorhandene Länge
83 return shorten(cl_LF_exp1(),len);
87 // TheLfloat(cl_LF_exp1())->len um mindestens einen konstanten Faktor
88 // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
89 var uintC newlen = len;
90 oldlen += floor(oldlen,2); // oldlen * 3/2
94 // gewünschte > vorhandene Länge -> muß nachberechnen:
95 cl_LF_exp1() = compute_exp1(newlen); // (exp 1)
96 return (len < newlen ? shorten(cl_LF_exp1(),len) : cl_LF_exp1());