// exp1().
// General includes.
-#include "cl_sysdep.h"
+#include "base/cl_sysdep.h"
// Specification.
-#include "cl_F_tran.h"
+#include "float/transcendental/cl_F_tran.h"
// Implementation.
#include "cln/lfloat.h"
-#include "cl_LF_tran.h"
-#include "cl_LF.h"
+#include "float/transcendental/cl_LF_tran.h"
+#include "float/lfloat/cl_LF.h"
#include "cln/integer.h"
-#include "cl_alloca.h"
#undef floor
#include <cmath>
// Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
// with appropriate N, and
// a(n) = 1, b(n) = 1, p(n) = 1, q(n) = n for n>0.
- var uintC actuallen = len+1; // 1 Schutz-Digit
- // How many terms to we need for M bits of precision? N terms suffice,
+ var uintC actuallen = len+1; // 1 guard digit
+ // How many terms do we need for M bits of precision? N terms suffice,
// provided that
// 1/N! < 2^-M
// <== N*(log(N)-1) > M*log(2)
var uintC N1 = (uintC)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0));
var uintC N2 = (uintC)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0))+1;
var uintC N = N2+2;
- CL_ALLOCA_STACK;
- var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
- var uintC n;
- for (n = 0; n < N; n++) {
- init1(cl_I, qv[n]) (n==0 ? 1 : n);
- }
- var cl_q_series series;
- series.qv = qv;
- series.qsv = (len >= 1000 && len <= 10000 ? qsv : 0); // tiny speedup
- var cl_LF fsum = eval_rational_series(N,series,actuallen);
- for (n = 0; n < N; n++) {
- qv[n].~cl_I();
- }
+ struct rational_series_stream : cl_q_series_stream {
+ var uintC n;
+ static cl_q_series_term computenext (cl_q_series_stream& thisss)
+ {
+ var rational_series_stream& thiss = (rational_series_stream&)thisss;
+ var uintC n = thiss.n;
+ var cl_q_series_term result;
+ result.q = (n==0 ? 1 : n);
+ thiss.n = n+1;
+ return result;
+ }
+ rational_series_stream()
+ : cl_q_series_stream (rational_series_stream::computenext),
+ n(0) {}
+ } series;
+ var cl_LF fsum = eval_rational_series<false>(N,series,actuallen);
return shorten(fsum,len); // verkürzen und fertig
}
// Bit complexity (N = len): O(log(N)*M(N)).
const cl_LF exp1 (uintC len)
{
- var uintC oldlen = TheLfloat(cl_LF_exp1)->len; // vorhandene Länge
+ var uintC oldlen = TheLfloat(cl_LF_exp1())->len; // vorhandene Länge
if (len < oldlen)
- return shorten(cl_LF_exp1,len);
+ return shorten(cl_LF_exp1(),len);
if (len == oldlen)
- return cl_LF_exp1;
+ return cl_LF_exp1();
- // TheLfloat(cl_LF_exp1)->len um mindestens einen konstanten Faktor
+ // TheLfloat(cl_LF_exp1())->len um mindestens einen konstanten Faktor
// > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
var uintC newlen = len;
oldlen += floor(oldlen,2); // oldlen * 3/2
newlen = oldlen;
// gewünschte > vorhandene Länge -> muß nachberechnen:
- cl_LF_exp1 = compute_exp1(newlen); // (exp 1)
- return (len < newlen ? shorten(cl_LF_exp1,len) : cl_LF_exp1);
+ cl_LF_exp1() = compute_exp1(newlen); // (exp 1)
+ return (len < newlen ? shorten(cl_LF_exp1(),len) : cl_LF_exp1());
}
} // namespace cln