Get rid of CL_REQUIRE/CL_PROVIDE(cl_F_ln10_var).

[cln.git] / src / float / transcendental / cl_LF_coshsinh_aux.cc

diff --git a/src/float/transcendental/cl_LF_coshsinh_aux.cc b/src/float/transcendental/cl_LF_coshsinh_aux.cc
index 33d794d1b86e9c35ec2d1e085337739a505790f8..fab71c5c24127025985dbe2667a0536845d688d4 100644 (file)
--- a/src/float/transcendental/cl_LF_coshsinh_aux.cc
+++ b/src/float/transcendental/cl_LF_coshsinh_aux.cc
@@ -9,30 +9,32 @@
 
 // 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"
+#include "cln/exception.h"
 
 #undef floor
-#include <math.h>
+#include <cmath>
 #define floor cln_floor
 
+namespace cln {
+
 // Computing cosh(x) = sqrt(1+sinh(x)^2) instead of computing separately
 // by a power series evaluation brings 20% speedup, even more for small lengths.
 #define TRIVIAL_SPEEDUP
 
-const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintL lq, uintC len)
+const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintE lq, uintC len)
 {
  {     Mutable(cl_I,p);
-       var uintL lp = integer_length(p); // now |p| < 2^lp.
-       if (!(lp <= lq)) cl_abort();
+       var uintE lp = integer_length(p); // now |p| < 2^lp.
+       if (!(lp <= lq)) throw runtime_exception();
        lp = lq - lp; // now |p/2^lq| < 2^-lp.
        // Minimize lq (saves computation time).
        {
-               var uintL lp2 = ord2(p);
+               var uintC lp2 = ord2(p);
                if (lp2 > 0) {
                        p = p >> lp2;
                        lq = lq - lp2;
@@ -50,8 +52,8 @@ const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintL lq, uintC 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)
@@ -63,16 +65,15 @@ const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintL lq, uintC len)
        // Third approximation:
        //   N2 = ceiling(M*log(2)/(log(N1)-1+lp*log(2))), slightly too large.
        //   N = N2+2, two more terms for safety.
-       var uintL N0 = intDsize*actuallen;
-       var uintL N1 = (uintL)(0.693147*intDsize*actuallen/(log((double)N0)-1.0+0.693148*lp));
-       var uintL N2 = (uintL)(0.693148*intDsize*actuallen/(log((double)N1)-1.0+0.693147*lp))+1;
-       var uintL N = N2+2;
+       var uintC N0 = intDsize*actuallen;
+       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;
        N = ceiling(N,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 uintL* qsv = (uintL*) cl_alloca(N*sizeof(uintL));
-       var uintL n;
+       var uintC n;
        var cl_I p2 = square(p);
        var cl_LF sinhsum;
        {
@@ -83,8 +84,8 @@ const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintL lq, uintC len)
                        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();
@@ -100,8 +101,8 @@ const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintL lq, uintC len)
                        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();
@@ -110,8 +111,9 @@ const cl_LF_cosh_sinh_t cl_coshsinh_aux (const cl_I& p, uintL lq, uintC len)
        #else // TRIVIAL_SPEEDUP
        var cl_LF coshsum = sqrt(cl_I_to_LF(1,actuallen) + square(sinhsum));
        #endif
-       return cl_LF_cosh_sinh_t(shorten(coshsum,len),shorten(sinhsum,len)); // verkürzen und fertig
+       return cl_LF_cosh_sinh_t(shorten(coshsum,len),shorten(sinhsum,len)); // verkürzen und fertig
 }}
 // Bit complexity (N = len, and if p has length O(log N) and ql = O(log N)):
 // O(log(N)*M(N)).
 
+}  // namespace cln