From 4109ef07d74b779bb4d4371c24a6364302f464bb Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Fri, 7 Sep 2007 19:45:30 +0000 Subject: [PATCH] More memory efficient Euler-Mascheroni constant: * src/float/transcendental/cl_LF_tran.h (cl_pqd_series_stream): New. * src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc: New file. * src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc: New file. * src/float/transcendental/cl_LF_eulerconst.cc: Compute series coefficients on demand, using a series stream object. --- ChangeLog | 10 +++ src/float/transcendental/cl_LF_eulerconst.cc | 33 ++++--- .../cl_LF_ratsumseries_stream_pqd.cc | 31 +++++++ .../cl_LF_ratsumseries_stream_pqd_aux.cc | 86 +++++++++++++++++++ src/float/transcendental/cl_LF_tran.h | 8 ++ 5 files changed, 154 insertions(+), 14 deletions(-) create mode 100644 src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc create mode 100644 src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc diff --git a/ChangeLog b/ChangeLog index aebc719..27391ec 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,13 @@ +2007-09-07 Richard B. Kreckel + + More memory efficient Euler-Mascheroni constant: + * src/float/transcendental/cl_LF_tran.h (cl_pqd_series_stream): New. + * src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc: New file. + * src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc: New + file. + * src/float/transcendental/cl_LF_eulerconst.cc: Compute series + coefficients on demand, using a series stream object. + 2007-08-02 Richard B. Kreckel * src/base/digitseq/cl_DS_div.cc (cl_recip_suitable): uintC arguments. diff --git a/src/float/transcendental/cl_LF_eulerconst.cc b/src/float/transcendental/cl_LF_eulerconst.cc index 1258f75..3e13c19 100644 --- a/src/float/transcendental/cl_LF_eulerconst.cc +++ b/src/float/transcendental/cl_LF_eulerconst.cc @@ -447,16 +447,26 @@ const cl_LF compute_eulerconst_besselintegral4 (uintC len) var uintC sx = (uintC)(0.25*0.693148*intDsize*actuallen)+1; var uintC N = (uintC)(3.591121477*sx); var cl_I x = square((cl_I)sx); - CL_ALLOCA_STACK; - var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term)); - var uintC n; - for (n = 0; n < N; n++) { - init1(cl_I, args[n].p) (x); - init1(cl_I, args[n].q) (square((cl_I)(n+1))); - init1(cl_I, args[n].d) (n+1); - } + struct rational_series_stream : cl_pqd_series_stream { + uintC n; + cl_I x; + static cl_pqd_series_term computenext (cl_pqd_series_stream& thisss) + { + var rational_series_stream& thiss = (rational_series_stream&)thisss; + var uintC n = thiss.n; + var cl_pqd_series_term result; + result.p = thiss.x; + result.q = square((cl_I)(n+1)); + result.d = n+1; + thiss.n = n+1; + return result; + } + rational_series_stream (const cl_I& _x) + : cl_pqd_series_stream (rational_series_stream::computenext), + n (0), x (_x) {} + } series(x); var cl_pqd_series_result sums; - eval_pqd_series_aux(N,args,sums); + eval_pqd_series_aux(N,series,sums); // Instead of computing fsum = 1 + T/Q and gsum = V/(D*Q) // and then dividing them, to compute gsum/fsum, we save two // divisions by computing V/(D*(Q+T)). @@ -464,11 +474,6 @@ const cl_LF compute_eulerconst_besselintegral4 (uintC len) cl_I_to_LF(sums.V,actuallen) / The(cl_LF)(sums.D * cl_I_to_LF(sums.Q+sums.T,actuallen)) - ln(cl_I_to_LF(sx,actuallen)); - for (n = 0; n < N; n++) { - args[n].p.~cl_I(); - args[n].q.~cl_I(); - args[n].d.~cl_I(); - } return shorten(result,len); // verkürzen und fertig } // Bit complexity (N = len): O(log(N)^2*M(N)). diff --git a/src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc b/src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc new file mode 100644 index 0000000..2f1aaaf --- /dev/null +++ b/src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc @@ -0,0 +1,31 @@ +// eval_pqd_series(). + +// General includes. +#include "cl_sysdep.h" + +// Specification. +#include "cl_LF_tran.h" + + +// Implementation. + +#include "cln/lfloat.h" +#include "cln/integer.h" +#include "cl_LF.h" + +namespace cln { + +const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_pqd_series_result sums; + eval_pqd_series_aux(N,args,sums); + // Instead of computing fsum = T/Q and gsum = V/(D*Q) + // and then dividing them, to compute gsum/fsum, we save two + // divisions by computing V/(D*T). + return + cl_I_to_LF(sums.V,len) / The(cl_LF)(sums.D * cl_I_to_LF(sums.T,len)); +} + +} // namespace cln diff --git a/src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc b/src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc new file mode 100644 index 0000000..cc64077 --- /dev/null +++ b/src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc @@ -0,0 +1,86 @@ +// eval_pqd_series_aux(). + +// General includes. +#include "cl_sysdep.h" + +// Specification. +#include "cl_LF_tran.h" + + +// Implementation. + +#include "cln/integer.h" +#include "cln/exception.h" + +namespace cln { + +void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result& Z, bool rightmost) +{ + // N = N2-N1 + switch (N) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pqd_series_term v0 = args.next(); // [N1] + if (!rightmost) { Z.P = v0.p; } + Z.Q = v0.q; + Z.T = v0.p; + if (!rightmost) { Z.C = 1; } + Z.D = v0.d; + Z.V = v0.p; + break; + } + case 2: { + var cl_pqd_series_term v0 = args.next(); // [N1] + var cl_pqd_series_term v1 = args.next(); // [N1+1] + var cl_I p01 = v0.p * v1.p; + if (!rightmost) { Z.P = p01; } + Z.Q = v0.q * v1.q; + var cl_I p0q1 = v0.p * v1.q + p01; + Z.T = p0q1; + if (!rightmost) { Z.C = v1.d + v0.d; } + Z.D = v0.d * v1.d; + Z.V = v1.d * p0q1 + v0.d * p01; + break; + } + case 3: { + var cl_pqd_series_term v0 = args.next(); // [N1] + var cl_pqd_series_term v1 = args.next(); // [N1+1] + var cl_pqd_series_term v2 = args.next(); // [N1+2] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + if (!rightmost) { Z.P = p012; } + Z.Q = v0.q * v1.q * v2.q; + var cl_I p0q1 = v0.p * v1.q + p01; + Z.T = v2.q * p0q1 + p012; + var cl_I d01 = v0.d * v1.d; + if (!rightmost) { Z.C = (v1.d + v0.d) * v2.d + d01; } + Z.D = d01 * v2.d; + Z.V = v2.d * (v2.q * (v1.d * p0q1 + v0.d * p01) + (v1.d + v0.d) * p012) + d01 * p012; + break; + } + default: { + var uintC Nm = N/2; // midpoint + // Compute left part. + var cl_pqd_series_result L; + eval_pqd_series_aux(Nm,args,L,false); + // Compute right part. + var cl_pqd_series_result R; + eval_pqd_series_aux(N-Nm,args,R,rightmost); + // Put together partial results. + if (!rightmost) { Z.P = L.P * R.P; } + Z.Q = L.Q * R.Q; + // Z.S = L.S + L.P/L.Q*R.S; + var cl_I tmp = L.P * R.T; + Z.T = R.Q * L.T + tmp; + if (!rightmost) { Z.C = L.C * R.D + L.D * R.C; } + Z.D = L.D * R.D; + // Z.U = L.U + L.C/L.D * L.P/L.Q * R.S + L.P/L.Q * R.U; + // Z.V = R.D * R.Q * L.V + R.D * L.C * L.P * R.T + L.D * L.P * R.V; + Z.V = R.D * (R.Q * L.V + L.C * tmp) + L.D * L.P * R.V; + break; + } + } +} + +} // namespace cln diff --git a/src/float/transcendental/cl_LF_tran.h b/src/float/transcendental/cl_LF_tran.h index 9a32a13..b975fa3 100644 --- a/src/float/transcendental/cl_LF_tran.h +++ b/src/float/transcendental/cl_LF_tran.h @@ -255,9 +255,17 @@ struct cl_pqd_series_result { cl_I D; cl_I V; }; +struct cl_pqd_series_stream { + cl_pqd_series_term (*nextfn)(cl_pqd_series_stream&); + cl_pqd_series_term next () { return nextfn(*this); } + // Constructor. + cl_pqd_series_stream( cl_pqd_series_term (*n)(cl_pqd_series_stream&)) : nextfn (n) {} +}; extern void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_result& Z, bool rightmost = true); +extern void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result& Z, bool rightmost = true); // Ditto, but returns U/S. extern const cl_LF eval_pqd_series (uintC N, cl_pqd_series_term* args, uintC len); +extern const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len); } // namespace cln