X-Git-Url: https://ginac.de/CLN/cln.git//cln.git?a=blobdiff_plain;ds=inline;f=src%2Ffloat%2Ftranscendental%2Fcl_LF_eulerconst.cc;h=0cc8bbfc00670ece9e8058c723751bf7ae00499a;hb=eedac861e082bb3e12392e0037842470310f80d6;hp=d6cdcf71e8256bc5e7d90318de70cfc3a589ad4e;hpb=dd9e0f894eec7e2a8cf85078330ddc0a6639090b;p=cln.git diff --git a/src/float/transcendental/cl_LF_eulerconst.cc b/src/float/transcendental/cl_LF_eulerconst.cc index d6cdcf7..0cc8bbf 100644 --- a/src/float/transcendental/cl_LF_eulerconst.cc +++ b/src/float/transcendental/cl_LF_eulerconst.cc @@ -1,4 +1,4 @@ -// cl_eulerconst(). +// eulerconst(). // General includes. #include "cl_sysdep.h" @@ -9,12 +9,14 @@ // 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 @@ -69,14 +71,14 @@ const cl_LF compute_eulerconst_expintegral (uintC len) // 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); @@ -84,15 +86,15 @@ const cl_LF compute_eulerconst_expintegral (uintC len) } 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(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)). @@ -144,15 +146,15 @@ const cl_LF compute_eulerconst_expintegral1 (uintC len) // 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!. @@ -172,7 +174,7 @@ const cl_LF compute_eulerconst_expintegral1 (uintC len) } } 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). @@ -184,12 +186,12 @@ const cl_LF compute_eulerconst_expintegral1 (uintC len) // 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); @@ -209,7 +211,7 @@ const cl_LF compute_eulerconst_expintegral2 (uintC len) 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)). @@ -274,16 +276,16 @@ const cl_LF compute_eulerconst_expintegral2 (uintC len) 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. @@ -303,7 +305,7 @@ const cl_LF compute_eulerconst_besselintegral1 (uintC len) } } 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). @@ -324,15 +326,15 @@ const cl_LF compute_eulerconst_besselintegral2 (uintC len) // 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); @@ -341,8 +343,8 @@ const cl_LF compute_eulerconst_besselintegral2 (uintC len) 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(N,fseries,actuallen); for (n = 0; n < N; n++) { pv[n].~cl_I(); qv[n].~cl_I(); @@ -363,15 +365,15 @@ const cl_LF compute_eulerconst_besselintegral2 (uintC len) } 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(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). @@ -406,14 +408,14 @@ struct cl_rational_series_for_g : cl_pqa_series_stream { }; 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); @@ -422,17 +424,17 @@ const cl_LF compute_eulerconst_besselintegral3 (uintC len) 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(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(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). @@ -442,33 +444,38 @@ const cl_LF compute_eulerconst_besselintegral3 (uintC len) // the sums. const cl_LF compute_eulerconst_besselintegral4 (uintC len) { - var uintC actuallen = len+2; // 2 Schutz-Digits - var 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 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)). @@ -498,22 +505,24 @@ const cl_LF compute_eulerconst (uintC len) 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