X-Git-Url: https://ginac.de/CLN/cln.git//cln.git?a=blobdiff_plain;f=src%2Ffloat%2Ftranscendental%2Fcl_LF_eulerconst.cc;h=3e13c1917c242b59a5d7cd34aa9ae6e00b7fbbc3;hb=4109ef07d74b779bb4d4371c24a6364302f464bb;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..3e13c19 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,13 @@
 
 // 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 "cl_alloca.h"
-#include "cl_abort.h"
+
+namespace cln {
 
 #if 0 // works, but besselintegral4 is always faster
 
@@ -70,13 +71,13 @@ const cl_LF compute_eulerconst_expintegral (uintC len)
 	// 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 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);
@@ -145,14 +146,14 @@ const cl_LF compute_eulerconst_expintegral1 (uintC len)
 	// 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 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!.
@@ -185,11 +186,11 @@ const cl_LF compute_eulerconst_expintegral1 (uintC len)
 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 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);
@@ -275,15 +276,15 @@ 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 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.
@@ -325,14 +326,14 @@ const cl_LF compute_eulerconst_besselintegral2 (uintC len)
 	// (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 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);
@@ -407,13 +408,13 @@ 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 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);
@@ -443,19 +444,29 @@ const cl_LF compute_eulerconst_besselintegral3 (uintC len)
 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 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_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);
-	}
+	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 (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,args,sums);
+	eval_pqd_series_aux(N,series,sums);
 	// 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)).
@@ -463,11 +474,6 @@ const cl_LF compute_eulerconst_besselintegral4 (uintC len)
 	  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
 }
 // Bit complexity (N = len): O(log(N)^2*M(N)).
@@ -498,7 +504,7 @@ 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
 	if (len < oldlen)
@@ -517,3 +523,5 @@ const cl_LF cl_eulerconst (uintC len)
 	cl_LF_eulerconst = compute_eulerconst(newlen);
 	return (len < newlen ? shorten(cl_LF_eulerconst,len) : cl_LF_eulerconst);
 }
+
+}  // namespace cln