12 #include "cln/lfloat.h"
13 #include "cl_LF_tran.h"
15 #include "cln/integer.h"
16 #include "cln/exception.h"
20 #define floor cln_floor
24 const cl_LF cl_exp_aux (const cl_I& p, uintE lq, uintC len)
27 var uintE lp = integer_length(p); // now |p| < 2^lp.
28 if (!(lp <= lq)) throw runtime_exception();
29 lp = lq - lp; // now |p/2^lq| < 2^-lp.
30 // Minimize lq (saves computation time).
32 var uintC lp2 = ord2(p);
38 // Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
39 // with appropriate N, and
40 // a(n) = 1, b(n) = 1, p(n) = p for n>0, q(n) = n*2^lq for n>0.
41 var uintC actuallen = len+1; // 1 guard digit
42 // How many terms do we need for M bits of precision? N terms suffice,
44 // 1/(2^(N*lp)*N!) < 2^-M
45 // <== N*(log(N)-1)+N*lp*log(2) > M*log(2)
46 // First approximation:
47 // N0 = M will suffice, so put N<=N0.
48 // Second approximation:
49 // N1 = floor(M*log(2)/(log(N0)-1+lp*log(2))), slightly too small,
51 // Third approximation:
52 // N2 = ceiling(M*log(2)/(log(N1)-1+lp*log(2))), slightly too large.
53 // N = N2+2, two more terms for safety.
54 var uintC N0 = intDsize*actuallen;
55 var uintC N1 = (uintC)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0+0.693148*lp));
56 var uintC N2 = (uintC)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0+0.693147*lp))+1;
58 struct rational_series_stream : cl_pq_series_stream {
62 static cl_pq_series_term computenext (cl_pq_series_stream& thisss)
64 var rational_series_stream& thiss = (rational_series_stream&)thisss;
65 var uintC n = thiss.n;
66 var cl_pq_series_term result;
72 result.q = (cl_I)n << thiss.lq;
77 rational_series_stream(const cl_I& p_, uintE lq_)
78 : cl_pq_series_stream (rational_series_stream::computenext),
79 n (0), p(p_), lq(lq_) {}
81 var cl_LF fsum = eval_rational_series<true>(N,series,actuallen);
82 return shorten(fsum,len); // verkürzen und fertig
84 // Bit complexity (N = len, and if p has length O(log N) and ql = O(log N)):