From 962f4b052d22895fcc0f3a7a09614bafe152ad71 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Thu, 13 Sep 2007 19:47:51 +0000 Subject: [PATCH] Truncated binary splitting for even more memory efficiency: * src/float/transcendental/cl_LF_tran.h: Added new overloads. See below. * src/float/transcendental/cl_LF_ratseries_stream_pq.cc: Removed and moved everything to... * src/float/transcendental/cl_LF_ratseries_pq.cc: ...here. Added an overload for truncated expansion. * src/float/transcendental/cl_LF_ratseries_stream_pqa.cc: Removed and moved everything to... * src/float/transcendental/cl_LF_ratseries_pqa.cc: ...here. Added an overload for truncated expansion. * src/float/transcendental/cl_LF_ratseries_stream_pqb.cc: Removed and moved everything to... * src/float/transcendental/cl_LF_ratseries_pqb.cc: ...here. Added an overload for truncated expansion. * src/float/transcendental/cl_LF_ratseries_stream_pqab.cc: Removed and moved everything to... * src/float/transcendental/cl_LF_ratseries_pqab.cc: ...here. Added an overload for truncated expansion. * src/float/transcendental/cl_LF_ratsumseries_pqcd_aux.cc: Added overloads for streamed and truncated expansion. * src/float/transcendental/cl_LF_ratsumseries_pqcd.cc: Likewise. * src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc: Removed and moved everything to... * src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc: ...here. Added an overload for truncated expansion. * src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc: Removed and moved everything to... * src/float/transcendental/cl_LF_ratsumseries_pqd.cc: ...here. Added an overload for truncated expansion. * src/float/transcendental/cl_LF_pi.cc: Use truncated series. * src/float/transcendental/cl_LF_catalanconst.cc: Likewise. * src/float/transcendental/cl_LF_eulerconst.cc: Likewise. * src/float/transcendental/cl_LF_zeta_int.cc: Likewise. * src/float/transcendental/cl_LF_zeta3.cc: Likewise. --- ChangeLog | 37 ++++ .../transcendental/cl_LF_catalanconst.cc | 6 +- src/float/transcendental/cl_LF_eulerconst.cc | 23 +- src/float/transcendental/cl_LF_pi.cc | 2 +- .../transcendental/cl_LF_ratseries_pq.cc | 173 +++++++++++++++ .../transcendental/cl_LF_ratseries_pqa.cc | 173 +++++++++++++++ .../transcendental/cl_LF_ratseries_pqab.cc | 190 +++++++++++++++++ .../transcendental/cl_LF_ratseries_pqb.cc | 191 +++++++++++++++++ .../cl_LF_ratseries_stream_pq.cc | 102 --------- .../cl_LF_ratseries_stream_pqa.cc | 102 --------- .../cl_LF_ratseries_stream_pqab.cc | 110 ---------- .../cl_LF_ratseries_stream_pqb.cc | 110 ---------- .../transcendental/cl_LF_ratsumseries_pqcd.cc | 29 ++- .../cl_LF_ratsumseries_pqcd_aux.cc | 163 +++++++++++++- .../transcendental/cl_LF_ratsumseries_pqd.cc | 29 ++- .../cl_LF_ratsumseries_pqd_aux.cc | 155 +++++++++++++- .../cl_LF_ratsumseries_stream_pqd.cc | 31 --- .../cl_LF_ratsumseries_stream_pqd_aux.cc | 86 -------- src/float/transcendental/cl_LF_tran.h | 199 +++++++++++------- src/float/transcendental/cl_LF_zeta3.cc | 2 +- src/float/transcendental/cl_LF_zeta_int.cc | 75 +++++-- 21 files changed, 1325 insertions(+), 663 deletions(-) delete mode 100644 src/float/transcendental/cl_LF_ratseries_stream_pq.cc delete mode 100644 src/float/transcendental/cl_LF_ratseries_stream_pqa.cc delete mode 100644 src/float/transcendental/cl_LF_ratseries_stream_pqab.cc delete mode 100644 src/float/transcendental/cl_LF_ratseries_stream_pqb.cc delete mode 100644 src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc delete mode 100644 src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc diff --git a/ChangeLog b/ChangeLog index 27391ec..3bf642f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,40 @@ +2007-09-13 Richard B. Kreckel + + Truncated binary splitting for even more memory efficiency: + * src/float/transcendental/cl_LF_tran.h: Added new overloads. See below. + * src/float/transcendental/cl_LF_ratseries_stream_pq.cc: Removed and + moved everything to... + * src/float/transcendental/cl_LF_ratseries_pq.cc: ...here. Added an + overload for truncated expansion. + * src/float/transcendental/cl_LF_ratseries_stream_pqa.cc: Removed and + moved everything to... + * src/float/transcendental/cl_LF_ratseries_pqa.cc: ...here. Added an + overload for truncated expansion. + * src/float/transcendental/cl_LF_ratseries_stream_pqb.cc: Removed and + moved everything to... + * src/float/transcendental/cl_LF_ratseries_pqb.cc: ...here. Added an + overload for truncated expansion. + * src/float/transcendental/cl_LF_ratseries_stream_pqab.cc: Removed and + moved everything to... + * src/float/transcendental/cl_LF_ratseries_pqab.cc: ...here. Added an + overload for truncated expansion. + * src/float/transcendental/cl_LF_ratsumseries_pqcd_aux.cc: Added + overloads for streamed and truncated expansion. + * src/float/transcendental/cl_LF_ratsumseries_pqcd.cc: Likewise. + * src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc: Removed + and moved everything to... + * src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc: ...here. Added + an overload for truncated expansion. + * src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc: Removed + and moved everything to... + * src/float/transcendental/cl_LF_ratsumseries_pqd.cc: ...here. Added an + overload for truncated expansion. + * src/float/transcendental/cl_LF_pi.cc: Use truncated series. + * src/float/transcendental/cl_LF_catalanconst.cc: Likewise. + * src/float/transcendental/cl_LF_eulerconst.cc: Likewise. + * src/float/transcendental/cl_LF_zeta_int.cc: Likewise. + * src/float/transcendental/cl_LF_zeta3.cc: Likewise. + 2007-09-07 Richard B. Kreckel More memory efficient Euler-Mascheroni constant: diff --git a/src/float/transcendental/cl_LF_catalanconst.cc b/src/float/transcendental/cl_LF_catalanconst.cc index 28fbd2b..e2acd23 100644 --- a/src/float/transcendental/cl_LF_catalanconst.cc +++ b/src/float/transcendental/cl_LF_catalanconst.cc @@ -80,7 +80,7 @@ const cl_LF compute_catalanconst_ramanujan_fast (uintC len) // p(n) = n for n>0, q(n) = 2*(2*n+1) for n>0. var uintC N = (intDsize/2)*actuallen; // 4^-N <= 2^(-intDsize*actuallen). - var cl_LF fsum = eval_rational_series(N,series,actuallen); + var cl_LF fsum = eval_rational_series(N,series,actuallen,actuallen); var cl_LF g = scale_float(The(cl_LF)(3*fsum) + The(cl_LF)(pi(actuallen)) @@ -217,7 +217,7 @@ const cl_LF compute_catalanconst_cvz2 (uintC len) ? square((cl_I)(2*n+1)) : -square((cl_I)(2*n+1))); } - var cl_pqd_series_result sums; + var cl_pqd_series_result sums; eval_pqd_series_aux(N,args,sums); // Here we need U/(1+S) = V/D(Q+T). var cl_LF result = @@ -266,7 +266,7 @@ const cl_LF compute_catalanconst_lupas (uintC len) } series; var uintC actuallen = len + 2; // 2 guard digits var uintC N = (intDsize/2)*actuallen; - var cl_LF fsum = eval_rational_series(N,series,actuallen); + var cl_LF fsum = eval_rational_series(N,series,actuallen,actuallen); var cl_LF g = fsum*cl_I_to_LF(19,actuallen)/cl_I_to_LF(18,actuallen); return shorten(g,len); } diff --git a/src/float/transcendental/cl_LF_eulerconst.cc b/src/float/transcendental/cl_LF_eulerconst.cc index 3e13c19..467dfb6 100644 --- a/src/float/transcendental/cl_LF_eulerconst.cc +++ b/src/float/transcendental/cl_LF_eulerconst.cc @@ -13,6 +13,7 @@ #include "cl_LF_tran.h" #include "cl_LF.h" #include "cln/integer.h" +#include "cln/real.h" #include "cl_alloca.h" namespace cln { @@ -443,10 +444,10 @@ const cl_LF compute_eulerconst_besselintegral3 (uintC len) // the sums. const cl_LF compute_eulerconst_besselintegral4 (uintC len) { - var uintC actuallen = len+2; // 2 Schutz-Digits + var uintC actuallen = len+2; // 2 guard digits 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); + var cl_I x = square(cl_I(sx)); struct rational_series_stream : cl_pqd_series_stream { uintC n; cl_I x; @@ -456,24 +457,24 @@ const cl_LF compute_eulerconst_besselintegral4 (uintC len) var uintC n = thiss.n; var cl_pqd_series_term result; result.p = thiss.x; - result.q = square((cl_I)(n+1)); + result.q = square(cl_I(n+1)); result.d = n+1; thiss.n = n+1; return result; } - rational_series_stream (const cl_I& _x) + rational_series_stream (uintC _n, 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,series,sums); + n (_n), x (_x) {} + } series(0,x); + var cl_pqd_series_result sums; + eval_pqd_series_aux(N,series,sums,actuallen); // 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)). var cl_LF result = - 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)); + cl_R_to_LF(sums.V,actuallen) + / The(cl_LF)(sums.D * cl_R_to_LF(sums.Q+sums.T,actuallen)) + - ln(cl_R_to_LF(sx,actuallen)); 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_pi.cc b/src/float/transcendental/cl_LF_pi.cc index 6660f62..6b0a263 100644 --- a/src/float/transcendental/cl_LF_pi.cc +++ b/src/float/transcendental/cl_LF_pi.cc @@ -230,7 +230,7 @@ const cl_LF compute_pi_ramanujan_163_fast (uintC len) var uintC N = (n_slope*actuallen)/32 + 1; // N > intDsize*log(2)/log(|J|) * actuallen, hence // |J|^-N < 2^(-intDsize*actuallen). - var cl_LF fsum = eval_rational_series(N,series,actuallen); + var cl_LF fsum = eval_rational_series(N,series,actuallen,actuallen); static const cl_I J3 = "262537412640768000"; // -1728*J var cl_LF pires = sqrt(cl_I_to_LF(J3,actuallen)) / fsum; return shorten(pires,len); // verkürzen und fertig diff --git a/src/float/transcendental/cl_LF_ratseries_pq.cc b/src/float/transcendental/cl_LF_ratseries_pq.cc index 4bc195d..2939900 100644 --- a/src/float/transcendental/cl_LF_ratseries_pq.cc +++ b/src/float/transcendental/cl_LF_ratseries_pq.cc @@ -11,6 +11,7 @@ #include "cln/lfloat.h" #include "cln/integer.h" +#include "cln/real.h" #include "cln/exception.h" #include "cl_LF.h" @@ -186,6 +187,178 @@ const cl_LF eval_rational_series (uintC N, const cl_pq_series& args, uintC len) return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS); } } + +static void eval_pq_series_aux (uintC N1, uintC N2, + cl_pq_series_stream& args, + cl_I* P, cl_I* Q, cl_I* T) +{ + switch (N2 - N1) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pq_series_term v0 = args.next(); // [N1] + if (P) { *P = v0.p; } + *Q = v0.q; + *T = v0.p; + break; + } + case 2: { + var cl_pq_series_term v0 = args.next(); // [N1] + var cl_pq_series_term v1 = args.next(); // [N1+1] + var cl_I p01 = v0.p * v1.p; + if (P) { *P = p01; } + *Q = v0.q * v1.q; + *T = v1.q * v0.p + + p01; + break; + } + case 3: { + var cl_pq_series_term v0 = args.next(); // [N1] + var cl_pq_series_term v1 = args.next(); // [N1+1] + var cl_pq_series_term v2 = args.next(); // [N1+2] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + if (P) { *P = p012; } + var cl_I q12 = v1.q * v2.q; + *Q = v0.q * q12; + *T = q12 * v0.p + + v2.q * p01 + + p012; + break; + } + case 4: { + var cl_pq_series_term v0 = args.next(); // [N1] + var cl_pq_series_term v1 = args.next(); // [N1+1] + var cl_pq_series_term v2 = args.next(); // [N1+2] + var cl_pq_series_term v3 = args.next(); // [N1+3] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + var cl_I p0123 = p012 * v3.p; + if (P) { *P = p0123; } + var cl_I q23 = v2.q * v3.q; + var cl_I q123 = v1.q * q23; + *Q = v0.q * q123; + *T = q123 * v0.p + + q23 * p01 + + v3.q * p012 + + p0123; + break; + } + default: { + var uintC Nm = (N1+N2)/2; // midpoint + // Compute left part. + var cl_I LP, LQ, LT; + eval_pq_series_aux(N1,Nm,args,&LP,&LQ,<); + // Compute right part. + var cl_I RP, RQ, RT; + eval_pq_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT); + // Put together partial results. + if (P) { *P = LP*RP; } + *Q = LQ*RQ; + // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT. + *T = RQ*LT + LP*RT; + break; + } + } +} + +const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_I Q, T; + eval_pq_series_aux(0,N,args,NULL,&Q,&T); + return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len); +} + +static void eval_pq_series_aux (uintC N1, uintC N2, + cl_pq_series_stream& args, + cl_R* P, cl_R* Q, cl_R* T, + uintC trunclen) +{ + switch (N2 - N1) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pq_series_term v0 = args.next(); // [N1] + if (P) { *P = v0.p; } + *Q = v0.q; + *T = v0.p; + break; + } + case 2: { + var cl_pq_series_term v0 = args.next(); // [N1] + var cl_pq_series_term v1 = args.next(); // [N1+1] + var cl_I p01 = v0.p * v1.p; + if (P) { *P = p01; } + *Q = v0.q * v1.q; + *T = v1.q * v0.p + + p01; + break; + } + case 3: { + var cl_pq_series_term v0 = args.next(); // [N1] + var cl_pq_series_term v1 = args.next(); // [N1+1] + var cl_pq_series_term v2 = args.next(); // [N1+2] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + if (P) { *P = p012; } + var cl_I q12 = v1.q * v2.q; + *Q = v0.q * q12; + *T = q12 * v0.p + + v2.q * p01 + + p012; + break; + } + case 4: { + var cl_pq_series_term v0 = args.next(); // [N1] + var cl_pq_series_term v1 = args.next(); // [N1+1] + var cl_pq_series_term v2 = args.next(); // [N1+2] + var cl_pq_series_term v3 = args.next(); // [N1+3] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + var cl_I p0123 = p012 * v3.p; + if (P) { *P = p0123; } + var cl_I q23 = v2.q * v3.q; + var cl_I q123 = v1.q * q23; + *Q = v0.q * q123; + *T = q123 * v0.p + + q23 * p01 + + v3.q * p012 + + p0123; + break; + } + default: { + var uintC Nm = (N1+N2)/2; // midpoint + // Compute left part. + var cl_R LP, LQ, LT; + eval_pq_series_aux(N1,Nm,args,&LP,&LQ,<,trunclen); + // Compute right part. + var cl_R RP, RQ, RT; + eval_pq_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT,trunclen); + // Put together partial results. + if (P) { + *P = LP*RP; + truncate_precision(*P,trunclen); + } + *Q = LQ*RQ; + truncate_precision(*Q,trunclen); + // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT. + *T = RQ*LT + LP*RT; + truncate_precision(*T,trunclen); + break; + } + } +} + +const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len, uintC trunclen) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_R Q, T; + eval_pq_series_aux(0,N,args,NULL,&Q,&T,trunclen); + return cl_R_to_LF(T,len) / cl_R_to_LF(Q,len); +} // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))): // O(log(N)^2*M(N)). diff --git a/src/float/transcendental/cl_LF_ratseries_pqa.cc b/src/float/transcendental/cl_LF_ratseries_pqa.cc index ee0dba1..f11e2a7 100644 --- a/src/float/transcendental/cl_LF_ratseries_pqa.cc +++ b/src/float/transcendental/cl_LF_ratseries_pqa.cc @@ -11,6 +11,7 @@ #include "cln/lfloat.h" #include "cln/integer.h" +#include "cln/real.h" #include "cln/exception.h" #include "cl_LF.h" @@ -186,6 +187,178 @@ const cl_LF eval_rational_series (uintC N, const cl_pqa_series& args, uintC len) return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS); } } + +static void eval_pqa_series_aux (uintC N1, uintC N2, + cl_pqa_series_stream& args, + cl_I* P, cl_I* Q, cl_I* T) +{ + switch (N2 - N1) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pqa_series_term v0 = args.next(); // [N1] + if (P) { *P = v0.p; } + *Q = v0.q; + *T = v0.a * v0.p; + break; + } + case 2: { + var cl_pqa_series_term v0 = args.next(); // [N1] + var cl_pqa_series_term v1 = args.next(); // [N1+1] + var cl_I p01 = v0.p * v1.p; + if (P) { *P = p01; } + *Q = v0.q * v1.q; + *T = v1.q * v0.a * v0.p + + v1.a * p01; + break; + } + case 3: { + var cl_pqa_series_term v0 = args.next(); // [N1] + var cl_pqa_series_term v1 = args.next(); // [N1+1] + var cl_pqa_series_term v2 = args.next(); // [N1+2] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + if (P) { *P = p012; } + var cl_I q12 = v1.q * v2.q; + *Q = v0.q * q12; + *T = q12 * v0.a * v0.p + + v2.q * v1.a * p01 + + v2.a * p012; + break; + } + case 4: { + var cl_pqa_series_term v0 = args.next(); // [N1] + var cl_pqa_series_term v1 = args.next(); // [N1+1] + var cl_pqa_series_term v2 = args.next(); // [N1+2] + var cl_pqa_series_term v3 = args.next(); // [N1+3] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + var cl_I p0123 = p012 * v3.p; + if (P) { *P = p0123; } + var cl_I q23 = v2.q * v3.q; + var cl_I q123 = v1.q * q23; + *Q = v0.q * q123; + *T = q123 * v0.a * v0.p + + q23 * v1.a * p01 + + v3.q * v2.a * p012 + + v3.a * p0123; + break; + } + default: { + var uintC Nm = (N1+N2)/2; // midpoint + // Compute left part. + var cl_I LP, LQ, LT; + eval_pqa_series_aux(N1,Nm,args,&LP,&LQ,<); + // Compute right part. + var cl_I RP, RQ, RT; + eval_pqa_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT); + // Put together partial results. + if (P) { *P = LP*RP; } + *Q = LQ*RQ; + // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT. + *T = RQ*LT + LP*RT; + break; + } + } +} + +const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_I Q, T; + eval_pqa_series_aux(0,N,args,NULL,&Q,&T); + return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len); +} + +static void eval_pqa_series_aux (uintC N1, uintC N2, + cl_pqa_series_stream& args, + cl_R* P, cl_R* Q, cl_R* T, + uintC trunclen) +{ + switch (N2 - N1) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pqa_series_term v0 = args.next(); // [N1] + if (P) { *P = v0.p; } + *Q = v0.q; + *T = v0.a * v0.p; + break; + } + case 2: { + var cl_pqa_series_term v0 = args.next(); // [N1] + var cl_pqa_series_term v1 = args.next(); // [N1+1] + var cl_I p01 = v0.p * v1.p; + if (P) { *P = p01; } + *Q = v0.q * v1.q; + *T = v1.q * v0.a * v0.p + + v1.a * p01; + break; + } + case 3: { + var cl_pqa_series_term v0 = args.next(); // [N1] + var cl_pqa_series_term v1 = args.next(); // [N1+1] + var cl_pqa_series_term v2 = args.next(); // [N1+2] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + if (P) { *P = p012; } + var cl_I q12 = v1.q * v2.q; + *Q = v0.q * q12; + *T = q12 * v0.a * v0.p + + v2.q * v1.a * p01 + + v2.a * p012; + break; + } + case 4: { + var cl_pqa_series_term v0 = args.next(); // [N1] + var cl_pqa_series_term v1 = args.next(); // [N1+1] + var cl_pqa_series_term v2 = args.next(); // [N1+2] + var cl_pqa_series_term v3 = args.next(); // [N1+3] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + var cl_I p0123 = p012 * v3.p; + if (P) { *P = p0123; } + var cl_I q23 = v2.q * v3.q; + var cl_I q123 = v1.q * q23; + *Q = v0.q * q123; + *T = q123 * v0.a * v0.p + + q23 * v1.a * p01 + + v3.q * v2.a * p012 + + v3.a * p0123; + break; + } + default: { + var uintC Nm = (N1+N2)/2; // midpoint + // Compute left part. + var cl_R LP, LQ, LT; + eval_pqa_series_aux(N1,Nm,args,&LP,&LQ,<,trunclen); + // Compute right part. + var cl_R RP, RQ, RT; + eval_pqa_series_aux(Nm,N2,args,(P?&RP:(cl_R*)0),&RQ,&RT,trunclen); + // Put together partial results. + if (P) { + *P = LP*RP; + truncate_precision(*P,trunclen); + } + *Q = LQ*RQ; + truncate_precision(*Q,trunclen); + // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT. + *T = RQ*LT + LP*RT; + truncate_precision(*T,trunclen); + break; + } + } +} + +const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len, uintC trunclen) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_R Q, T; + eval_pqa_series_aux(0,N,args,NULL,&Q,&T,trunclen); + return cl_R_to_LF(T,len) / cl_R_to_LF(Q,len); +} // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))): // O(log(N)^2*M(N)). diff --git a/src/float/transcendental/cl_LF_ratseries_pqab.cc b/src/float/transcendental/cl_LF_ratseries_pqab.cc index 2710042..8c624fd 100644 --- a/src/float/transcendental/cl_LF_ratseries_pqab.cc +++ b/src/float/transcendental/cl_LF_ratseries_pqab.cc @@ -11,6 +11,7 @@ #include "cln/lfloat.h" #include "cln/integer.h" +#include "cln/real.h" #include "cln/exception.h" #include "cl_LF.h" @@ -202,6 +203,195 @@ const cl_LF eval_rational_series (uintC N, const cl_pqab_series& args, uintC len return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(B*Q,len),QS); } } + +static void eval_pqab_series_aux (uintC N1, uintC N2, + cl_pqab_series_stream& args, + cl_I* P, cl_I* Q, cl_I* B, cl_I* T) +{ + switch (N2 - N1) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pqab_series_term v0 = args.next(); // [N1] + if (P) { *P = v0.p; } + *Q = v0.q; + *B = v0.b; + *T = v0.a * v0.p; + break; + } + case 2: { + var cl_pqab_series_term v0 = args.next(); // [N1] + var cl_pqab_series_term v1 = args.next(); // [N1+1] + var cl_I p01 = v0.p * v1.p; + if (P) { *P = p01; } + *Q = v0.q * v1.q; + *B = v0.b * v1.b; + *T = v1.b * v1.q * v0.a * v0.p + + v0.b * v1.a * p01; + break; + } + case 3: { + var cl_pqab_series_term v0 = args.next(); // [N1] + var cl_pqab_series_term v1 = args.next(); // [N1+1] + var cl_pqab_series_term v2 = args.next(); // [N1+2] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + if (P) { *P = p012; } + var cl_I q12 = v1.q * v2.q; + *Q = v0.q * q12; + var cl_I b12 = v1.b * v2.b; + *B = v0.b * b12; + *T = b12 * q12 * v0.a * v0.p + + v0.b * (v2.b * v2.q * v1.a * p01 + + v1.b * v2.a * p012); + break; + } + case 4: { + var cl_pqab_series_term v0 = args.next(); // [N1] + var cl_pqab_series_term v1 = args.next(); // [N1+1] + var cl_pqab_series_term v2 = args.next(); // [N1+2] + var cl_pqab_series_term v3 = args.next(); // [N1+3] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + var cl_I p0123 = p012 * v3.p; + if (P) { *P = p0123; } + var cl_I q23 = v2.q * v3.q; + var cl_I q123 = v1.q * q23; + *Q = v0.q * q123; + var cl_I b01 = v0.b * v1.b; + var cl_I b23 = v2.b * v3.b; + *B = b01 * b23; + *T = b23 * (v1.b * q123 * v0.a * v0.p + + v0.b * q23 * v1.a * p01) + + b01 * (v3.b * v3.q * v2.a * p012 + + v2.b * v3.a * p0123); + break; + } + default: { + var uintC Nm = (N1+N2)/2; // midpoint + // Compute left part. + var cl_I LP, LQ, LB, LT; + eval_pqab_series_aux(N1,Nm,args,&LP,&LQ,&LB,<); + // Compute right part. + var cl_I RP, RQ, RB, RT; + eval_pqab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT); + // Put together partial results. + if (P) { *P = LP*RP; } + *Q = LQ*RQ; + *B = LB*RB; + // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT. + *T = RB*RQ*LT + LB*LP*RT; + break; + } + } +} + +const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_I Q, B, T; + eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T); + return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len); +} + +static void eval_pqab_series_aux (uintC N1, uintC N2, + cl_pqab_series_stream& args, + cl_R* P, cl_R* Q, cl_R* B, cl_R* T, + uintC trunclen) +{ + switch (N2 - N1) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pqab_series_term v0 = args.next(); // [N1] + if (P) { *P = v0.p; } + *Q = v0.q; + *B = v0.b; + *T = v0.a * v0.p; + break; + } + case 2: { + var cl_pqab_series_term v0 = args.next(); // [N1] + var cl_pqab_series_term v1 = args.next(); // [N1+1] + var cl_I p01 = v0.p * v1.p; + if (P) { *P = p01; } + *Q = v0.q * v1.q; + *B = v0.b * v1.b; + *T = v1.b * v1.q * v0.a * v0.p + + v0.b * v1.a * p01; + break; + } + case 3: { + var cl_pqab_series_term v0 = args.next(); // [N1] + var cl_pqab_series_term v1 = args.next(); // [N1+1] + var cl_pqab_series_term v2 = args.next(); // [N1+2] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + if (P) { *P = p012; } + var cl_I q12 = v1.q * v2.q; + *Q = v0.q * q12; + var cl_I b12 = v1.b * v2.b; + *B = v0.b * b12; + *T = b12 * q12 * v0.a * v0.p + + v0.b * (v2.b * v2.q * v1.a * p01 + + v1.b * v2.a * p012); + break; + } + case 4: { + var cl_pqab_series_term v0 = args.next(); // [N1] + var cl_pqab_series_term v1 = args.next(); // [N1+1] + var cl_pqab_series_term v2 = args.next(); // [N1+2] + var cl_pqab_series_term v3 = args.next(); // [N1+3] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + var cl_I p0123 = p012 * v3.p; + if (P) { *P = p0123; } + var cl_I q23 = v2.q * v3.q; + var cl_I q123 = v1.q * q23; + *Q = v0.q * q123; + var cl_I b01 = v0.b * v1.b; + var cl_I b23 = v2.b * v3.b; + *B = b01 * b23; + *T = b23 * (v1.b * q123 * v0.a * v0.p + + v0.b * q23 * v1.a * p01) + + b01 * (v3.b * v3.q * v2.a * p012 + + v2.b * v3.a * p0123); + break; + } + default: { + var uintC Nm = (N1+N2)/2; // midpoint + // Compute left part. + var cl_R LP, LQ, LB, LT; + eval_pqab_series_aux(N1,Nm,args,&LP,&LQ,&LB,<,trunclen); + // Compute right part. + var cl_R RP, RQ, RB, RT; + eval_pqab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT,trunclen); + // Put together partial results. + if (P) { + *P = LP*RP; + truncate_precision(*P,trunclen); + } + *Q = LQ*RQ; + truncate_precision(*Q,trunclen); + *B = LB*RB; + truncate_precision(*B,trunclen); + // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT. + *T = RB*RQ*LT + LB*LP*RT; + truncate_precision(*T,trunclen); + break; + } + } +} + +const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len, uintC trunclen) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_R Q, B, T; + eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T,trunclen); + return cl_R_to_LF(T,len) / cl_R_to_LF(B*Q,len); +} // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))): // O(log(N)^2*M(N)). diff --git a/src/float/transcendental/cl_LF_ratseries_pqb.cc b/src/float/transcendental/cl_LF_ratseries_pqb.cc index 2ef3f0c..84dbf79 100644 --- a/src/float/transcendental/cl_LF_ratseries_pqb.cc +++ b/src/float/transcendental/cl_LF_ratseries_pqb.cc @@ -11,6 +11,7 @@ #include "cln/lfloat.h" #include "cln/integer.h" +#include "cln/real.h" #include "cln/exception.h" #include "cl_LF.h" @@ -202,6 +203,196 @@ const cl_LF eval_rational_series (uintC N, const cl_pqb_series& args, uintC len) return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(B*Q,len),QS); } } + +static void eval_pqb_series_aux (uintC N1, uintC N2, + cl_pqb_series_stream& args, + cl_I* P, cl_I* Q, cl_I* B, cl_I* T) +{ + switch (N2 - N1) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pqb_series_term v0 = args.next(); // [N1] + if (P) { *P = v0.p; } + *Q = v0.q; + *B = v0.b; + *T = v0.p; + break; + } + case 2: { + var cl_pqb_series_term v0 = args.next(); // [N1] + var cl_pqb_series_term v1 = args.next(); // [N1+1] + var cl_I p01 = v0.p * v1.p; + if (P) { *P = p01; } + *Q = v0.q * v1.q; + *B = v0.b * v1.b; + *T = v1.b * v1.q * v0.p + + v0.b * p01; + break; + } + case 3: { + var cl_pqb_series_term v0 = args.next(); // [N1] + var cl_pqb_series_term v1 = args.next(); // [N1+1] + var cl_pqb_series_term v2 = args.next(); // [N1+2] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + if (P) { *P = p012; } + var cl_I q12 = v1.q * v2.q; + *Q = v0.q * q12; + var cl_I b12 = v1.b * v2.b; + *B = v0.b * b12; + *T = b12 * q12 * v0.p + + v0.b * (v2.b * v2.q * p01 + + v1.b * p012); + break; + } + case 4: { + var cl_pqb_series_term v0 = args.next(); // [N1] + var cl_pqb_series_term v1 = args.next(); // [N1+1] + var cl_pqb_series_term v2 = args.next(); // [N1+2] + var cl_pqb_series_term v3 = args.next(); // [N1+3] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + var cl_I p0123 = p012 * v3.p; + if (P) { *P = p0123; } + var cl_I q23 = v2.q * v3.q; + var cl_I q123 = v1.q * q23; + *Q = v0.q * q123; + var cl_I b01 = v0.b * v1.b; + var cl_I b23 = v2.b * v3.b; + *B = b01 * b23; + *T = b23 * (v1.b * q123 * v0.p + + v0.b * q23 * p01) + + b01 * (v3.b * v3.q * p012 + + v2.b * p0123); + break; + } + default: { + var uintC Nm = (N1+N2)/2; // midpoint + // Compute left part. + var cl_I LP, LQ, LB, LT; + eval_pqb_series_aux(N1,Nm,args,&LP,&LQ,&LB,<); + // Compute right part. + var cl_I RP, RQ, RB, RT; + eval_pqb_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT); + // Put together partial results. + if (P) { *P = LP*RP; } + *Q = LQ*RQ; + *B = LB*RB; + // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT. + *T = RB*RQ*LT + LB*LP*RT; + break; + } + } +} + +const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_I Q, B, T; + eval_pqb_series_aux(0,N,args,NULL,&Q,&B,&T); + return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len); +} + +static void eval_pqb_series_aux (uintC N1, uintC N2, + cl_pqb_series_stream& args, + cl_R* P, cl_R* Q, cl_R* B, cl_R* T, + uintC trunclen) +{ + switch (N2 - N1) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pqb_series_term v0 = args.next(); // [N1] + if (P) { *P = v0.p; } + *Q = v0.q; + *B = v0.b; + *T = v0.p; + break; + } + case 2: { + var cl_pqb_series_term v0 = args.next(); // [N1] + var cl_pqb_series_term v1 = args.next(); // [N1+1] + var cl_I p01 = v0.p * v1.p; + if (P) { *P = p01; } + *Q = v0.q * v1.q; + *B = v0.b * v1.b; + *T = v1.b * v1.q * v0.p + + v0.b * p01; + break; + } + case 3: { + var cl_pqb_series_term v0 = args.next(); // [N1] + var cl_pqb_series_term v1 = args.next(); // [N1+1] + var cl_pqb_series_term v2 = args.next(); // [N1+2] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + if (P) { *P = p012; } + var cl_I q12 = v1.q * v2.q; + *Q = v0.q * q12; + var cl_I b12 = v1.b * v2.b; + *B = v0.b * b12; + *T = b12 * q12 * v0.p + + v0.b * (v2.b * v2.q * p01 + + v1.b * p012); + break; + } + case 4: { + var cl_pqb_series_term v0 = args.next(); // [N1] + var cl_pqb_series_term v1 = args.next(); // [N1+1] + var cl_pqb_series_term v2 = args.next(); // [N1+2] + var cl_pqb_series_term v3 = args.next(); // [N1+3] + var cl_I p01 = v0.p * v1.p; + var cl_I p012 = p01 * v2.p; + var cl_I p0123 = p012 * v3.p; + if (P) { *P = p0123; } + var cl_I q23 = v2.q * v3.q; + var cl_I q123 = v1.q * q23; + *Q = v0.q * q123; + var cl_I b01 = v0.b * v1.b; + var cl_I b23 = v2.b * v3.b; + *B = b01 * b23; + *T = b23 * (v1.b * q123 * v0.p + + v0.b * q23 * p01) + + b01 * (v3.b * v3.q * p012 + + v2.b * p0123); + break; + } + default: { + var uintC Nm = (N1+N2)/2; // midpoint + // Compute left part. + var cl_R LP, LQ, LB, LT; + eval_pqb_series_aux(N1,Nm,args,&LP,&LQ,&LB,<,trunclen); + // Compute right part. + var cl_R RP, RQ, RB, RT; + eval_pqb_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT,trunclen); + // Put together partial results. + if (P) { + *P = LP*RP; + truncate_precision(*P,trunclen); + } + *Q = LQ*RQ; + truncate_precision(*Q,trunclen); + *B = LB*RB; + truncate_precision(*B,trunclen); + // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT. + *T = RB*RQ*LT + LB*LP*RT; + truncate_precision(*T,trunclen); + break; + } + } +} + +const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len, uintC trunclen) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_R Q, B, T; + eval_pqb_series_aux(0,N,args,NULL,&Q,&B,&T,trunclen); + return cl_R_to_LF(T,len) / cl_R_to_LF(B*Q,len); +} + // Bit complexity (if p(n), q(n), a(n), b(n) have length O(log(n))): // O(log(N)^2*M(N)). diff --git a/src/float/transcendental/cl_LF_ratseries_stream_pq.cc b/src/float/transcendental/cl_LF_ratseries_stream_pq.cc deleted file mode 100644 index 91b3391..0000000 --- a/src/float/transcendental/cl_LF_ratseries_stream_pq.cc +++ /dev/null @@ -1,102 +0,0 @@ -// eval_rational_series(). - -// General includes. -#include "cl_sysdep.h" - -// Specification. -#include "cl_LF_tran.h" - - -// Implementation. - -#include "cln/lfloat.h" -#include "cln/integer.h" -#include "cln/exception.h" -#include "cl_LF.h" - -namespace cln { - -static void eval_pq_series_aux (uintC N1, uintC N2, - cl_pq_series_stream& args, - cl_I* P, cl_I* Q, cl_I* T) -{ - switch (N2 - N1) { - case 0: - throw runtime_exception(); break; - case 1: { - var cl_pq_series_term v0 = args.next(); // [N1] - if (P) { *P = v0.p; } - *Q = v0.q; - *T = v0.p; - break; - } - case 2: { - var cl_pq_series_term v0 = args.next(); // [N1] - var cl_pq_series_term v1 = args.next(); // [N1+1] - var cl_I p01 = v0.p * v1.p; - if (P) { *P = p01; } - *Q = v0.q * v1.q; - *T = v1.q * v0.p - + p01; - break; - } - case 3: { - var cl_pq_series_term v0 = args.next(); // [N1] - var cl_pq_series_term v1 = args.next(); // [N1+1] - var cl_pq_series_term v2 = args.next(); // [N1+2] - var cl_I p01 = v0.p * v1.p; - var cl_I p012 = p01 * v2.p; - if (P) { *P = p012; } - var cl_I q12 = v1.q * v2.q; - *Q = v0.q * q12; - *T = q12 * v0.p - + v2.q * p01 - + p012; - break; - } - case 4: { - var cl_pq_series_term v0 = args.next(); // [N1] - var cl_pq_series_term v1 = args.next(); // [N1+1] - var cl_pq_series_term v2 = args.next(); // [N1+2] - var cl_pq_series_term v3 = args.next(); // [N1+3] - var cl_I p01 = v0.p * v1.p; - var cl_I p012 = p01 * v2.p; - var cl_I p0123 = p012 * v3.p; - if (P) { *P = p0123; } - var cl_I q23 = v2.q * v3.q; - var cl_I q123 = v1.q * q23; - *Q = v0.q * q123; - *T = q123 * v0.p - + q23 * p01 - + v3.q * p012 - + p0123; - break; - } - default: { - var uintC Nm = (N1+N2)/2; // midpoint - // Compute left part. - var cl_I LP, LQ, LT; - eval_pq_series_aux(N1,Nm,args,&LP,&LQ,<); - // Compute right part. - var cl_I RP, RQ, RT; - eval_pq_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT); - // Put together partial results. - if (P) { *P = LP*RP; } - *Q = LQ*RQ; - // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT. - *T = RQ*LT + LP*RT; - break; - } - } -} - -const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len) -{ - if (N==0) - return cl_I_to_LF(0,len); - var cl_I Q, T; - eval_pq_series_aux(0,N,args,NULL,&Q,&T); - return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len); -} - -} // namespace cln diff --git a/src/float/transcendental/cl_LF_ratseries_stream_pqa.cc b/src/float/transcendental/cl_LF_ratseries_stream_pqa.cc deleted file mode 100644 index 3a2f3f0..0000000 --- a/src/float/transcendental/cl_LF_ratseries_stream_pqa.cc +++ /dev/null @@ -1,102 +0,0 @@ -// eval_rational_series(). - -// General includes. -#include "cl_sysdep.h" - -// Specification. -#include "cl_LF_tran.h" - - -// Implementation. - -#include "cln/lfloat.h" -#include "cln/integer.h" -#include "cln/exception.h" -#include "cl_LF.h" - -namespace cln { - -static void eval_pqa_series_aux (uintC N1, uintC N2, - cl_pqa_series_stream& args, - cl_I* P, cl_I* Q, cl_I* T) -{ - switch (N2 - N1) { - case 0: - throw runtime_exception(); break; - case 1: { - var cl_pqa_series_term v0 = args.next(); // [N1] - if (P) { *P = v0.p; } - *Q = v0.q; - *T = v0.a * v0.p; - break; - } - case 2: { - var cl_pqa_series_term v0 = args.next(); // [N1] - var cl_pqa_series_term v1 = args.next(); // [N1+1] - var cl_I p01 = v0.p * v1.p; - if (P) { *P = p01; } - *Q = v0.q * v1.q; - *T = v1.q * v0.a * v0.p - + v1.a * p01; - break; - } - case 3: { - var cl_pqa_series_term v0 = args.next(); // [N1] - var cl_pqa_series_term v1 = args.next(); // [N1+1] - var cl_pqa_series_term v2 = args.next(); // [N1+2] - var cl_I p01 = v0.p * v1.p; - var cl_I p012 = p01 * v2.p; - if (P) { *P = p012; } - var cl_I q12 = v1.q * v2.q; - *Q = v0.q * q12; - *T = q12 * v0.a * v0.p - + v2.q * v1.a * p01 - + v2.a * p012; - break; - } - case 4: { - var cl_pqa_series_term v0 = args.next(); // [N1] - var cl_pqa_series_term v1 = args.next(); // [N1+1] - var cl_pqa_series_term v2 = args.next(); // [N1+2] - var cl_pqa_series_term v3 = args.next(); // [N1+3] - var cl_I p01 = v0.p * v1.p; - var cl_I p012 = p01 * v2.p; - var cl_I p0123 = p012 * v3.p; - if (P) { *P = p0123; } - var cl_I q23 = v2.q * v3.q; - var cl_I q123 = v1.q * q23; - *Q = v0.q * q123; - *T = q123 * v0.a * v0.p - + q23 * v1.a * p01 - + v3.q * v2.a * p012 - + v3.a * p0123; - break; - } - default: { - var uintC Nm = (N1+N2)/2; // midpoint - // Compute left part. - var cl_I LP, LQ, LT; - eval_pqa_series_aux(N1,Nm,args,&LP,&LQ,<); - // Compute right part. - var cl_I RP, RQ, RT; - eval_pqa_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RT); - // Put together partial results. - if (P) { *P = LP*RP; } - *Q = LQ*RQ; - // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT. - *T = RQ*LT + LP*RT; - break; - } - } -} - -const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len) -{ - if (N==0) - return cl_I_to_LF(0,len); - var cl_I Q, T; - eval_pqa_series_aux(0,N,args,NULL,&Q,&T); - return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len); -} - -} // namespace cln diff --git a/src/float/transcendental/cl_LF_ratseries_stream_pqab.cc b/src/float/transcendental/cl_LF_ratseries_stream_pqab.cc deleted file mode 100644 index 8203c9f..0000000 --- a/src/float/transcendental/cl_LF_ratseries_stream_pqab.cc +++ /dev/null @@ -1,110 +0,0 @@ -// eval_rational_series(). - -// General includes. -#include "cl_sysdep.h" - -// Specification. -#include "cl_LF_tran.h" - - -// Implementation. - -#include "cln/lfloat.h" -#include "cln/integer.h" -#include "cln/exception.h" -#include "cl_LF.h" - -namespace cln { - -static void eval_pqab_series_aux (uintC N1, uintC N2, - cl_pqab_series_stream& args, - cl_I* P, cl_I* Q, cl_I* B, cl_I* T) -{ - switch (N2 - N1) { - case 0: - throw runtime_exception(); break; - case 1: { - var cl_pqab_series_term v0 = args.next(); // [N1] - if (P) { *P = v0.p; } - *Q = v0.q; - *B = v0.b; - *T = v0.a * v0.p; - break; - } - case 2: { - var cl_pqab_series_term v0 = args.next(); // [N1] - var cl_pqab_series_term v1 = args.next(); // [N1+1] - var cl_I p01 = v0.p * v1.p; - if (P) { *P = p01; } - *Q = v0.q * v1.q; - *B = v0.b * v1.b; - *T = v1.b * v1.q * v0.a * v0.p - + v0.b * v1.a * p01; - break; - } - case 3: { - var cl_pqab_series_term v0 = args.next(); // [N1] - var cl_pqab_series_term v1 = args.next(); // [N1+1] - var cl_pqab_series_term v2 = args.next(); // [N1+2] - var cl_I p01 = v0.p * v1.p; - var cl_I p012 = p01 * v2.p; - if (P) { *P = p012; } - var cl_I q12 = v1.q * v2.q; - *Q = v0.q * q12; - var cl_I b12 = v1.b * v2.b; - *B = v0.b * b12; - *T = b12 * q12 * v0.a * v0.p - + v0.b * (v2.b * v2.q * v1.a * p01 - + v1.b * v2.a * p012); - break; - } - case 4: { - var cl_pqab_series_term v0 = args.next(); // [N1] - var cl_pqab_series_term v1 = args.next(); // [N1+1] - var cl_pqab_series_term v2 = args.next(); // [N1+2] - var cl_pqab_series_term v3 = args.next(); // [N1+3] - var cl_I p01 = v0.p * v1.p; - var cl_I p012 = p01 * v2.p; - var cl_I p0123 = p012 * v3.p; - if (P) { *P = p0123; } - var cl_I q23 = v2.q * v3.q; - var cl_I q123 = v1.q * q23; - *Q = v0.q * q123; - var cl_I b01 = v0.b * v1.b; - var cl_I b23 = v2.b * v3.b; - *B = b01 * b23; - *T = b23 * (v1.b * q123 * v0.a * v0.p - + v0.b * q23 * v1.a * p01) - + b01 * (v3.b * v3.q * v2.a * p012 - + v2.b * v3.a * p0123); - break; - } - default: { - var uintC Nm = (N1+N2)/2; // midpoint - // Compute left part. - var cl_I LP, LQ, LB, LT; - eval_pqab_series_aux(N1,Nm,args,&LP,&LQ,&LB,<); - // Compute right part. - var cl_I RP, RQ, RB, RT; - eval_pqab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT); - // Put together partial results. - if (P) { *P = LP*RP; } - *Q = LQ*RQ; - *B = LB*RB; - // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT. - *T = RB*RQ*LT + LB*LP*RT; - break; - } - } -} - -const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len) -{ - if (N==0) - return cl_I_to_LF(0,len); - var cl_I Q, B, T; - eval_pqab_series_aux(0,N,args,NULL,&Q,&B,&T); - return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len); -} - -} // namespace cln diff --git a/src/float/transcendental/cl_LF_ratseries_stream_pqb.cc b/src/float/transcendental/cl_LF_ratseries_stream_pqb.cc deleted file mode 100644 index 12925b7..0000000 --- a/src/float/transcendental/cl_LF_ratseries_stream_pqb.cc +++ /dev/null @@ -1,110 +0,0 @@ -// eval_rational_series(). - -// General includes. -#include "cl_sysdep.h" - -// Specification. -#include "cl_LF_tran.h" - - -// Implementation. - -#include "cln/lfloat.h" -#include "cln/integer.h" -#include "cln/exception.h" -#include "cl_LF.h" - -namespace cln { - -static void eval_pqb_series_aux (uintC N1, uintC N2, - cl_pqb_series_stream& args, - cl_I* P, cl_I* Q, cl_I* B, cl_I* T) -{ - switch (N2 - N1) { - case 0: - throw runtime_exception(); break; - case 1: { - var cl_pqb_series_term v0 = args.next(); // [N1] - if (P) { *P = v0.p; } - *Q = v0.q; - *B = v0.b; - *T = v0.p; - break; - } - case 2: { - var cl_pqb_series_term v0 = args.next(); // [N1] - var cl_pqb_series_term v1 = args.next(); // [N1+1] - var cl_I p01 = v0.p * v1.p; - if (P) { *P = p01; } - *Q = v0.q * v1.q; - *B = v0.b * v1.b; - *T = v1.b * v1.q * v0.p - + v0.b * p01; - break; - } - case 3: { - var cl_pqb_series_term v0 = args.next(); // [N1] - var cl_pqb_series_term v1 = args.next(); // [N1+1] - var cl_pqb_series_term v2 = args.next(); // [N1+2] - var cl_I p01 = v0.p * v1.p; - var cl_I p012 = p01 * v2.p; - if (P) { *P = p012; } - var cl_I q12 = v1.q * v2.q; - *Q = v0.q * q12; - var cl_I b12 = v1.b * v2.b; - *B = v0.b * b12; - *T = b12 * q12 * v0.p - + v0.b * (v2.b * v2.q * p01 - + v1.b * p012); - break; - } - case 4: { - var cl_pqb_series_term v0 = args.next(); // [N1] - var cl_pqb_series_term v1 = args.next(); // [N1+1] - var cl_pqb_series_term v2 = args.next(); // [N1+2] - var cl_pqb_series_term v3 = args.next(); // [N1+3] - var cl_I p01 = v0.p * v1.p; - var cl_I p012 = p01 * v2.p; - var cl_I p0123 = p012 * v3.p; - if (P) { *P = p0123; } - var cl_I q23 = v2.q * v3.q; - var cl_I q123 = v1.q * q23; - *Q = v0.q * q123; - var cl_I b01 = v0.b * v1.b; - var cl_I b23 = v2.b * v3.b; - *B = b01 * b23; - *T = b23 * (v1.b * q123 * v0.p - + v0.b * q23 * p01) - + b01 * (v3.b * v3.q * p012 - + v2.b * p0123); - break; - } - default: { - var uintC Nm = (N1+N2)/2; // midpoint - // Compute left part. - var cl_I LP, LQ, LB, LT; - eval_pqb_series_aux(N1,Nm,args,&LP,&LQ,&LB,<); - // Compute right part. - var cl_I RP, RQ, RB, RT; - eval_pqb_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RB,&RT); - // Put together partial results. - if (P) { *P = LP*RP; } - *Q = LQ*RQ; - *B = LB*RB; - // S = LS + LP/LQ * RS, so T = RB*RQ*LT + LB*LP*RT. - *T = RB*RQ*LT + LB*LP*RT; - break; - } - } -} - -const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len) -{ - if (N==0) - return cl_I_to_LF(0,len); - var cl_I Q, B, T; - eval_pqb_series_aux(0,N,args,NULL,&Q,&B,&T); - return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len); -} - -} // namespace cln diff --git a/src/float/transcendental/cl_LF_ratsumseries_pqcd.cc b/src/float/transcendental/cl_LF_ratsumseries_pqcd.cc index 009beea..68b93d4 100644 --- a/src/float/transcendental/cl_LF_ratsumseries_pqcd.cc +++ b/src/float/transcendental/cl_LF_ratsumseries_pqcd.cc @@ -11,6 +11,7 @@ #include "cln/lfloat.h" #include "cln/integer.h" +#include "cln/real.h" #include "cl_LF.h" namespace cln { @@ -19,7 +20,7 @@ const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_term* args, uintC len) { if (N==0) return cl_I_to_LF(0,len); - var cl_pqcd_series_result sums; + var cl_pqcd_series_result sums; eval_pqcd_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 @@ -28,4 +29,30 @@ const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_term* args, uintC len) cl_I_to_LF(sums.V,len) / The(cl_LF)(sums.D * cl_I_to_LF(sums.T,len)); } +const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_pqcd_series_result sums; + eval_pqcd_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)); +} + +const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len, uintC trunclen) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_pqcd_series_result sums; + eval_pqcd_series_aux(N,args,sums,trunclen); + // 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_R_to_LF(sums.V,len) / The(cl_LF)(sums.D * cl_R_to_LF(sums.T,len)); +} + } // namespace cln diff --git a/src/float/transcendental/cl_LF_ratsumseries_pqcd_aux.cc b/src/float/transcendental/cl_LF_ratsumseries_pqcd_aux.cc index 97cc535..91c0900 100644 --- a/src/float/transcendental/cl_LF_ratsumseries_pqcd_aux.cc +++ b/src/float/transcendental/cl_LF_ratsumseries_pqcd_aux.cc @@ -10,11 +10,12 @@ // Implementation. #include "cln/integer.h" +#include "cln/real.h" #include "cln/exception.h" namespace cln { -void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_result& Z, bool rightmost) +void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_result& Z, bool rightmost) { // N = N2-N1 switch (N) { @@ -59,10 +60,10 @@ void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_re default: { var uintC Nm = N/2; // midpoint // Compute left part. - var cl_pqcd_series_result L; + var cl_pqcd_series_result L; eval_pqcd_series_aux(Nm,args+0,L,false); // Compute right part. - var cl_pqcd_series_result R; + var cl_pqcd_series_result R; eval_pqcd_series_aux(N-Nm,args+Nm,R,rightmost); // Put together partial results. if (!rightmost) { Z.P = L.P * R.P; } @@ -80,4 +81,160 @@ void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_re } } +void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result& Z, bool rightmost) +{ + // N = N2-N1 + switch (N) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pqcd_series_term v0 = args.next(); // [N1] + if (!rightmost) { Z.P = v0.p; } + Z.Q = v0.q; + Z.T = v0.p; + if (!rightmost) { Z.C = v0.c; } + Z.D = v0.d; + Z.V = v0.c * v0.p; + break; + } + case 2: { + var cl_pqcd_series_term v0 = args.next(); // [N1] + var cl_pqcd_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; + var cl_I c0d1 = v0.c * v1.d; + var cl_I c1d0 = v1.c * v0.d; + if (!rightmost) { Z.C = c0d1 + c1d0; } + Z.D = v0.d * v1.d; + Z.V = c0d1 * p0q1 + c1d0 * p01; + break; + } + case 3: { + var cl_pqcd_series_term v0 = args.next(); // [N1] + var cl_pqcd_series_term v1 = args.next(); // [N1+1] + var cl_pqcd_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 c0d1 = v0.c * v1.d; + var cl_I c1d0 = v1.c * v0.d; + var cl_I d01 = v0.d * v1.d; + if (!rightmost) { Z.C = (c0d1 + c1d0) * v2.d + v2.c * d01; } + Z.D = d01 * v2.d; + Z.V = v2.d * (v2.q * (c0d1 * p0q1 + c1d0 * p01) + (c0d1 + c1d0) * p012) + v2.c * d01 * p012; + break; + } + default: { + var uintC Nm = N/2; // midpoint + // Compute left part. + var cl_pqcd_series_result L; + eval_pqcd_series_aux(Nm,args,L,false); + // Compute right part. + var cl_pqcd_series_result R; + eval_pqcd_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; + } + } +} + +void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result& Z, uintC trunclen, bool rightmost) +{ + // N = N2-N1 + switch (N) { + case 0: + throw runtime_exception(); break; + case 1: { + var cl_pqcd_series_term v0 = args.next(); // [N1] + if (!rightmost) { Z.P = v0.p; } + Z.Q = v0.q; + Z.T = v0.p; + if (!rightmost) { Z.C = v0.c; } + Z.D = v0.d; + Z.V = v0.c * v0.p; + break; + } + case 2: { + var cl_pqcd_series_term v0 = args.next(); // [N1] + var cl_pqcd_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; + var cl_I c0d1 = v0.c * v1.d; + var cl_I c1d0 = v1.c * v0.d; + if (!rightmost) { Z.C = c0d1 + c1d0; } + Z.D = v0.d * v1.d; + Z.V = c0d1 * p0q1 + c1d0 * p01; + break; + } + case 3: { + var cl_pqcd_series_term v0 = args.next(); // [N1] + var cl_pqcd_series_term v1 = args.next(); // [N1+1] + var cl_pqcd_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 c0d1 = v0.c * v1.d; + var cl_I c1d0 = v1.c * v0.d; + var cl_I d01 = v0.d * v1.d; + if (!rightmost) { Z.C = (c0d1 + c1d0) * v2.d + v2.c * d01; } + Z.D = d01 * v2.d; + Z.V = v2.d * (v2.q * (c0d1 * p0q1 + c1d0 * p01) + (c0d1 + c1d0) * p012) + v2.c * d01 * p012; + break; + } + default: { + var uintC Nm = N/2; // midpoint + // Compute left part. + var cl_pqcd_series_result L; + eval_pqcd_series_aux(Nm,args,L,trunclen,false); + // Compute right part. + var cl_pqcd_series_result R; + eval_pqcd_series_aux(N-Nm,args,R,trunclen,rightmost); + // Put together partial results. + if (!rightmost) { + Z.P = L.P * R.P; + truncate_precision(Z.P,trunclen); + } + Z.Q = L.Q * R.Q; + truncate_precision(Z.Q,trunclen); + // Z.S = L.S + L.P/L.Q*R.S; + var cl_R tmp = L.P * R.T; + Z.T = R.Q * L.T + tmp; + truncate_precision(Z.T,trunclen); + if (!rightmost) { + Z.C = L.C * R.D + L.D * R.C; + truncate_precision(Z.C,trunclen); + } + Z.D = L.D * R.D; + truncate_precision(Z.D,trunclen); + // 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; + truncate_precision(Z.V,trunclen); + break; + } + } +} + } // namespace cln diff --git a/src/float/transcendental/cl_LF_ratsumseries_pqd.cc b/src/float/transcendental/cl_LF_ratsumseries_pqd.cc index 7d70ad7..a47db3a 100644 --- a/src/float/transcendental/cl_LF_ratsumseries_pqd.cc +++ b/src/float/transcendental/cl_LF_ratsumseries_pqd.cc @@ -11,6 +11,7 @@ #include "cln/lfloat.h" #include "cln/integer.h" +#include "cln/real.h" #include "cl_LF.h" namespace cln { @@ -19,7 +20,7 @@ const cl_LF eval_pqd_series (uintC N, cl_pqd_series_term* args, uintC len) { if (N==0) return cl_I_to_LF(0,len); - var cl_pqd_series_result sums; + 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 @@ -28,4 +29,30 @@ const cl_LF eval_pqd_series (uintC N, cl_pqd_series_term* args, uintC len) cl_I_to_LF(sums.V,len) / The(cl_LF)(sums.D * cl_I_to_LF(sums.T,len)); } +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)); +} + +const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len, uintC trunclen) +{ + if (N==0) + return cl_I_to_LF(0,len); + var cl_pqd_series_result sums; + eval_pqd_series_aux(N,args,sums,trunclen); + // 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_R_to_LF(sums.V,len) / The(cl_LF)(sums.D * cl_R_to_LF(sums.T,len)); +} + } // namespace cln diff --git a/src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc b/src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc index 3373fbc..ed3e865 100644 --- a/src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc +++ b/src/float/transcendental/cl_LF_ratsumseries_pqd_aux.cc @@ -10,11 +10,12 @@ // Implementation. #include "cln/integer.h" +#include "cln/real.h" #include "cln/exception.h" namespace cln { -void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_result& Z, bool rightmost) +void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_result& Z, bool rightmost) { // N = N2-N1 switch (N) { @@ -55,10 +56,10 @@ void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_resul default: { var uintC Nm = N/2; // midpoint // Compute left part. - var cl_pqd_series_result L; + var cl_pqd_series_result L; eval_pqd_series_aux(Nm,args+0,L,false); // Compute right part. - var cl_pqd_series_result R; + var cl_pqd_series_result R; eval_pqd_series_aux(N-Nm,args+Nm,R,rightmost); // Put together partial results. if (!rightmost) { Z.P = L.P * R.P; } @@ -76,4 +77,152 @@ void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_resul } } +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; + } + } +} + +void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result& Z, uintC trunclen, 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,trunclen,false); + // Compute right part. + var cl_pqd_series_result R; + eval_pqd_series_aux(N-Nm,args,R,trunclen,rightmost); + // Put together partial results. + if (!rightmost) { + Z.P = L.P * R.P; + truncate_precision(Z.P,trunclen); + } + Z.Q = L.Q * R.Q; + truncate_precision(Z.Q,trunclen); + // Z.S = L.S + L.P/L.Q*R.S; + var cl_R tmp = L.P * R.T; + Z.T = R.Q * L.T + tmp; + truncate_precision(Z.T,trunclen); + if (!rightmost) { + Z.C = L.C * R.D + L.D * R.C; + truncate_precision(Z.C,trunclen); + } + Z.D = L.D * R.D; + truncate_precision(Z.D,trunclen); + // 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; + truncate_precision(Z.V,trunclen); + break; + } + } +} + } // namespace cln diff --git a/src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc b/src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc deleted file mode 100644 index 2f1aaaf..0000000 --- a/src/float/transcendental/cl_LF_ratsumseries_stream_pqd.cc +++ /dev/null @@ -1,31 +0,0 @@ -// 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 deleted file mode 100644 index cc64077..0000000 --- a/src/float/transcendental/cl_LF_ratsumseries_stream_pqd_aux.cc +++ /dev/null @@ -1,86 +0,0 @@ -// 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 b975fa3..3484f9e 100644 --- a/src/float/transcendental/cl_LF_tran.h +++ b/src/float/transcendental/cl_LF_tran.h @@ -4,14 +4,24 @@ #define _CL_LF_TRAN_H #include "cln/integer.h" +#include "cln/integer_ring.h" #include "cln/lfloat.h" +#include "cl_LF.h" namespace cln { // Subroutine for evaluating // sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n))) // where all the entries are small integers (ideally polynomials in n). -// This is fast because it groups factors together before multiplying. +// Some of the factors (a,b,p,q) may be omitted. They are then understood to +// be 1. This is fast because it groups factors together before multiplying. +// Result will be a cl_LF with len digits. + +// There are various alternative implementations of the same algorithm that +// differ in the way the series is represented and, as a consequence, in memory +// consumption. +// +// 1st implementation (series is precomputed entirely) // Arguments: // Vectors p[0..N-1], q[0..N-1], a[0..N-1], b[0..N-1], N. // Some of the vectors (a,b,p,q) can be a NULL pointer, all of its entries @@ -20,7 +30,24 @@ namespace cln { // split off q[n] into q[n]*2^qs[n]. qs may be NULL, in that case no shift // optimizations will be used. (They are worth it only if a significant // amount of multiplication work can be saved by shifts.) -// Result will be a cl_LF with len digits. +// +// 2nd implemenation (series is computed on demand, as a stream) +// In this alternate implementation the series is not represented as a couple +// of arrays, but as a method returning each tuple (p(n),q(n),a(n),b(n)) +// in turn. This is preferrable if the a(n) are big, in order to avoid too +// much memory usage at the same time. +// The next() function is called N times and is expected to return +// (p(n),q(n),a(n),b(n)) for n=0..N-1 in that order. +// +// 3rd implemenation (series is computed on demand and truncated early) +// This is like the second implementation, but it coerces the integer factors +// to cl_LF of a given length (trunclen) as soon as the integer factor's size +// exceeds the size to store the cl_LF. For this to make sense, trunclen must +// not be smaller than len. In practice, this can shave off substantially from +// the memory consumption but it also bears a potential for rounding errors. +// A minimum trunclen that guarantees correctness must be evaluated on a +// case-by-case basis. + struct cl_rational_series { // To be set explicitly. @@ -42,7 +69,21 @@ struct cl_pqab_series { const cl_I* bv; uintC* qsv; }; +struct cl_pqab_series_term { + cl_I p; + cl_I q; + cl_I a; + cl_I b; +}; +struct cl_pqab_series_stream { + cl_pqab_series_term (*nextfn)(cl_pqab_series_stream&); + cl_pqab_series_term next () { return nextfn(*this); } + // Constructor. + cl_pqab_series_stream (cl_pqab_series_term (*n)(cl_pqab_series_stream&)) : nextfn (n) {} +}; extern const cl_LF eval_rational_series (uintC N, const cl_pqab_series& args, uintC len); +extern const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len); +extern const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len, uintC trunclen); struct cl_pqb_series { const cl_I* pv; @@ -50,7 +91,20 @@ struct cl_pqb_series { const cl_I* bv; uintC* qsv; }; +struct cl_pqb_series_term { + cl_I p; + cl_I q; + cl_I b; +}; +struct cl_pqb_series_stream { + cl_pqb_series_term (*nextfn)(cl_pqb_series_stream&); + cl_pqb_series_term next () { return nextfn(*this); } + // Constructor. + cl_pqb_series_stream (cl_pqb_series_term (*n)(cl_pqb_series_stream&)) : nextfn (n) {} +}; extern const cl_LF eval_rational_series (uintC N, const cl_pqb_series& args, uintC len); +extern const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len); +extern const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len, uintC trunclen); struct cl_pqa_series { const cl_I* pv; @@ -58,14 +112,39 @@ struct cl_pqa_series { const cl_I* av; uintC* qsv; }; +struct cl_pqa_series_term { + cl_I p; + cl_I q; + cl_I a; +}; +struct cl_pqa_series_stream { + cl_pqa_series_term (*nextfn)(cl_pqa_series_stream&); + cl_pqa_series_term next () { return nextfn(*this); } + // Constructor. + cl_pqa_series_stream (cl_pqa_series_term (*n)(cl_pqa_series_stream&)) : nextfn (n) {} +}; extern const cl_LF eval_rational_series (uintC N, const cl_pqa_series& args, uintC len); +extern const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len); +extern const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len, uintC trunclen); struct cl_pq_series { const cl_I* pv; cl_I* qv; uintC* qsv; }; +struct cl_pq_series_term { + cl_I p; + cl_I q; +}; +struct cl_pq_series_stream { + cl_pq_series_term (*nextfn)(cl_pq_series_stream&); + cl_pq_series_term next () { return nextfn(*this); } + // Constructor. + cl_pq_series_stream (cl_pq_series_term (*n)(cl_pq_series_stream&)) : nextfn (n) {} +}; extern const cl_LF eval_rational_series (uintC N, const cl_pq_series& args, uintC len); +extern const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len); +extern const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len, uintC trunclen); struct cl_pab_series { const cl_I* pv; @@ -140,67 +219,6 @@ struct cl__series { extern const cl_LF eval_rational_series (uintC N, const cl__series& args, uintC len); -// In this alternate implementation the series is not represented as a couple -// of arrays, but as a method returning each tuple (p(n),q(n),a(n),b(n)) -// in turn. This is preferrable if the a(n) are big, in order to avoid too -// much memory usage at the same time. -// Some of the factors (a,b) may be omitted. They are then understood to be 1. -// The next function is called N times and is expected to return -// (p(n),q(n),a(n),b(n)) for n=0..N-1 in that order. - -struct cl_pqab_series_term { - cl_I p; - cl_I q; - cl_I a; - cl_I b; -}; -struct cl_pqab_series_stream { - cl_pqab_series_term (*nextfn)(cl_pqab_series_stream&); - cl_pqab_series_term next () { return nextfn(*this); } - // Constructor. - cl_pqab_series_stream (cl_pqab_series_term (*n)(cl_pqab_series_stream&)) : nextfn (n) {} -}; -extern const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len); - -struct cl_pqb_series_term { - cl_I p; - cl_I q; - cl_I b; -}; -struct cl_pqb_series_stream { - cl_pqb_series_term (*nextfn)(cl_pqb_series_stream&); - cl_pqb_series_term next () { return nextfn(*this); } - // Constructor. - cl_pqb_series_stream (cl_pqb_series_term (*n)(cl_pqb_series_stream&)) : nextfn (n) {} -}; -extern const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len); - -struct cl_pqa_series_term { - cl_I p; - cl_I q; - cl_I a; -}; -struct cl_pqa_series_stream { - cl_pqa_series_term (*nextfn)(cl_pqa_series_stream&); - cl_pqa_series_term next () { return nextfn(*this); } - // Constructor. - cl_pqa_series_stream (cl_pqa_series_term (*n)(cl_pqa_series_stream&)) : nextfn (n) {} -}; -extern const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len); - -struct cl_pq_series_term { - cl_I p; - cl_I q; -}; -struct cl_pq_series_stream { - cl_pq_series_term (*nextfn)(cl_pq_series_stream&); - cl_pq_series_term next () { return nextfn(*this); } - // Constructor. - cl_pq_series_stream (cl_pq_series_term (*n)(cl_pq_series_stream&)) : nextfn (n) {} -}; -extern const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len); - - // [Generalization.] // Subroutine: // Evaluates S = sum(N1 <= n < N2, (p(N1)...p(n))/(q(N1)...q(n))) @@ -219,17 +237,28 @@ struct cl_pqcd_series_term { cl_I c; cl_I d; }; +template struct cl_pqcd_series_result { - cl_I P; - cl_I Q; - cl_I T; - cl_I C; - cl_I D; - cl_I V; + cl_T P; + cl_T Q; + cl_T T; + cl_T C; + cl_T D; + cl_T V; +}; +struct cl_pqcd_series_stream { + cl_pqcd_series_term (*nextfn)(cl_pqcd_series_stream&); + cl_pqcd_series_term next () { return nextfn(*this); } + // Constructor. + cl_pqcd_series_stream( cl_pqcd_series_term (*n)(cl_pqcd_series_stream&)) : nextfn (n) {} }; -extern void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_result& Z, bool rightmost = true); +extern void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_result& Z, bool rightmost = true); +extern void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result& Z, bool rightmost = true); +extern void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result& Z, uintC trunclen, bool rightmost = true); // Ditto, but returns U/S. extern const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_term* args, uintC len); +extern const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len); +extern const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len, uintC trunclen); // [Special case c(n)=1.] // Subroutine: @@ -247,13 +276,14 @@ struct cl_pqd_series_term { cl_I q; cl_I d; }; +template struct cl_pqd_series_result { - cl_I P; - cl_I Q; - cl_I T; - cl_I C; - cl_I D; - cl_I V; + cl_T P; + cl_T Q; + cl_T T; + cl_T C; + cl_T D; + cl_T V; }; struct cl_pqd_series_stream { cl_pqd_series_term (*nextfn)(cl_pqd_series_stream&); @@ -261,11 +291,24 @@ struct cl_pqd_series_stream { // 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); +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); +extern void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result& Z, uintC trunclen, 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); +extern const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len, uintC trunclen); + +// Helper function to convert integer of length > trunclen to long float of +// length = trunclen. +inline void +truncate_precision(cl_R& x, uintC trunclen) +{ + if (instanceof(x,cl_I_ring) && + integer_length(the(x))>trunclen*intDsize) { + x = cl_I_to_LF(the(x),trunclen); + } +} } // namespace cln diff --git a/src/float/transcendental/cl_LF_zeta3.cc b/src/float/transcendental/cl_LF_zeta3.cc index 31914d3..1c09773 100644 --- a/src/float/transcendental/cl_LF_zeta3.cc +++ b/src/float/transcendental/cl_LF_zeta3.cc @@ -62,7 +62,7 @@ const cl_LF zeta3 (uintC len) var uintC actuallen = len+2; // 2 Schutz-Digits var uintC N = ceiling(actuallen*intDsize,10); // 1024^-N <= 2^(-intDsize*actuallen). - var cl_LF sum = eval_rational_series(N,series,actuallen); + var cl_LF sum = eval_rational_series(N,series,actuallen,actuallen); return scale_float(shorten(sum,len),-1); } // Bit complexity (N := len): O(log(N)^2*M(N)). diff --git a/src/float/transcendental/cl_LF_zeta_int.cc b/src/float/transcendental/cl_LF_zeta_int.cc index 3f80cdd..503b835 100644 --- a/src/float/transcendental/cl_LF_zeta_int.cc +++ b/src/float/transcendental/cl_LF_zeta_int.cc @@ -91,35 +91,68 @@ const cl_LF compute_zeta_cvz2 (int s, uintC len) // Method: // zeta(s) = 1/(1-2^(1-s)) sum(n=0..infty, (-1)^n/(n+1)^s), // with Cohen-Villegas-Zagier convergence acceleration, and - // evaluated using the binary splitting algorithm. - var uintC actuallen = len+2; // 2 Schutz-Digits + // evaluated using the binary splitting algorithm with truncation. + var uintC actuallen = len+2; // 2 guard digits var uintC N = (uintC)(0.39321985*intDsize*actuallen)+1; - 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) (2*(cl_I)(N-n)*(cl_I)(N+n)); - init1(cl_I, args[n].q) ((cl_I)(2*n+1)*(cl_I)(n+1)); - init1(cl_I, args[n].d) (evenp(n) - ? expt_pos(n+1,s) - : -expt_pos(n+1,s)); - } - var cl_pqd_series_result sums; - eval_pqd_series_aux(N,args,sums); + struct rational_series_stream : cl_pqd_series_stream { + uintC n; + int s; + uintC N; + 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 uintC s = thiss.s; + var uintC N = thiss.N; + var cl_pqd_series_term result; + result.p = 2*(cl_I)(N-n)*(cl_I)(N+n); + result.q = (cl_I)(2*n+1)*(cl_I)(n+1); + result.d = evenp(n) ? expt_pos(n+1,s) : -expt_pos(n+1,s); + thiss.n = n+1; + return result; + } + rational_series_stream (int _s, uintC _N) + : cl_pqd_series_stream (rational_series_stream::computenext), + n (0), s (_s), N (_N) {} + } series(s,N); + var cl_pqd_series_result sums; + eval_pqd_series_aux(N,series,sums,actuallen); // Here we need U/(1+S) = V/D(Q+T). var cl_LF result = cl_I_to_LF(sums.V,actuallen) / The(cl_LF)(sums.D * cl_I_to_LF(sums.Q+sums.T,actuallen)); - for (n = 0; n < N; n++) { - args[n].p.~cl_I(); - args[n].q.~cl_I(); - args[n].d.~cl_I(); - } result = shorten(result,len); // verkürzen und fertig // Zum Schluss mit 2^(s-1)/(2^(s-1)-1) multiplizieren: return scale_float(result,s-1) / (ash(1,s-1)-1); } // Bit complexity (N = len): O(log(N)^2*M(N)). +// Timings of the above algorithm in seconds, on a P-4, 3GHz, running Linux. +// s 5 15 +// N sum_exp sum_cvz1 sum_cvz2 sum_exp sum_cvz1 sum_cvz2 +// 125 0.60 0.04 0.06 1.88 0.04 0.20 +// 250 1.60 0.13 0.19 4.82 0.15 0.58 +// 500 4.3 0.48 0.60 12.2 0.55 1.67 +// 1000 11.0 1.87 1.63 31.7 2.11 4.60 +// 2000 28.0 7.4 4.23 111 8.2 11.3 +// 4000 70.2 30.6 10.6 50 44 +// 8000 142 26.8 169 75 +// asymp. FAST N^2 FAST FAST N^2 FAST +// +// s 35 75 +// N sum_exp sum_cvz1 sum_cvz2 sum_exp sum_cvz1 sum_cvz2 +// 125 4.70 0.05 0.53 11.3 0.07 1.35 +// 250 12.5 0.19 1.62 28.7 0.25 3.74 +// 500 31.3 0.69 4.40 70.2 0.96 10.2 +// 1000 88.8 2.70 11.4 191 3.76 25.4 +// 2000 10.9 28.9 15.6 64.3 +// 4000 46 73 64.4 170 +// 8000 215 178 295 397 +// 16000 898 419 1290 972 +// asymp. FAST N^2 FAST FAST N^2 FAST +// +// The break-even point between cvz1 and cvz2 seems to grow linearly with s. + // Timings of the above algorithm, on an i486 33 MHz, running Linux. // s 5 15 // N sum_exp sum_cvz1 sum_cvz2 sum_exp sum_cvz1 sum_cvz2 @@ -137,8 +170,10 @@ const cl_LF compute_zeta_cvz2 (int s, uintC len) const cl_LF zeta (int s, uintC len) { if (!(s > 1)) - throw runtime_exception(); - if (len < 280*(uintL)s) + throw runtime_exception("zeta(s) with illegal s<2."); + if (s==3) + return zeta3(len); + if (len < 220*(uintC)s) return compute_zeta_cvz1(s,len); else return compute_zeta_cvz2(s,len);