-// cl_eulerconst().
+// eulerconst().
// General includes.
#include "cl_sysdep.h"
// Implementation.
-#include "cl_lfloat.h"
+#include "cln/lfloat.h"
#include "cl_LF_tran.h"
#include "cl_LF.h"
-#include "cl_integer.h"
+#include "cln/integer.h"
+#include "cln/real.h"
#include "cl_alloca.h"
-#include "cl_abort.h"
+
+namespace cln {
#if 0 // works, but besselintegral4 is always faster
// 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 uintL z = (uintL)(0.693148*intDsize*actuallen)+1;
- var uintL N = (uintL)(3.591121477*z);
+ 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_I* bv = (cl_I*) cl_alloca(N*sizeof(cl_I));
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 uintL n;
+ var uintC n;
for (n = 0; n < N; n++) {
init1(cl_I, bv[n]) (n+1);
init1(cl_I, pv[n]) (n==0 ? (cl_I)z : -(cl_I)z);
}
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();
qv[n].~cl_I();
}
fsum = fsum - ln(cl_I_to_LF(z,actuallen)); // log(z) subtrahieren
- return shorten(fsum,len); // verkürzen und fertig
+ return shorten(fsum,len); // verkürzen und fertig
}
// Bit complexity (N = len): O(log(N)^2*M(N)).
// 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 uintL x = (uintL)(0.693148*intDsize*actuallen)+1;
- var uintL N = (uintL)(2.718281828*x);
+ 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);
var cl_LF fterm = one;
var cl_LF fsum = fterm;
var cl_LF gterm = cl_I_to_LF(0,actuallen);
var cl_LF gsum = gterm;
- var uintL n;
+ var uintC n;
// After n loops
// fterm = x^n/n!, fsum = 1 + x/1! + ... + x^n/n!,
// gterm = H_n*x^n/n!, gsum = H_1*x/1! + ... + H_n*x^n/n!.
}
}
var cl_LF result = gsum/fsum - ln(cl_I_to_LF(x,actuallen));
- return shorten(result,len); // verkürzen und fertig
+ return shorten(result,len); // verkürzen und fertig
}
// Bit complexity (N = len): O(N^2).
// the sums.
const cl_LF compute_eulerconst_expintegral2 (uintC len)
{
- var uintC actuallen = len+2; // 2 Schutz-Digits
- var uintL x = (uintL)(0.693148*intDsize*actuallen)+1;
- var uintL N = (uintL)(2.718281828*x);
+ 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;
var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
- var uintL n;
+ var uintC n;
for (n = 0; n < N; n++) {
init1(cl_I, args[n].p) (x);
init1(cl_I, args[n].q) (n+1);
args[n].q.~cl_I();
args[n].d.~cl_I();
}
- return shorten(result,len); // verkürzen und fertig
+ return shorten(result,len); // verkürzen und fertig
}
// Bit complexity (N = len): O(log(N)^2*M(N)).
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 uintL sx = (uintL)(0.25*0.693148*intDsize*actuallen)+1;
- var uintL N = (uintL)(3.591121477*sx);
+ 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);
- var cl_LF eps = scale_float(cl_I_to_LF(1,LF_minlen),-(sintL)(sx*2.88539+10));
+ var cl_LF eps = scale_float(cl_I_to_LF(1,LF_minlen),-(sintC)(sx*2.88539+10));
var cl_LF fterm = cl_I_to_LF(1,actuallen);
var cl_LF fsum = fterm;
var cl_LF gterm = cl_I_to_LF(0,actuallen);
var cl_LF gsum = gterm;
- var uintL n;
+ var uintC n;
// After n loops
// fterm = x^n/n!^2, fsum = 1 + x/1!^2 + ... + x^n/n!^2,
// gterm = H_n*x^n/n!^2, gsum = H_1*x/1!^2 + ... + H_n*x^n/n!^2.
}
}
var cl_LF result = gsum/fsum - ln(cl_I_to_LF(sx,actuallen));
- return shorten(result,len); // verkürzen und fertig
+ return shorten(result,len); // verkürzen und fertig
}
// Bit complexity (N = len): O(N^2).
// 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 uintL sx = (uintL)(0.25*0.693148*intDsize*actuallen)+1;
- var uintL N = (uintL)(3.591121477*sx);
+ 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);
CL_ALLOCA_STACK;
var cl_I* av = (cl_I*) cl_alloca(N*sizeof(cl_I));
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 uintL n;
+ var uintC n;
// Evaluate f(x).
init1(cl_I, pv[0]) (1);
init1(cl_I, qv[0]) (1);
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();
qv[n].~cl_I();
}
var cl_LF result = gsum/fsum - ln(cl_I_to_LF(sx,actuallen));
- return shorten(result,len); // verkürzen und fertig
+ return shorten(result,len); // verkürzen und fertig
}
// Bit complexity (N = len): O(N^2).
// Memory consumption: O(N^2).
};
const cl_LF compute_eulerconst_besselintegral3 (uintC len)
{
- var uintC actuallen = len+1; // 1 Schutz-Digit
- var uintL sx = (uintL)(0.25*0.693148*intDsize*actuallen)+1;
- var uintL N = (uintL)(3.591121477*sx);
+ 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);
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 uintL n;
+ var uintC n;
// Evaluate f(x).
init1(cl_I, pv[0]) (1);
init1(cl_I, qv[0]) (1);
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
+ return shorten(result,len); // verkürzen und fertig
}
// Bit complexity (N = len): O(N^2).
// the sums.
const cl_LF compute_eulerconst_besselintegral4 (uintC len)
{
- var uintC actuallen = len+2; // 2 Schutz-Digits
- var uintL sx = (uintL)(0.25*0.693148*intDsize*actuallen)+1;
- var uintL N = (uintL)(3.591121477*sx);
- var cl_I x = square((cl_I)sx);
- CL_ALLOCA_STACK;
- var cl_pqd_series_term* args = (cl_pqd_series_term*) cl_alloca(N*sizeof(cl_pqd_series_term));
- var uintL n;
- for (n = 0; n < N; n++) {
- init1(cl_I, args[n].p) (x);
- init1(cl_I, args[n].q) (square((cl_I)(n+1)));
- init1(cl_I, args[n].d) (n+1);
- }
- var cl_pqd_series_result sums;
- eval_pqd_series_aux(N,args,sums);
+ 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));
+ struct rational_series_stream : cl_pqd_series_stream {
+ uintC n;
+ cl_I x;
+ static cl_pqd_series_term computenext (cl_pqd_series_stream& thisss)
+ {
+ var rational_series_stream& thiss = (rational_series_stream&)thisss;
+ var uintC n = thiss.n;
+ var cl_pqd_series_term result;
+ result.p = thiss.x;
+ result.q = square(cl_I(n+1));
+ result.d = n+1;
+ thiss.n = n+1;
+ return result;
+ }
+ rational_series_stream (uintC n_, const cl_I& x_)
+ : cl_pqd_series_stream (rational_series_stream::computenext),
+ 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));
- for (n = 0; n < N; n++) {
- args[n].p.~cl_I();
- args[n].q.~cl_I();
- args[n].d.~cl_I();
- }
- return shorten(result,len); // verkürzen und fertig
+ 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)).
return compute_eulerconst_besselintegral1(len);
}
-const cl_LF cl_eulerconst (uintC len)
+const cl_LF eulerconst (uintC len)
{
- var uintC oldlen = TheLfloat(cl_LF_eulerconst)->len; // vorhandene Länge
+ var uintC oldlen = TheLfloat(cl_LF_eulerconst)->len; // vorhandene Länge
if (len < oldlen)
return shorten(cl_LF_eulerconst,len);
if (len == oldlen)
return cl_LF_eulerconst;
// TheLfloat(cl_LF_eulerconst)->len um mindestens einen konstanten Faktor
- // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
+ // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
var uintC newlen = len;
oldlen += floor(oldlen,2); // oldlen * 3/2
if (newlen < oldlen)
newlen = oldlen;
- // gewünschte > vorhandene Länge -> muß nachberechnen:
+ // gewünschte > vorhandene Länge -> muß nachberechnen:
cl_LF_eulerconst = compute_eulerconst(newlen);
return (len < newlen ? shorten(cl_LF_eulerconst,len) : cl_LF_eulerconst);
}
+
+} // namespace cln
git://www.ginac.de / cln.git / blobdiff