+2008-01-11 Richard B. Kreckel <kreckel@ginac.de>
+
+ Make some functions more memory efficient:
+ * src/float/transcendental/cl_LF_tran.h (eval_rational_series): The
+ evaluation of streamed rational series may profit from shift-counting Q,
+ too. Introduce a template parameter to determine whether shift-counting
+ is to be used or not.
+ * src/float/transcendental/cl_LF_ratseries_pqb.cc: Introduce template
+ parameter.
+ * src/float/transcendental/cl_LF_ratseries_pqa.cc: Likewise.
+ * src/float/transcendental/cl_LF_ratseries_pqab.cc: Likewise.
+ * src/float/transcendental/cl_LF_ratseries_qa.cc: Likewise.
+ * src/float/transcendental/cl_LF_ratseries_qab.cc: Likewise.
+ * src/float/transcendental/cl_LF_ratseries_q.cc: Likewise, added
+ overload for streamed expansion.
+ * src/float/transcendental/cl_LF_ratseries_qb.cc: Likewise.
+ * src/float/transcendental/cl_LF_ratseries_pq.cc: Introduce template
+ parameter, added overload for streamed expansion using shift-counts.
+ * src/float/transcendental/cl_LF_zeta3.cc: Adapt to above changes.
+ * src/float/transcendental/cl_LF_pi.cc: Likewise.
+ * src/float/transcendental/cl_LF_eulerconst.cc: Likewise.
+ * src/float/transcendental/cl_LF_catalanconst.cc: Likewise.
+ * src/float/transcendental/cl_LF_cossin_aux.cc: Likewise.
+ * src/float/transcendental/cl_LF_coshsinh_aux.cc: Likewise.
+ * src/float/transcendental/cl_LF_atanh_recip.cc: Use streamed series.
+ * src/float/transcendental/cl_LF_atan_recip.cc: Likewise.
+ * src/float/transcendental/cl_LF_exp1.cc: Likewise.
+ * src/float/transcendental/cl_LF_exp_aux.cc: Likewise.
+ * src/float/transcendental/cl_LF_ratseries.cc: Removed.
+
2008-01-06 Ralf Wildenhues <Ralf.Wildenhues@gmx.de>
Richard B. Kreckel <kreckel@ginac.de>
#include "cln/lfloat.h"
#include "cl_LF.h"
#include "cl_LF_tran.h"
-#include "cl_alloca.h"
#undef floor
#include <cmath>
var uintC actuallen = len + 1;
var cl_I m2 = m*m+1;
var uintC N = (uintC)(0.69314718*intDsize*actuallen/::log(double_approx(m2))) + 1;
- CL_ALLOCA_STACK;
- var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var uintC n;
- new (&pv[0]) cl_I (m);
- new (&qv[0]) cl_I (m2);
- for (n = 1; n < N; n++) {
- new (&pv[n]) cl_I (2*n);
- new (&qv[n]) cl_I ((2*n+1)*m2);
- }
- var cl_pq_series series;
- series.pv = pv; series.qv = qv; series.qsv = NULL;
- var cl_LF result = eval_rational_series(N,series,actuallen);
- for (n = 0; n < N; n++) {
- pv[n].~cl_I();
- qv[n].~cl_I();
- }
+ struct rational_series_stream : cl_pq_series_stream {
+ var uintC n;
+ var cl_I m;
+ var cl_I m2;
+ static cl_pq_series_term computenext (cl_pq_series_stream& thisss)
+ {
+ var rational_series_stream& thiss = (rational_series_stream&)thisss;
+ var uintC n = thiss.n;
+ var cl_pq_series_term result;
+ if (n==0) {
+ result.p = thiss.m;
+ result.q = thiss.m2;
+ } else {
+ result.p = 2*n;
+ result.q = (2*n+1)*thiss.m2;
+ }
+ thiss.n = n+1;
+ return result;
+ }
+ rational_series_stream(const cl_I& m_, const cl_I& m2_)
+ : cl_pq_series_stream (rational_series_stream::computenext),
+ n(0), m(m_), m2(m2_) {}
+ } series(m,m2);
+ var cl_LF result = eval_rational_series<false>(N,series,actuallen);
return shorten(result,len);
}
// Bit complexity (N = len): O(log(N)^2*M(N)).
#include "cln/lfloat.h"
#include "cl_LF.h"
#include "cl_LF_tran.h"
-#include "cl_alloca.h"
#undef floor
#include <cmath>
const cl_LF cl_atanh_recip (cl_I m, uintC len)
{
var uintC actuallen = len + 1;
- var cl_I m2 = m*m;
var uintC N = (uintC)(0.69314718*intDsize/2*actuallen/::log(double_approx(m))) + 1;
- CL_ALLOCA_STACK;
- var cl_I* bv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var uintC n;
- for (n = 0; n < N; n++) {
- new (&bv[n]) cl_I ((cl_I)(2*n+1));
- new (&qv[n]) cl_I (n==0 ? m : m2);
- }
- var cl_qb_series series;
- series.bv = bv;
- series.qv = qv; series.qsv = NULL;
- var cl_LF result = eval_rational_series(N,series,actuallen);
- for (n = 0; n < N; n++) {
- bv[n].~cl_I();
- qv[n].~cl_I();
- }
+ struct rational_series_stream : cl_qb_series_stream {
+ var uintC n;
+ var cl_I m;
+ var cl_I m2;
+ static cl_qb_series_term computenext (cl_qb_series_stream& thisss)
+ {
+ var rational_series_stream& thiss = (rational_series_stream&)thisss;
+ var uintC n = thiss.n;
+ var cl_qb_series_term result;
+ result.b = 2*n+1;
+ result.q = (n==0 ? thiss.m : thiss.m2);
+ thiss.n = n+1;
+ return result;
+ }
+ rational_series_stream(const cl_I& m_)
+ : cl_qb_series_stream (rational_series_stream::computenext),
+ n(0), m(m_), m2(square(m_)) {}
+ } series(m);
+ var cl_LF result = eval_rational_series<false>(N,series,actuallen);
return shorten(result,len);
}
// Bit complexity (N = len): O(log(N)^2*M(N)).
// Every summand gives 0.6 new decimal digits in precision.
// The sum is best evaluated using fixed-point arithmetic,
// so that the precision is reduced for the later summands.
- var uintC actuallen = len + 2; // 2 Schutz-Digits
+ var uintC actuallen = len + 2; // 2 guard digits
var sintC scale = intDsize*actuallen;
var cl_I sum = 0;
var cl_I n = 0;
// 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,actuallen);
+ var cl_LF fsum = eval_rational_series<false>(N,series,actuallen,actuallen);
var cl_LF g =
scale_float(The(cl_LF)(3*fsum)
+ The(cl_LF)(pi(actuallen))
const cl_LF compute_catalanconst_expintegral1 (uintC len)
{
// We compute f(x) classically and g(x) using the partial sums of f(x).
- var uintC actuallen = len+2; // 2 Schutz-Digits
+ var uintC actuallen = len+2; // 2 guard digits
var uintC x = (uintC)(0.693148*intDsize*actuallen)+1;
var uintC N = (uintC)(2.718281828*x);
var cl_LF fterm = cl_I_to_LF(1,actuallen);
// the sums.
const cl_LF compute_catalanconst_expintegral2 (uintC len)
{
- var uintC actuallen = len+2; // 2 Schutz-Digits
+ var uintC actuallen = len+2; // 2 guard digits
var uintC x = (uintC)(0.693148*intDsize*actuallen)+1;
var uintC N = (uintC)(2.718281828*x);
CL_ALLOCA_STACK;
// Using Cohen-Villegas-Zagier acceleration, but without binary splitting.
const cl_LF compute_catalanconst_cvz1 (uintC len)
{
- var uintC actuallen = len+2; // 2 Schutz-Digits
+ var uintC actuallen = len+2; // 2 guard digits
var uintC N = (uintC)(0.39321985*intDsize*actuallen)+1;
#if 0
var cl_LF fterm = cl_I_to_LF(2*(cl_I)N*(cl_I)N,actuallen);
// Using Cohen-Villegas-Zagier acceleration, with binary splitting.
const cl_LF compute_catalanconst_cvz2 (uintC len)
{
- var uintC actuallen = len+2; // 2 Schutz-Digits
+ 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));
} 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,actuallen);
+ var cl_LF fsum = eval_rational_series<false>(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);
}
// a(n) = 1, b(n) = 1,
// p(0) = p, p(n) = p^2 for n>0,
// q(0) = 2^lq, q(n) = (2*n)*(2*n+1)*(2^lq)^2 for n>0.
- var uintC actuallen = len+1; // 1 Schutz-Digit
- // How many terms to we need for M bits of precision? N/2 terms suffice,
+ var uintC actuallen = len+1; // 1 guard digit
+ // How many terms do we need for M bits of precision? N/2 terms suffice,
// provided that
// 1/(2^(N*lp)*N!) < 2^-M
// <== N*(log(N)-1)+N*lp*log(2) > M*log(2)
CL_ALLOCA_STACK;
var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
var uintC n;
var cl_I p2 = square(p);
var cl_LF sinhsum;
init1(cl_I, qv[n]) (((cl_I)n*(cl_I)(2*n+1)) << (2*lq+1));
}
var cl_pq_series series;
- series.pv = pv; series.qv = qv; series.qsv = qsv;
- sinhsum = eval_rational_series(N,series,actuallen);
+ series.pv = pv; series.qv = qv;
+ sinhsum = eval_rational_series<true>(N,series,actuallen);
for (n = 0; n < N; n++) {
pv[n].~cl_I();
qv[n].~cl_I();
init1(cl_I, qv[n]) (((cl_I)n*(cl_I)(2*n-1)) << (2*lq+1));
}
var cl_pq_series series;
- series.pv = pv; series.qv = qv; series.qsv = qsv;
- coshsum = eval_rational_series(N,series,actuallen);
+ series.pv = pv; series.qv = qv;
+ coshsum = eval_rational_series<true>(N,series,actuallen);
for (n = 0; n < N; n++) {
pv[n].~cl_I();
qv[n].~cl_I();
// a(n) = 1, b(n) = 1,
// p(0) = p, p(n) = -p^2 for n>0,
// q(0) = 2^lq, q(n) = (2*n)*(2*n+1)*(2^lq)^2 for n>0.
- var uintC actuallen = len+1; // 1 Schutz-Digit
- // How many terms to we need for M bits of precision? N/2 terms suffice,
+ var uintC actuallen = len+1; // 1 guard digit
+ // How many terms do we need for M bits of precision? N/2 terms suffice,
// provided that
// 1/(2^(N*lp)*N!) < 2^-M
// <== N*(log(N)-1)+N*lp*log(2) > M*log(2)
CL_ALLOCA_STACK;
var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
var uintC n;
var cl_I p2 = -square(p);
var cl_LF sinsum;
init1(cl_I, qv[n]) (((cl_I)n*(cl_I)(2*n+1)) << (2*lq+1));
}
var cl_pq_series series;
- series.pv = pv; series.qv = qv; series.qsv = qsv;
- sinsum = eval_rational_series(N,series,actuallen);
+ series.pv = pv; series.qv = qv;
+ sinsum = eval_rational_series<true>(N,series,actuallen);
for (n = 0; n < N; n++) {
pv[n].~cl_I();
qv[n].~cl_I();
init1(cl_I, qv[n]) (((cl_I)n*(cl_I)(2*n-1)) << (2*lq+1));
}
var cl_pq_series series;
- series.pv = pv; series.qv = qv; series.qsv = qsv;
- cossum = eval_rational_series(N,series,actuallen);
+ series.pv = pv; series.qv = qv;
+ cossum = eval_rational_series<true>(N,series,actuallen);
for (n = 0; n < N; n++) {
pv[n].~cl_I();
qv[n].~cl_I();
// If we computed this with floating-point numbers, we would have
// to more than double the floating-point precision because of the large
// extinction which takes place. But luckily we compute with integers.
- var uintC actuallen = len+1; // 1 Schutz-Digit
+ var uintC actuallen = len+1; // 1 guard digit
var uintC z = (uintC)(0.693148*intDsize*actuallen)+1;
var uintC N = (uintC)(3.591121477*z);
CL_ALLOCA_STACK;
}
var cl_pqb_series series;
series.bv = bv;
- series.pv = pv; series.qv = qv; series.qsv = NULL;
- var cl_LF fsum = eval_rational_series(N,series,actuallen);
+ series.pv = pv; series.qv = qv;
+ var cl_LF fsum = eval_rational_series<false>(N,series,actuallen);
for (n = 0; n < N; n++) {
bv[n].~cl_I();
pv[n].~cl_I();
// Finally we compute the sums of the series f(x) and g(x) with N terms
// each.
// We compute f(x) classically and g(x) using the partial sums of f(x).
- var uintC actuallen = len+2; // 2 Schutz-Digits
+ var uintC actuallen = len+2; // 2 guard digits
var uintC x = (uintC)(0.693148*intDsize*actuallen)+1;
var uintC N = (uintC)(2.718281828*x);
var cl_LF one = cl_I_to_LF(1,actuallen);
// the sums.
const cl_LF compute_eulerconst_expintegral2 (uintC len)
{
- var uintC actuallen = len+2; // 2 Schutz-Digits
+ var uintC actuallen = len+2; // 2 guard digits
var uintC x = (uintC)(0.693148*intDsize*actuallen)+1;
var uintC N = (uintC)(2.718281828*x);
CL_ALLOCA_STACK;
const cl_LF compute_eulerconst_besselintegral1 (uintC len)
{
// We compute f(x) classically and g(x) using the partial sums of f(x).
- var uintC actuallen = len+1; // 1 Schutz-Digit
+ var uintC actuallen = len+1; // 1 guard digit
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);
// WARNING: The memory used by this algorithm grown quadratically in N.
// (Because HD_n grows like exp(n), hence HN_n grows like exp(n) as
// well, and we store all HN_n values in an array!)
- var uintC actuallen = len+1; // 1 Schutz-Digit
+ var uintC actuallen = len+1; // 1 guard digit
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);
init1(cl_I, qv[n]) ((cl_I)n*(cl_I)n);
}
var cl_pq_series fseries;
- fseries.pv = pv; fseries.qv = qv; fseries.qsv = NULL;
- var cl_LF fsum = eval_rational_series(N,fseries,actuallen);
+ fseries.pv = pv; fseries.qv = qv;
+ var cl_LF fsum = eval_rational_series<false>(N,fseries,actuallen);
for (n = 0; n < N; n++) {
pv[n].~cl_I();
qv[n].~cl_I();
}
var cl_pqa_series gseries;
gseries.av = av;
- gseries.pv = pv; gseries.qv = qv; gseries.qsv = NULL;
- var cl_LF gsum = eval_rational_series(N,gseries,actuallen);
+ gseries.pv = pv; gseries.qv = qv;
+ var cl_LF gsum = eval_rational_series<false>(N,gseries,actuallen);
for (n = 0; n < N; n++) {
av[n].~cl_I();
pv[n].~cl_I();
};
const cl_LF compute_eulerconst_besselintegral3 (uintC len)
{
- var uintC actuallen = len+1; // 1 Schutz-Digit
+ var uintC actuallen = len+1; // 1 guard digit
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);
init1(cl_I, qv[n]) ((cl_I)n*(cl_I)n);
}
var cl_pq_series fseries;
- fseries.pv = pv; fseries.qv = qv; fseries.qsv = NULL;
- var cl_LF fsum = eval_rational_series(N,fseries,actuallen);
+ fseries.pv = pv; fseries.qv = qv;
+ var cl_LF fsum = eval_rational_series<false>(N,fseries,actuallen);
for (n = 0; n < N; n++) {
pv[n].~cl_I();
qv[n].~cl_I();
}
// Evaluate g(x).
var cl_rational_series_for_g gseries = cl_rational_series_for_g(x);
- var cl_LF gsum = eval_rational_series(N,gseries,actuallen);
+ var cl_LF gsum = eval_rational_series<false>(N,gseries,actuallen);
var cl_LF result = gsum/fsum - ln(cl_I_to_LF(sx,actuallen));
return shorten(result,len); // verkürzen und fertig
}
#include "cl_LF_tran.h"
#include "cl_LF.h"
#include "cln/integer.h"
-#include "cl_alloca.h"
#undef floor
#include <cmath>
// Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
// with appropriate N, and
// a(n) = 1, b(n) = 1, p(n) = 1, q(n) = n for n>0.
- var uintC actuallen = len+1; // 1 Schutz-Digit
- // How many terms to we need for M bits of precision? N terms suffice,
+ var uintC actuallen = len+1; // 1 guard digit
+ // How many terms do we need for M bits of precision? N terms suffice,
// provided that
// 1/N! < 2^-M
// <== N*(log(N)-1) > M*log(2)
var uintC N1 = (uintC)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0));
var uintC N2 = (uintC)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0))+1;
var uintC N = N2+2;
- CL_ALLOCA_STACK;
- var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
- var uintC n;
- for (n = 0; n < N; n++) {
- init1(cl_I, qv[n]) (n==0 ? 1 : n);
- }
- var cl_q_series series;
- series.qv = qv;
- series.qsv = (len >= 1000 && len <= 10000 ? qsv : 0); // tiny speedup
- var cl_LF fsum = eval_rational_series(N,series,actuallen);
- for (n = 0; n < N; n++) {
- qv[n].~cl_I();
- }
+ struct rational_series_stream : cl_q_series_stream {
+ var uintC n;
+ static cl_q_series_term computenext (cl_q_series_stream& thisss)
+ {
+ var rational_series_stream& thiss = (rational_series_stream&)thisss;
+ var uintC n = thiss.n;
+ var cl_q_series_term result;
+ result.q = (n==0 ? 1 : n);
+ thiss.n = n+1;
+ return result;
+ }
+ rational_series_stream()
+ : cl_q_series_stream (rational_series_stream::computenext),
+ n(0) {}
+ } series;
+ var cl_LF fsum = eval_rational_series<false>(N,series,actuallen);
return shorten(fsum,len); // verkürzen und fertig
}
// Bit complexity (N = len): O(log(N)*M(N)).
#include "cl_LF_tran.h"
#include "cl_LF.h"
#include "cln/integer.h"
-#include "cl_alloca.h"
#include "cln/exception.h"
#undef floor
// Evaluate a sum(0 <= n < N, a(n)/b(n) * (p(0)...p(n))/(q(0)...q(n)))
// with appropriate N, and
// a(n) = 1, b(n) = 1, p(n) = p for n>0, q(n) = n*2^lq for n>0.
- var uintC actuallen = len+1; // 1 Schutz-Digit
- // How many terms to we need for M bits of precision? N terms suffice,
+ var uintC actuallen = len+1; // 1 guard digit
+ // How many terms do we need for M bits of precision? N terms suffice,
// provided that
// 1/(2^(N*lp)*N!) < 2^-M
// <== N*(log(N)-1)+N*lp*log(2) > M*log(2)
var uintC N1 = (uintC)(0.693147*intDsize*actuallen/(::log((double)N0)-1.0+0.693148*lp));
var uintC N2 = (uintC)(0.693148*intDsize*actuallen/(::log((double)N1)-1.0+0.693147*lp))+1;
var uintC N = N2+2;
- CL_ALLOCA_STACK;
- var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
- var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
- var uintC n;
- init1(cl_I, pv[0]) (1);
- init1(cl_I, qv[0]) (1);
- for (n = 1; n < N; n++) {
- init1(cl_I, pv[n]) (p);
- init1(cl_I, qv[n]) ((cl_I)n << lq);
- }
- var cl_pq_series series;
- series.pv = pv; series.qv = qv; series.qsv = qsv;
- var cl_LF fsum = eval_rational_series(N,series,actuallen);
- for (n = 0; n < N; n++) {
- pv[n].~cl_I();
- qv[n].~cl_I();
- }
+ struct rational_series_stream : cl_pq_series_stream {
+ var uintC n;
+ var cl_I p;
+ var uintE lq;
+ static cl_pq_series_term computenext (cl_pq_series_stream& thisss)
+ {
+ var rational_series_stream& thiss = (rational_series_stream&)thisss;
+ var uintC n = thiss.n;
+ var cl_pq_series_term result;
+ if (n==0) {
+ result.p = 1;
+ result.q = 1;
+ } else {
+ result.p = thiss.p;
+ result.q = (cl_I)n << thiss.lq;
+ }
+ thiss.n = n+1;
+ return result;
+ }
+ rational_series_stream(const cl_I& p_, uintE lq_)
+ : cl_pq_series_stream (rational_series_stream::computenext),
+ n (0), p(p_), lq(lq_) {}
+ } series(p, lq);
+ var cl_LF fsum = eval_rational_series<true>(N,series,actuallen);
return shorten(fsum,len); // verkürzen und fertig
}}
// Bit complexity (N = len, and if p has length O(log N) and ql = O(log N)):
// )
// (/ (expt a 2) t)
// )
- var uintC actuallen = len + 1; // 1 Schutz-Digit
+ var uintC actuallen = len + 1; // 1 guard digit
var uintE uexp_limit = LF_exp_mid - intDsize*len;
// Ein Long-Float ist genau dann betragsmäßig <2^-n, wenn
// sein Exponent < LF_exp_mid-n = uexp_limit ist.
// = 2^(k+1) * [wa_k^4 - ((wa_k^2+wb_k^2)/2)^2].
// Hence,
// pi = AGM(1,1/sqrt(2))^2 * 1/(1/2 - sum(k even, 2^k*[....])).
- var uintC actuallen = len + 1; // 1 Schutz-Digit
+ var uintC actuallen = len + 1; // 1 guard digit
var uintE uexp_limit = LF_exp_mid - intDsize*len;
var cl_LF one = cl_I_to_LF(1,actuallen);
var cl_LF a = one;
// in precision.
// The sum is best evaluated using fixed-point arithmetic,
// so that the precision is reduced for the later summands.
- var uintC actuallen = len + 4; // 4 Schutz-Digits
+ var uintC actuallen = len + 4; // 4 guard digits
var sintC scale = intDsize*actuallen;
static const cl_I A = "163096908";
static const cl_I B = "6541681608";
: cl_pqa_series_stream (rational_series_stream::computenext),
n (0) {}
} series;
- var uintC actuallen = len + 4; // 4 Schutz-Digits
+ var uintC actuallen = len + 4; // 4 guard digits
static const cl_I A = "163096908";
static const cl_I B = "6541681608";
static const cl_I J1 = "10939058860032000"; // 72*abs(J)
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,actuallen);
+ var cl_LF fsum = eval_rational_series<false>(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
-// eval_rational_series().
+// eval_rational_series<bool>().
// General includes.
#include "cl_sysdep.h"
#include "cln/real.h"
#include "cln/exception.h"
#include "cl_LF.h"
+#include "cl_alloca.h"
namespace cln {
}
}
+template<>
+const cl_LF eval_rational_series<false> (uintC N, const cl_pq_series& 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_pqs_series_aux (uintC N1, uintC N2,
- const cl_pq_series& args,
+ const cl_pq_series& args, const uintC* qsv,
cl_I* P, cl_I* Q, uintC* QS, cl_I* T)
{
switch (N2 - N1) {
case 1:
if (P) { *P = args.pv[N1]; }
*Q = args.qv[N1];
- *QS = args.qsv[N1];
+ *QS = qsv[N1];
*T = args.pv[N1];
break;
case 2: {
var cl_I p01 = args.pv[N1] * args.pv[N1+1];
if (P) { *P = p01; }
*Q = args.qv[N1] * args.qv[N1+1];
- *QS = args.qsv[N1] + args.qsv[N1+1];
- *T = ((args.qv[N1+1] * args.pv[N1]) << args.qsv[N1+1])
+ *QS = qsv[N1] + qsv[N1+1];
+ *T = ((args.qv[N1+1] * args.pv[N1]) << qsv[N1+1])
+ p01;
break;
}
if (P) { *P = p012; }
var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
*Q = args.qv[N1] * q12;
- *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2];
- *T = ((q12 * args.pv[N1]) << (args.qsv[N1+1] + args.qsv[N1+2]))
- + ((args.qv[N1+2] * p01) << args.qsv[N1+2])
+ *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2];
+ *T = ((q12 * args.pv[N1]) << (qsv[N1+1] + qsv[N1+2]))
+ + ((args.qv[N1+2] * p01) << qsv[N1+2])
+ p012;
break;
}
var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
var cl_I q123 = args.qv[N1+1] * q23;
*Q = args.qv[N1] * q123;
- *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2] + args.qsv[N1+3];
- *T = ((((((q123 * args.pv[N1]) << args.qsv[N1+1])
- + q23 * p01) << args.qsv[N1+2])
- + args.qv[N1+3] * p012) << args.qsv[N1+3])
+ *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2] + qsv[N1+3];
+ *T = ((((((q123 * args.pv[N1]) << qsv[N1+1])
+ + q23 * p01) << qsv[N1+2])
+ + args.qv[N1+3] * p012) << qsv[N1+3])
+ p0123;
break;
}
// Compute left part.
var cl_I LP, LQ, LT;
var uintC LQS;
- eval_pqs_series_aux(N1,Nm,args,&LP,&LQ,&LQS,<);
+ eval_pqs_series_aux(N1,Nm,args,qsv,&LP,&LQ,&LQS,<);
// Compute right part.
var cl_I RP, RQ, RT;
var uintC RQS;
- eval_pqs_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
+ eval_pqs_series_aux(Nm,N2,args,qsv,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
// Put together partial results.
if (P) { *P = LP*RP; }
*Q = LQ*RQ;
}
}
-const cl_LF eval_rational_series (uintC N, const cl_pq_series& args, uintC len)
+template<>
+const cl_LF eval_rational_series<true> (uintC N, const cl_pq_series& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
var cl_I Q, T;
- if (!args.qsv) {
- eval_pq_series_aux(0,N,args,NULL,&Q,&T);
- return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
- } else {
- // Precomputation of the shift counts:
- // Split qv[n] into qv[n]*2^qsv[n].
- {
- var cl_I* qp = args.qv;
- var uintC* qsp = args.qsv;
- for (var uintC n = 0; n < N; n++, qp++, qsp++) {
- // Pull out maximal power of 2 out of *qp = args.qv[n].
- var uintC qs = 0;
- if (!zerop(*qp)) {
- qs = ord2(*qp);
- if (qs > 0)
- *qp = *qp >> qs;
- }
- *qsp = qs;
- }
- }
- // Main computation.
- var uintC QS;
- eval_pqs_series_aux(0,N,args,NULL,&Q,&QS,&T);
- return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS);
+ // Precomputation of the shift counts:
+ // Split qv[n] into qv[n]*2^qsv[n].
+ CL_ALLOCA_STACK;
+ var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
+ var cl_I* qp = args.qv;
+ var uintC* qsp = qsv;
+ for (var uintC n = 0; n < N; n++, qp++, qsp++) {
+ *qsp = pullout_shiftcount(*qp);
}
+ // Main computation.
+ var uintC QS;
+ eval_pqs_series_aux(0,N,args,qsv,NULL,&Q,&QS,&T);
+ 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,
}
}
-const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, cl_pq_series_stream& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
}
+static void eval_pqs_series_aux (uintC N1, uintC N2,
+ cl_pq_series_stream& args,
+ cl_I* P, cl_I* Q, uintC* QS, cl_I* T)
+{
+ switch (N2 - N1) {
+ case 0:
+ throw runtime_exception(); break;
+ case 1: {
+ var cl_pq_series_term v0 = args.next(); // [N1]
+ var uintC qs0 = pullout_shiftcount(v0.q);
+ if (P) { *P = v0.p; }
+ *Q = v0.q;
+ *QS = qs0;
+ *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 uintC qs0 = pullout_shiftcount(v0.q);
+ var uintC qs1 = pullout_shiftcount(v1.q);
+ var cl_I p01 = v0.p * v1.p;
+ if (P) { *P = p01; }
+ *Q = v0.q * v1.q;
+ *QS = qs0 + qs1;
+ *T = ((v1.q * v0.p) << qs1)
+ + 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 uintC qs0 = pullout_shiftcount(v0.q);
+ var uintC qs1 = pullout_shiftcount(v1.q);
+ var uintC qs2 = pullout_shiftcount(v2.q);
+ 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;
+ *QS = qs0 + qs1 + qs2;
+ *T = ((q12 * v0.p) << (qs1 + qs2))
+ + ((v2.q * p01) << qs2)
+ + 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 uintC qs0 = pullout_shiftcount(v0.q);
+ var uintC qs1 = pullout_shiftcount(v1.q);
+ var uintC qs2 = pullout_shiftcount(v2.q);
+ var uintC qs3 = pullout_shiftcount(v3.q);
+ 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;
+ *QS = qs0 + qs1 + qs2 + qs3;
+ *T = ((((((q123 * v0.p) << qs1)
+ + q23 * p01) << qs2)
+ + v3.q * p012) << qs3)
+ + p0123;
+ break;
+ }
+ default: {
+ var uintC Nm = (N1+N2)/2; // midpoint
+ // Compute left part.
+ var cl_I LP, LQ, LT;
+ var uintC LQS;
+ eval_pqs_series_aux(N1,Nm,args,&LP,&LQ,&LQS,<);
+ // Compute right part.
+ var cl_I RP, RQ, RT;
+ var uintC RQS;
+ eval_pqs_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
+ // Put together partial results.
+ if (P) { *P = LP*RP; }
+ *Q = LQ*RQ;
+ *QS = LQS+RQS;
+ // S = LS + LP/LQ * RS, so T = RQ*LT + LP*RT.
+ *T = ((RQ*LT) << RQS) + LP*RT;
+ break;
+ }
+ }
+}
+
+template<>
+const cl_LF eval_rational_series<true> (uintC N, cl_pq_series_stream& args, uintC len)
+{
+ if (N==0)
+ return cl_I_to_LF(0,len);
+ var cl_I Q, T;
+ var uintC QS;
+ eval_pqs_series_aux(0,N,args,NULL,&Q,&QS,&T);
+ 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_R* P, cl_R* Q, cl_R* T,
}
}
-const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len, uintC trunclen)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, cl_pq_series_stream& args, uintC len, uintC trunclen)
{
if (N==0)
return cl_I_to_LF(0,len);
-// eval_rational_series().
+// eval_rational_series<bool>().
// General includes.
#include "cl_sysdep.h"
#include "cln/real.h"
#include "cln/exception.h"
#include "cl_LF.h"
+#include "cl_alloca.h"
namespace cln {
}
}
+template<>
+const cl_LF eval_rational_series<false> (uintC N, const cl_pqa_series& 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_pqsa_series_aux (uintC N1, uintC N2,
- const cl_pqa_series& args,
+ const cl_pqa_series& args, const uintC* qsv,
cl_I* P, cl_I* Q, uintC* QS, cl_I* T)
{
switch (N2 - N1) {
case 1:
if (P) { *P = args.pv[N1]; }
*Q = args.qv[N1];
- *QS = args.qsv[N1];
+ *QS = qsv[N1];
*T = args.av[N1] * args.pv[N1];
break;
case 2: {
var cl_I p01 = args.pv[N1] * args.pv[N1+1];
if (P) { *P = p01; }
*Q = args.qv[N1] * args.qv[N1+1];
- *QS = args.qsv[N1] + args.qsv[N1+1];
- *T = ((args.qv[N1+1] * args.av[N1] * args.pv[N1]) << args.qsv[N1+1])
+ *QS = qsv[N1] + qsv[N1+1];
+ *T = ((args.qv[N1+1] * args.av[N1] * args.pv[N1]) << qsv[N1+1])
+ args.av[N1+1] * p01;
break;
}
if (P) { *P = p012; }
var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
*Q = args.qv[N1] * q12;
- *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2];
- *T = ((q12 * args.av[N1] * args.pv[N1]) << (args.qsv[N1+1] + args.qsv[N1+2]))
- + ((args.qv[N1+2] * args.av[N1+1] * p01) << args.qsv[N1+2])
+ *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2];
+ *T = ((q12 * args.av[N1] * args.pv[N1]) << (qsv[N1+1] + qsv[N1+2]))
+ + ((args.qv[N1+2] * args.av[N1+1] * p01) << qsv[N1+2])
+ args.av[N1+2] * p012;
break;
}
var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
var cl_I q123 = args.qv[N1+1] * q23;
*Q = args.qv[N1] * q123;
- *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2] + args.qsv[N1+3];
- *T = ((((((q123 * args.av[N1] * args.pv[N1]) << args.qsv[N1+1])
- + q23 * args.av[N1+1] * p01) << args.qsv[N1+2])
- + args.qv[N1+3] * args.av[N1+2] * p012) << args.qsv[N1+3])
+ *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2] + qsv[N1+3];
+ *T = ((((((q123 * args.av[N1] * args.pv[N1]) << qsv[N1+1])
+ + q23 * args.av[N1+1] * p01) << qsv[N1+2])
+ + args.qv[N1+3] * args.av[N1+2] * p012) << qsv[N1+3])
+ args.av[N1+3] * p0123;
break;
}
// Compute left part.
var cl_I LP, LQ, LT;
var uintC LQS;
- eval_pqsa_series_aux(N1,Nm,args,&LP,&LQ,&LQS,<);
+ eval_pqsa_series_aux(N1,Nm,args,qsv,&LP,&LQ,&LQS,<);
// Compute right part.
var cl_I RP, RQ, RT;
var uintC RQS;
- eval_pqsa_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
+ eval_pqsa_series_aux(Nm,N2,args,qsv,(P?&RP:(cl_I*)0),&RQ,&RQS,&RT);
// Put together partial results.
if (P) { *P = LP*RP; }
*Q = LQ*RQ;
}
}
-const cl_LF eval_rational_series (uintC N, const cl_pqa_series& args, uintC len)
+template<>
+const cl_LF eval_rational_series<true> (uintC N, const cl_pqa_series& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
var cl_I Q, T;
- if (!args.qsv) {
- eval_pqa_series_aux(0,N,args,NULL,&Q,&T);
- return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
- } else {
- // Precomputation of the shift counts:
- // Split qv[n] into qv[n]*2^qsv[n].
- {
- var cl_I* qp = args.qv;
- var uintC* qsp = args.qsv;
- for (var uintC n = 0; n < N; n++, qp++, qsp++) {
- // Pull out maximal power of 2 out of *qp = args.qv[n].
- var uintC qs = 0;
- if (!zerop(*qp)) {
- qs = ord2(*qp);
- if (qs > 0)
- *qp = *qp >> qs;
- }
- *qsp = qs;
- }
- }
- // Main computation.
- var uintC QS;
- eval_pqsa_series_aux(0,N,args,NULL,&Q,&QS,&T);
- return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(Q,len),QS);
+ // Precomputation of the shift counts:
+ // Split qv[n] into qv[n]*2^qsv[n].
+ CL_ALLOCA_STACK;
+ var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
+ var cl_I* qp = args.qv;
+ var uintC* qsp = qsv;
+ for (var uintC n = 0; n < N; n++, qp++, qsp++) {
+ *qsp = pullout_shiftcount(*qp);
}
+ // Main computation.
+ var uintC QS;
+ eval_pqsa_series_aux(0,N,args,qsv,NULL,&Q,&QS,&T);
+ 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,
}
}
-const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, cl_pqa_series_stream& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
}
}
-const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len, uintC trunclen)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, cl_pqa_series_stream& args, uintC len, uintC trunclen)
{
if (N==0)
return cl_I_to_LF(0,len);
-// eval_rational_series().
+// eval_rational_series<bool>().
// General includes.
#include "cl_sysdep.h"
#include "cln/real.h"
#include "cln/exception.h"
#include "cl_LF.h"
+#include "cl_alloca.h"
namespace cln {
}
}
+template<>
+const cl_LF eval_rational_series<false> (uintC N, const cl_pqab_series& 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_pqsab_series_aux (uintC N1, uintC N2,
- const cl_pqab_series& args,
+ const cl_pqab_series& args, const uintC* qsv,
cl_I* P, cl_I* Q, uintC* QS, cl_I* B, cl_I* T)
{
switch (N2 - N1) {
case 1:
if (P) { *P = args.pv[N1]; }
*Q = args.qv[N1];
- *QS = args.qsv[N1];
+ *QS = qsv[N1];
*B = args.bv[N1];
*T = args.av[N1] * args.pv[N1];
break;
var cl_I p01 = args.pv[N1] * args.pv[N1+1];
if (P) { *P = p01; }
*Q = args.qv[N1] * args.qv[N1+1];
- *QS = args.qsv[N1] + args.qsv[N1+1];
+ *QS = qsv[N1] + qsv[N1+1];
*B = args.bv[N1] * args.bv[N1+1];
- *T = ((args.bv[N1+1] * args.qv[N1+1] * args.av[N1] * args.pv[N1]) << args.qsv[N1+1])
+ *T = ((args.bv[N1+1] * args.qv[N1+1] * args.av[N1] * args.pv[N1]) << qsv[N1+1])
+ args.bv[N1] * args.av[N1+1] * p01;
break;
}
if (P) { *P = p012; }
var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
*Q = args.qv[N1] * q12;
- *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2];
+ *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2];
var cl_I b12 = args.bv[N1+1] * args.bv[N1+2];
*B = args.bv[N1] * b12;
- *T = ((b12 * q12 * args.av[N1] * args.pv[N1]) << (args.qsv[N1+1] + args.qsv[N1+2]))
- + args.bv[N1] * (((args.bv[N1+2] * args.qv[N1+2] * args.av[N1+1] * p01) << args.qsv[N1+2])
+ *T = ((b12 * q12 * args.av[N1] * args.pv[N1]) << (qsv[N1+1] + qsv[N1+2]))
+ + args.bv[N1] * (((args.bv[N1+2] * args.qv[N1+2] * args.av[N1+1] * p01) << qsv[N1+2])
+ args.bv[N1+1] * args.av[N1+2] * p012);
break;
}
var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
var cl_I q123 = args.qv[N1+1] * q23;
*Q = args.qv[N1] * q123;
- *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2] + args.qsv[N1+3];
+ *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2] + qsv[N1+3];
var cl_I b01 = args.bv[N1] * args.bv[N1+1];
var cl_I b23 = args.bv[N1+2] * args.bv[N1+3];
*B = b01 * b23;
- *T = ((b23 * (((args.bv[N1+1] * q123 * args.av[N1] * args.pv[N1]) << args.qsv[N1+1])
- + args.bv[N1] * q23 * args.av[N1+1] * p01)) << (args.qsv[N1+2] + args.qsv[N1+3]))
- + b01 * (((args.bv[N1+3] * args.qv[N1+3] * args.av[N1+2] * p012) << args.qsv[N1+3])
+ *T = ((b23 * (((args.bv[N1+1] * q123 * args.av[N1] * args.pv[N1]) << qsv[N1+1])
+ + args.bv[N1] * q23 * args.av[N1+1] * p01)) << (qsv[N1+2] + qsv[N1+3]))
+ + b01 * (((args.bv[N1+3] * args.qv[N1+3] * args.av[N1+2] * p012) << qsv[N1+3])
+ args.bv[N1+2] * args.av[N1+3] * p0123);
break;
}
// Compute left part.
var cl_I LP, LQ, LB, LT;
var uintC LQS;
- eval_pqsab_series_aux(N1,Nm,args,&LP,&LQ,&LQS,&LB,<);
+ eval_pqsab_series_aux(N1,Nm,args,qsv,&LP,&LQ,&LQS,&LB,<);
// Compute right part.
var cl_I RP, RQ, RB, RT;
var uintC RQS;
- eval_pqsab_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RB,&RT);
+ eval_pqsab_series_aux(Nm,N2,args,qsv,(P?&RP:(cl_I*)0),&RQ,&RQS,&RB,&RT);
// Put together partial results.
if (P) { *P = LP*RP; }
*Q = LQ*RQ;
}
}
-const cl_LF eval_rational_series (uintC N, const cl_pqab_series& args, uintC len)
+template<>
+const cl_LF eval_rational_series<true> (uintC N, const cl_pqab_series& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
var cl_I Q, B, T;
- if (!args.qsv) {
- 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);
- } else {
- // Precomputation of the shift counts:
- // Split qv[n] into qv[n]*2^qsv[n].
- {
- var cl_I* qp = args.qv;
- var uintC* qsp = args.qsv;
- for (var uintC n = 0; n < N; n++, qp++, qsp++) {
- // Pull out maximal power of 2 out of *qp = args.qv[n].
- var uintC qs = 0;
- if (!zerop(*qp)) {
- qs = ord2(*qp);
- if (qs > 0)
- *qp = *qp >> qs;
- }
- *qsp = qs;
- }
- }
- // Main computation.
- var uintC QS;
- eval_pqsab_series_aux(0,N,args,NULL,&Q,&QS,&B,&T);
- return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(B*Q,len),QS);
+ // Precomputation of the shift counts:
+ // Split qv[n] into qv[n]*2^qsv[n].
+ CL_ALLOCA_STACK;
+ var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
+ var cl_I* qp = args.qv;
+ var uintC* qsp = qsv;
+ for (var uintC n = 0; n < N; n++, qp++, qsp++) {
+ *qsp = pullout_shiftcount(*qp);
}
+ // Main computation.
+ var uintC QS;
+ eval_pqsab_series_aux(0,N,args,qsv,NULL,&Q,&QS,&B,&T);
+ 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,
}
}
-const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, cl_pqab_series_stream& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
}
}
-const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len, uintC trunclen)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, cl_pqab_series_stream& args, uintC len, uintC trunclen)
{
if (N==0)
return cl_I_to_LF(0,len);
-// eval_rational_series().
+// eval_rational_series<bool>().
// General includes.
#include "cl_sysdep.h"
#include "cln/real.h"
#include "cln/exception.h"
#include "cl_LF.h"
+#include "cl_alloca.h"
namespace cln {
}
}
+template<>
+const cl_LF eval_rational_series<false> (uintC N, const cl_pqb_series& 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_pqsb_series_aux (uintC N1, uintC N2,
- const cl_pqb_series& args,
+ const cl_pqb_series& args, const uintC* qsv,
cl_I* P, cl_I* Q, uintC* QS, cl_I* B, cl_I* T)
{
switch (N2 - N1) {
case 1:
if (P) { *P = args.pv[N1]; }
*Q = args.qv[N1];
- *QS = args.qsv[N1];
+ *QS = qsv[N1];
*B = args.bv[N1];
*T = args.pv[N1];
break;
var cl_I p01 = args.pv[N1] * args.pv[N1+1];
if (P) { *P = p01; }
*Q = args.qv[N1] * args.qv[N1+1];
- *QS = args.qsv[N1] + args.qsv[N1+1];
+ *QS = qsv[N1] + qsv[N1+1];
*B = args.bv[N1] * args.bv[N1+1];
- *T = ((args.bv[N1+1] * args.qv[N1+1] * args.pv[N1]) << args.qsv[N1+1])
+ *T = ((args.bv[N1+1] * args.qv[N1+1] * args.pv[N1]) << qsv[N1+1])
+ args.bv[N1] * p01;
break;
}
if (P) { *P = p012; }
var cl_I q12 = args.qv[N1+1] * args.qv[N1+2];
*Q = args.qv[N1] * q12;
- *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2];
+ *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2];
var cl_I b12 = args.bv[N1+1] * args.bv[N1+2];
*B = args.bv[N1] * b12;
- *T = ((b12 * q12 * args.pv[N1]) << (args.qsv[N1+1] + args.qsv[N1+2]))
- + args.bv[N1] * (((args.bv[N1+2] * args.qv[N1+2] * p01) << args.qsv[N1+2])
+ *T = ((b12 * q12 * args.pv[N1]) << (qsv[N1+1] + qsv[N1+2]))
+ + args.bv[N1] * (((args.bv[N1+2] * args.qv[N1+2] * p01) << qsv[N1+2])
+ args.bv[N1+1] * p012);
break;
}
var cl_I q23 = args.qv[N1+2] * args.qv[N1+3];
var cl_I q123 = args.qv[N1+1] * q23;
*Q = args.qv[N1] * q123;
- *QS = args.qsv[N1] + args.qsv[N1+1] + args.qsv[N1+2] + args.qsv[N1+3];
+ *QS = qsv[N1] + qsv[N1+1] + qsv[N1+2] + qsv[N1+3];
var cl_I b01 = args.bv[N1] * args.bv[N1+1];
var cl_I b23 = args.bv[N1+2] * args.bv[N1+3];
*B = b01 * b23;
- *T = ((b23 * (((args.bv[N1+1] * q123 * args.pv[N1]) << args.qsv[N1+1])
- + args.bv[N1] * q23 * p01)) << (args.qsv[N1+2] + args.qsv[N1+3]))
- + b01 * (((args.bv[N1+3] * args.qv[N1+3] * p012) << args.qsv[N1+3])
+ *T = ((b23 * (((args.bv[N1+1] * q123 * args.pv[N1]) << qsv[N1+1])
+ + args.bv[N1] * q23 * p01)) << (qsv[N1+2] + qsv[N1+3]))
+ + b01 * (((args.bv[N1+3] * args.qv[N1+3] * p012) << qsv[N1+3])
+ args.bv[N1+2] * p0123);
break;
}
// Compute left part.
var cl_I LP, LQ, LB, LT;
var uintC LQS;
- eval_pqsb_series_aux(N1,Nm,args,&LP,&LQ,&LQS,&LB,<);
+ eval_pqsb_series_aux(N1,Nm,args,qsv,&LP,&LQ,&LQS,&LB,<);
// Compute right part.
var cl_I RP, RQ, RB, RT;
var uintC RQS;
- eval_pqsb_series_aux(Nm,N2,args,(P?&RP:(cl_I*)0),&RQ,&RQS,&RB,&RT);
+ eval_pqsb_series_aux(Nm,N2,args,qsv,(P?&RP:(cl_I*)0),&RQ,&RQS,&RB,&RT);
// Put together partial results.
if (P) { *P = LP*RP; }
*Q = LQ*RQ;
}
}
-const cl_LF eval_rational_series (uintC N, const cl_pqb_series& args, uintC len)
+template<>
+const cl_LF eval_rational_series<true> (uintC N, const cl_pqb_series& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
var cl_I Q, B, T;
- if (!args.qsv) {
- 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);
- } else {
- // Precomputation of the shift counts:
- // Split qv[n] into qv[n]*2^qsv[n].
- {
- var cl_I* qp = args.qv;
- var uintC* qsp = args.qsv;
- for (var uintC n = 0; n < N; n++, qp++, qsp++) {
- // Pull out maximal power of 2 out of *qp = args.qv[n].
- var uintC qs = 0;
- if (!zerop(*qp)) {
- qs = ord2(*qp);
- if (qs > 0)
- *qp = *qp >> qs;
- }
- *qsp = qs;
- }
- }
- // Main computation.
- var uintC QS;
- eval_pqsb_series_aux(0,N,args,NULL,&Q,&QS,&B,&T);
- return cl_I_to_LF(T,len) / scale_float(cl_I_to_LF(B*Q,len),QS);
+ // Precomputation of the shift counts:
+ // Split qv[n] into qv[n]*2^qsv[n].
+ CL_ALLOCA_STACK;
+ var uintC* qsv = (uintC*) cl_alloca(N*sizeof(uintC));
+ var cl_I* qp = args.qv;
+ var uintC* qsp = qsv;
+ for (var uintC n = 0; n < N; n++, qp++, qsp++) {
+ *qsp = pullout_shiftcount(*qp);
}
+ // Main computation.
+ var uintC QS;
+ eval_pqsb_series_aux(0,N,args,qsv,NULL,&Q,&QS,&B,&T);
+ 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,
}
}
-const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, cl_pqb_series_stream& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
-// eval_rational_series().
+// eval_rational_series<bool>().
// General includes.
#include "cl_sysdep.h"
}
}
-const cl_LF eval_rational_series (uintC N, const cl_q_series& args, uintC len)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, const cl_q_series& args, uintC len)
+{
+ if (N==0)
+ return cl_I_to_LF(0,len);
+ var cl_I Q, T;
+ eval_q_series_aux(0,N,args,&Q,&T);
+ return cl_I_to_LF(T,len) / cl_I_to_LF(Q,len);
+}
+
+static void eval_q_series_aux (uintC N1, uintC N2,
+ cl_q_series_stream& args,
+ cl_I* Q, cl_I* T)
+{
+ switch (N2 - N1) {
+ case 0:
+ throw runtime_exception(); break;
+ case 1: {
+ var cl_q_series_term v0 = args.next(); // [N1]
+ *Q = v0.q;
+ *T = 1;
+ break;
+ }
+ case 2: {
+ var cl_q_series_term v0 = args.next(); // [N1]
+ var cl_q_series_term v1 = args.next(); // [N1+1]
+ *Q = v0.q * v1.q;
+ *T = v1.q + 1;
+ break;
+ }
+ case 3: {
+ var cl_q_series_term v0 = args.next(); // [N1]
+ var cl_q_series_term v1 = args.next(); // [N1+1]
+ var cl_q_series_term v2 = args.next(); // [N1+2]
+ var cl_I q12 = v1.q * v2.q;
+ *Q = v0.q * q12;
+ *T = q12 + v2.q + 1;
+ break;
+ }
+ case 4: {
+ var cl_q_series_term v0 = args.next(); // [N1]
+ var cl_q_series_term v1 = args.next(); // [N1+1]
+ var cl_q_series_term v2 = args.next(); // [N1+2]
+ var cl_q_series_term v3 = args.next(); // [N1+3]
+ var cl_I q23 = v2.q * v3.q;
+ var cl_I q123 = v1.q * q23;
+ *Q = v0.q * q123;
+ *T = q123 + q23 + v3.q + 1;
+ break;
+ }
+ default: {
+ var uintC Nm = (N1+N2)/2; // midpoint
+ // Compute left part.
+ var cl_I LQ, LT;
+ eval_q_series_aux(N1,Nm,args,&LQ,<);
+ // Compute right part.
+ var cl_I RQ, RT;
+ eval_q_series_aux(Nm,N2,args,&RQ,&RT);
+ // Put together partial results.
+ *Q = LQ*RQ;
+ // S = LS + 1/LQ * RS, so T = RQ*LT + RT.
+ *T = RQ*LT + RT;
+ break;
+ }
+ }
+}
+
+template<>
+const cl_LF eval_rational_series<false> (uintC N, cl_q_series_stream& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
-// eval_rational_series().
+// eval_rational_series<bool>().
// General includes.
#include "cl_sysdep.h"
}
}
-const cl_LF eval_rational_series (uintC N, const cl_qa_series& args, uintC len)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, const cl_qa_series& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
-// eval_rational_series().
+// eval_rational_series<bool>().
// General includes.
#include "cl_sysdep.h"
}
}
-const cl_LF eval_rational_series (uintC N, const cl_qab_series& args, uintC len)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, const cl_qab_series& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
-// eval_rational_series().
+// eval_rational_series<bool>().
// General includes.
#include "cl_sysdep.h"
}
}
-const cl_LF eval_rational_series (uintC N, const cl_qb_series& args, uintC len)
+template<>
+const cl_LF eval_rational_series<false> (uintC N, const cl_qb_series& args, uintC len)
{
if (N==0)
return cl_I_to_LF(0,len);
eval_qb_series_aux(0,N,args,&Q,&B,&T);
return cl_I_to_LF(T,len) / cl_I_to_LF(B*Q,len);
}
+
+static void eval_qb_series_aux (uintC N1, uintC N2,
+ cl_qb_series_stream& args,
+ cl_I* Q, cl_I* B, cl_I* T)
+{
+ switch (N2 - N1) {
+ case 0:
+ throw runtime_exception(); break;
+ case 1: {
+ var cl_qb_series_term v0 = args.next(); // [N1]
+ *Q = v0.q;
+ *B = v0.b;
+ *T = 1;
+ break;
+ }
+ case 2: {
+ var cl_qb_series_term v0 = args.next(); // [N1]
+ var cl_qb_series_term v1 = args.next(); // [N1+1]
+ *Q = v0.q * v1.q;
+ *B = v0.b * v1.b;
+ *T = v1.b * v1.q + v0.b;
+ break;
+ }
+ case 3: {
+ var cl_qb_series_term v0 = args.next(); // [N1]
+ var cl_qb_series_term v1 = args.next(); // [N1+1]
+ var cl_qb_series_term v2 = args.next(); // [N1+2]
+ 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.b * (v2.b * v2.q + v1.b);
+ break;
+ }
+ case 4: {
+ var cl_qb_series_term v0 = args.next(); // [N1]
+ var cl_qb_series_term v1 = args.next(); // [N1+1]
+ var cl_qb_series_term v2 = args.next(); // [N1+2]
+ var cl_qb_series_term v3 = args.next(); // [N1+3]
+ 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.b * q23)
+ + b01 * (v3.b * v3.q + v2.b);
+ break;
+ }
+ default: {
+ var uintC Nm = (N1+N2)/2; // midpoint
+ // Compute left part.
+ var cl_I LQ, LB, LT;
+ eval_qb_series_aux(N1,Nm,args,&LQ,&LB,<);
+ // Compute right part.
+ var cl_I RQ, RB, RT;
+ eval_qb_series_aux(Nm,N2,args,&RQ,&RB,&RT);
+ // Put together partial results.
+ *Q = LQ*RQ;
+ *B = LB*RB;
+ // S = LS + 1/LQ * RS, so T = RB*RQ*LT + LB*RT.
+ *T = RB*RQ*LT + LB*RT;
+ break;
+ }
+ }
+}
+
+template<>
+const cl_LF eval_rational_series<false> (uintC N, cl_qb_series_stream& args, uintC len)
+{
+ if (N==0)
+ return cl_I_to_LF(0,len);
+ var cl_I Q, B, T;
+ eval_qb_series_aux(0,N,args,&Q,&B,&T);
+ return cl_I_to_LF(T,len) / cl_I_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)).
// 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.
+//
+// As a variation, it is sometimes advantageous to factor out from q[0..N-1]
+// powers of two and put them back in again at the end of the computation by
+// left-shifting. These are distinguished by a boolean template parameter.
+// (Some combinations might not be implemented yet.)
-struct cl_rational_series {
- // To be set explicitly.
- const cl_I* pv;
- cl_I* qv;
- const cl_I* av;
- const cl_I* bv;
- uintC* qsv;
-};
-extern const cl_LF eval_rational_series (uintC N, const cl_rational_series& args, uintC len);
-
-// In each the special cases below, none of (a,b,p,q) can be NULL. But qs can
-// still be given or NULL.
+// In each of the special cases below, none of (a,b,p,q) can be NULL.
struct cl_pqab_series {
const cl_I* pv;
cl_I* qv;
const cl_I* av;
const cl_I* bv;
- uintC* qsv;
};
struct cl_pqab_series_term {
cl_I p;
// 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);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pqab_series& args, uintC len);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqab_series_stream& args, uintC len);
+template<bool shift_q> 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;
cl_I* qv;
const cl_I* bv;
- uintC* qsv;
};
struct cl_pqb_series_term {
cl_I p;
// 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);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pqb_series& args, uintC len);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqb_series_stream& args, uintC len);
+template<bool shift_q> 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;
cl_I* qv;
const cl_I* av;
- uintC* qsv;
};
struct cl_pqa_series_term {
cl_I p;
// 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);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pqa_series& args, uintC len);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pqa_series_stream& args, uintC len);
+template<bool shift_q> 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;
// 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);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_pq_series& args, uintC len);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_pq_series_stream& args, uintC len);
+template<bool shift_q> 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;
const cl_I* av;
const cl_I* bv;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_pab_series& args, uintC len);
+const cl_LF eval_rational_series (uintC N, const cl_pab_series& args, uintC len);
struct cl_pb_series {
const cl_I* pv;
const cl_I* bv;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_pb_series& args, uintC len);
+const cl_LF eval_rational_series (uintC N, const cl_pb_series& args, uintC len);
struct cl_pa_series {
const cl_I* pv;
const cl_I* av;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_pa_series& args, uintC len);
+const cl_LF eval_rational_series (uintC N, const cl_pa_series& args, uintC len);
struct cl_p_series {
const cl_I* pv;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_p_series& args, uintC len);
+const cl_LF eval_rational_series (uintC N, const cl_p_series& args, uintC len);
struct cl_qab_series {
cl_I* qv;
const cl_I* av;
const cl_I* bv;
- uintC* qsv;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_qab_series& args, uintC len);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_qab_series& args, uintC len);
struct cl_qb_series {
cl_I* qv;
const cl_I* bv;
- uintC* qsv;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_qb_series& args, uintC len);
+struct cl_qb_series_term {
+ cl_I q;
+ cl_I b;
+};
+struct cl_qb_series_stream {
+ cl_qb_series_term (*nextfn)(cl_qb_series_stream&);
+ cl_qb_series_term next () { return nextfn(*this); }
+ // Constructor.
+ cl_qb_series_stream (cl_qb_series_term (*n)(cl_qb_series_stream&)) : nextfn (n) {}
+};
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_qb_series& args, uintC len);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_qb_series_stream& args, uintC len);
struct cl_qa_series {
cl_I* qv;
const cl_I* av;
- uintC* qsv;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_qa_series& args, uintC len);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_qa_series& args, uintC len);
struct cl_q_series {
cl_I* qv;
- uintC* qsv;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_q_series& args, uintC len);
+struct cl_q_series_term {
+ cl_I q;
+};
+struct cl_q_series_stream {
+ cl_q_series_term (*nextfn)(cl_q_series_stream&);
+ cl_q_series_term next () { return nextfn(*this); }
+ // Constructor.
+ cl_q_series_stream (cl_q_series_term (*n)(cl_q_series_stream&)) : nextfn (n) {}
+};
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, const cl_q_series& args, uintC len);
+template<bool shift_q> const cl_LF eval_rational_series (uintC N, cl_q_series_stream& args, uintC len);
struct cl_ab_series {
const cl_I* av;
const cl_I* bv;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_ab_series& args, uintC len);
+const cl_LF eval_rational_series (uintC N, const cl_ab_series& args, uintC len);
struct cl_b_series {
const cl_I* bv;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_b_series& args, uintC len);
+const cl_LF eval_rational_series (uintC N, const cl_b_series& args, uintC len);
struct cl_a_series {
const cl_I* av;
};
-extern const cl_LF eval_rational_series (uintC N, const cl_a_series& args, uintC len);
+const cl_LF eval_rational_series (uintC N, const cl_a_series& args, uintC len);
struct cl__series {
};
-extern const cl_LF eval_rational_series (uintC N, const cl__series& args, uintC len);
+const cl_LF eval_rational_series (uintC N, const cl__series& args, uintC len);
// [Generalization.]
// 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<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);
+void eval_pqcd_series_aux (uintC N, cl_pqcd_series_term* args, cl_pqcd_series_result<cl_I>& Z, bool rightmost = true);
+void eval_pqcd_series_aux (uintC N, cl_pqcd_series_stream& args, cl_pqcd_series_result<cl_I>& Z, bool rightmost = true);
+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);
+const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_term* args, uintC len);
+const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len);
+const cl_LF eval_pqcd_series (uintC N, cl_pqcd_series_stream& args, uintC len, uintC trunclen);
// [Special case c(n)=1.]
// Subroutine:
// 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<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);
+void eval_pqd_series_aux (uintC N, cl_pqd_series_term* args, cl_pqd_series_result<cl_I>& Z, bool rightmost = true);
+void eval_pqd_series_aux (uintC N, cl_pqd_series_stream& args, cl_pqd_series_result<cl_I>& Z, bool rightmost = true);
+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);
+const cl_LF eval_pqd_series (uintC N, cl_pqd_series_term* args, uintC len);
+const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len);
+const cl_LF eval_pqd_series (uintC N, cl_pqd_series_stream& args, uintC len, uintC trunclen);
+
+// Helper function to divide q by the largest s such that 2^s divides q, returns s.
+inline uintC
+pullout_shiftcount(cl_I& q)
+{
+ var uintC qs = 0;
+ if (!zerop(q)) {
+ qs = ord2(q);
+ if (qs > 0)
+ q = q >> qs;
+ }
+ return qs;
+}
// Helper function to convert integer of length > trunclen to long float of
// length = trunclen.
// with appropriate N, and
// a(n) = 205*n^2+250*n+77, b(n) = 1,
// p(0) = 1, p(n) = -n^5 for n>0, q(n) = 32*(2n+1)^5.
- var uintC actuallen = len+2; // 2 Schutz-Digits
+ var uintC actuallen = len+2; // 2 guard digits
var uintC N = ceiling(actuallen*intDsize,10);
// 1024^-N <= 2^(-intDsize*actuallen).
- var cl_LF sum = eval_rational_series(N,series,actuallen,actuallen);
+ var cl_LF sum = eval_rational_series<false>(N,series,actuallen,actuallen);
return scale_float(shorten(sum,len),-1);
}
// Bit complexity (N := len): O(log(N)^2*M(N)).
// zeta(s) = 1/(1-2^(1-s)) sum(n=0..infty, (-1)^n/(n+1)^s),
// with convergence acceleration through exp(x), and evaluated
// using the binary-splitting algorithm.
- var uintC actuallen = len+2; // 2 Schutz-Digits
+ var uintC actuallen = len+2; // 2 guard digits
var uintC x = (uintC)(0.693148*intDsize*actuallen)+1;
var uintC N = (uintC)(2.718281828*x);
CL_ALLOCA_STACK;
// Method:
// zeta(s) = 1/(1-2^(1-s)) sum(n=0..infty, (-1)^n/(n+1)^s),
// with Cohen-Villegas-Zagier convergence acceleration.
- var uintC actuallen = len+2; // 2 Schutz-Digits
+ var uintC actuallen = len+2; // 2 guard digits
var uintC N = (uintC)(0.39321985*intDsize*actuallen)+1;
var cl_I fterm = 2*(cl_I)N*(cl_I)N;
var cl_I fsum = fterm;
projects / cln.git / commitdiff