* 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-13 Richard B. Kreckel <kreckel@ginac.de>
+
+ 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 <kreckel@ginac.de>
More memory efficient Euler-Mascheroni constant:
// 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))
? square((cl_I)(2*n+1))
: -square((cl_I)(2*n+1)));
}
- var cl_pqd_series_result sums;
+ var cl_pqd_series_result<cl_I> sums;
eval_pqd_series_aux(N,args,sums);
// Here we need U/(1+S) = V/D(Q+T).
var cl_LF result =
} 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);
}
#include "cl_LF_tran.h"
#include "cl_LF.h"
#include "cln/integer.h"
+#include "cln/real.h"
#include "cl_alloca.h"
namespace cln {
// 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;
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<cl_R> 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)).
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
#include "cln/lfloat.h"
#include "cln/integer.h"
+#include "cln/real.h"
#include "cln/exception.h"
#include "cl_LF.h"
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)).
#include "cln/lfloat.h"
#include "cln/integer.h"
+#include "cln/real.h"
#include "cln/exception.h"
#include "cl_LF.h"
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)).
#include "cln/lfloat.h"
#include "cln/integer.h"
+#include "cln/real.h"
#include "cln/exception.h"
#include "cl_LF.h"
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)).
#include "cln/lfloat.h"
#include "cln/integer.h"
+#include "cln/real.h"
#include "cln/exception.h"
#include "cl_LF.h"
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)).
+++ /dev/null
-// 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
+++ /dev/null
-// 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
+++ /dev/null
-// 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
+++ /dev/null
-// 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
#include "cln/lfloat.h"
#include "cln/integer.h"
+#include "cln/real.h"
#include "cl_LF.h"
namespace cln {
{
if (N==0)
return cl_I_to_LF(0,len);
- var cl_pqcd_series_result sums;
+ var cl_pqcd_series_result<cl_I> 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
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<cl_I> 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<cl_R> 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
// 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<cl_I>& Z, bool rightmost)
{
// N = N2-N1
switch (N) {
default: {
var uintC Nm = N/2; // midpoint
// Compute left part.
- var cl_pqcd_series_result L;
+ var cl_pqcd_series_result<cl_I> L;
eval_pqcd_series_aux(Nm,args+0,L,false);
// Compute right part.
- var cl_pqcd_series_result R;
+ var cl_pqcd_series_result<cl_I> R;
eval_pqcd_series_aux(N-Nm,args+Nm,R,rightmost);
// Put together partial results.
if (!rightmost) { Z.P = L.P * R.P; }
}
}
+void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_I>& 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<cl_I> L;
+ eval_pqcd_series_aux(Nm,args,L,false);
+ // Compute right part.
+ var cl_pqcd_series_result<cl_I> 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<cl_R>& 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<cl_R> L;
+ eval_pqcd_series_aux(Nm,args,L,trunclen,false);
+ // Compute right part.
+ var cl_pqcd_series_result<cl_R> 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
#include "cln/lfloat.h"
#include "cln/integer.h"
+#include "cln/real.h"
#include "cl_LF.h"
namespace cln {
{
if (N==0)
return cl_I_to_LF(0,len);
- var cl_pqd_series_result sums;
+ var cl_pqd_series_result<cl_I> 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
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<cl_I> 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<cl_R> 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
// 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<cl_I>& Z, bool rightmost)
{
// N = N2-N1
switch (N) {
default: {
var uintC Nm = N/2; // midpoint
// Compute left part.
- var cl_pqd_series_result L;
+ var cl_pqd_series_result<cl_I> L;
eval_pqd_series_aux(Nm,args+0,L,false);
// Compute right part.
- var cl_pqd_series_result R;
+ var cl_pqd_series_result<cl_I> R;
eval_pqd_series_aux(N-Nm,args+Nm,R,rightmost);
// Put together partial results.
if (!rightmost) { Z.P = L.P * R.P; }
}
}
+void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_I>& 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<cl_I> L;
+ eval_pqd_series_aux(Nm,args,L,false);
+ // Compute right part.
+ var cl_pqd_series_result<cl_I> 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<cl_R>& 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<cl_R> L;
+ eval_pqd_series_aux(Nm,args,L,trunclen,false);
+ // Compute right part.
+ var cl_pqd_series_result<cl_R> 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
+++ /dev/null
-// 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
+++ /dev/null
-// 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
#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
// 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.
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;
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;
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;
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)))
cl_I c;
cl_I d;
};
+template<class cl_T>
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<cl_I>& Z, bool rightmost = true);
+extern void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_I>& Z, bool rightmost = true);
+extern void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_R>& 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:
cl_I q;
cl_I d;
};
+template<class cl_T>
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&);
// 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<cl_I>& Z, bool rightmost = true);
+extern void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_I>& Z, bool rightmost = true);
+extern void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_R>& 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<cl_I>(x))>trunclen*intDsize) {
+ x = cl_I_to_LF(the<cl_I>(x),trunclen);
+ }
+}
} // namespace cln
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)).
// 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<cl_I> 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
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);
projects / cln.git / commitdiff