]> www.ginac.de Git - cln.git/blobdiff - src/float/transcendental/cl_LF_exp1.cc
Update known-to-work-with compilers.
[cln.git] / src / float / transcendental / cl_LF_exp1.cc
index 0f55fe5596dd5a9ae66624c3b7bb16dbfd1310d6..c1f1acee5a98de5203bf0114a56e5f1fedbc079c 100644 (file)
@@ -1,19 +1,18 @@
 // exp1().
 
 // General includes.
-#include "cl_sysdep.h"
+#include "base/cl_sysdep.h"
 
 // Specification.
-#include "cl_F_tran.h"
+#include "float/transcendental/cl_F_tran.h"
 
 
 // Implementation.
 
 #include "cln/lfloat.h"
-#include "cl_LF_tran.h"
-#include "cl_LF.h"
+#include "float/transcendental/cl_LF_tran.h"
+#include "float/lfloat/cl_LF.h"
 #include "cln/integer.h"
-#include "cl_alloca.h"
 
 #undef floor
 #include <cmath>
@@ -26,8 +25,8 @@ const cl_LF compute_exp1 (uintC len)
        // 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)
@@ -42,20 +41,22 @@ const cl_LF compute_exp1 (uintC len)
        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)).
@@ -77,13 +78,13 @@ const cl_LF compute_exp1 (uintC len)
 
 const cl_LF exp1 (uintC len)
 {
-       var uintC oldlen = TheLfloat(cl_LF_exp1)->len; // vorhandene Länge
+       var uintC oldlen = TheLfloat(cl_LF_exp1())->len; // vorhandene Länge
        if (len < oldlen)
-               return shorten(cl_LF_exp1,len);
+               return shorten(cl_LF_exp1(),len);
        if (len == oldlen)
-               return cl_LF_exp1;
+               return cl_LF_exp1();
 
-       // TheLfloat(cl_LF_exp1)->len um mindestens einen konstanten Faktor
+       // TheLfloat(cl_LF_exp1())->len um mindestens einen konstanten Faktor
        // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
        var uintC newlen = len;
        oldlen += floor(oldlen,2); // oldlen * 3/2
@@ -91,8 +92,8 @@ const cl_LF exp1 (uintC len)
                newlen = oldlen;
 
        // gewünschte > vorhandene Länge -> muß nachberechnen:
-       cl_LF_exp1 = compute_exp1(newlen); // (exp 1)
-       return (len < newlen ? shorten(cl_LF_exp1,len) : cl_LF_exp1);
+       cl_LF_exp1() = compute_exp1(newlen); // (exp 1)
+       return (len < newlen ? shorten(cl_LF_exp1(),len) : cl_LF_exp1());
 }
 
 }  // namespace cln