]> www.ginac.de Git - cln.git/blobdiff - src/float/transcendental/cl_F_sinhx.cc
Fix linking problems on some platforms caused by inline/non-inline versions
[cln.git] / src / float / transcendental / cl_F_sinhx.cc
index fc658d09a1ad9483c298787c6a0cef7ea6e9bf6d..820049f0dcbb5dc3856e4d12117669067ef5d0a6 100644 (file)
@@ -16,8 +16,7 @@
 #include "cl_LF.h"
 #include "cln/integer.h"
 
-#undef MAYBE_INLINE
-#define MAYBE_INLINE inline
+#include "cl_inline.h"
 #include "cl_LF_zerop.cc"
 #include "cl_LF_exponent.cc"
 
@@ -42,7 +41,7 @@ const cl_F sinhxbyx_naive (const cl_F& x)
 //   berechne rekursiv z:=(sinh(y)/y)^2 und liefere z*(1+y^2*z).
 // [Die Grenze sqrt(d) ergibt sich so:
 //  Man braucht bei der Potenzreihe mit x=2^-k etwa j Glieder, mit
-//  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand beträgt etwa 2.8*(j/2)
+//  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand beträgt etwa 2.8*(j/2)
 //  Multiplikationen von d-Bit-Zahlen. Bei Halbierungen bis x=2^-k ist der
 //  Gesamtaufwand etwa 2*(k+e)+1.4*j(k). Dieses minimieren nach k: Soll sein
 //  -1.4 = d/dk j(k) = (d/dj k(j))^-1 = - j^2/(d+j)*ln 2, also j^2=2(d+j),
@@ -51,16 +50,16 @@ const cl_F sinhxbyx_naive (const cl_F& x)
 
        if (zerop(x))
                return cl_float(1,x);
-       var uintL d = float_digits(x);
-       var sintL e = float_exponent(x);
-       if (e <= (1-(sintL)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
+       var uintC d = float_digits(x);
+       var sintE e = float_exponent(x);
+       if (e <= (1-(sintC)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
                return cl_float(1,x); // ja -> 1.0 als Ergebnis
  {     Mutable(cl_F,x);
        // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
-       // angewandt werden. Wähle limit_slope = 13/32 = 0.4.
-       var sintL e_limit = -1-floor(isqrt(d)*13,32); // -1-floor(sqrt(d))
+       // angewandt werden. Wähle limit_slope = 13/32 = 0.4.
+       var sintL e_limit = -1-floor(isqrtC(d)*13,32); // -1-floor(sqrt(d))
        if (e > e_limit) {
-               // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
+               // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
                x = scale_float(x,e_limit-e);
                // Neuer Exponent = e_limit.
        }
@@ -82,7 +81,7 @@ const cl_F sinhxbyx_naive (const cl_F& x)
        while (e > e_limit) {
                z = z + x2 * square(z);
                x2 = scale_float(x2,2); // x^2 := x^2*4
-               e_limit++;
+               e--;
        }
        return z;
 }}
@@ -105,28 +104,28 @@ const cl_LF sinhx_naive (const cl_LF& x)
 //   berechne rekursiv z:=sinh(y)^2 und liefere 4*z*(1+z) = (1+2*z)^2-1.
 // [Die Grenze sqrt(d) ergibt sich so:
 //  Man braucht bei der Potenzreihe mit x=2^-k etwa j Glieder, mit
-//  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand beträgt etwa 2.8*(j/2)
+//  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand beträgt etwa 2.8*(j/2)
 //  Multiplikationen von d-Bit-Zahlen. Bei Halbierungen bis x=2^-k ist der
 //  Gesamtaufwand etwa 2*(k+e)+1.4*j(k). Dieses minimieren nach k: Soll sein
 //  -1.4 = d/dk j(k) = (d/dj k(j))^-1 = - j^2/(d+j)*ln 2, also j^2=2(d+j),
 //  grob j=sqrt(2d) und damit k=sqrt(d).]
 // Aufwand: asymptotisch d^2.5 .
 
-       if (zerop(x))
+       if (zerop_inline(x))
                return x;
-       var uintL actuallen = TheLfloat(x)->len;
-       var uintL d = float_digits(x);
-       var sintL e = float_exponent(x);
-       if (e <= (1-(sintL)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
-               return x; // ja -> x als Ergebnis
+       var uintC actuallen = TheLfloat(x)->len;
+       var uintC d = float_digits(x);
+       var sintE e = float_exponent_inline(x);
+       if (e <= (1-(sintC)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
+               return square(x); // ja -> x^2 als Ergebnis
  {     Mutable(cl_LF,x);
-       var sintL ee = e;
+       var sintE ee = e;
        // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
-       // angewandt werden. Ein guter Wert für naive1 ist limit_slope = 0.6,
-       // für naive3 aber limit_slope = 0.5.
-       var sintL e_limit = -1-floor(isqrt(d),2); // -1-floor(sqrt(d))
+       // angewandt werden. Ein guter Wert für naive1 ist limit_slope = 0.6,
+       // für naive3 aber limit_slope = 0.5.
+       var sintL e_limit = -1-floor(isqrtC(d),2); // -1-floor(sqrt(d))
        if (e > e_limit) {
-               // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
+               // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
                x = scale_float(x,e_limit-e);
                ee = e_limit; // Neuer Exponent = e_limit.
        }
@@ -147,7 +146,7 @@ const cl_LF sinhx_naive (const cl_LF& x)
                        b = round1(round1(The(cl_LF)(b*a)),(cl_I)((i+1)*(i+2)));
                        i = i+2;
                }
-               powser_value = scale_float(cl_float(sum,x),-d);
+               powser_value = scale_float(cl_float(sum,x),-(sintC)d);
        } else if (actuallen <= 7) { // Break-even-Point before extendsqrt: N<=6
                // naive2:
                // floating-point representation
@@ -166,7 +165,7 @@ const cl_LF sinhx_naive (const cl_LF& x)
                // naive3:
                // floating-point representation with smooth precision reduction
                var cl_LF b = x; // b := x
-               var cl_LF eps = scale_float(b,-(sintL)d-10);
+               var cl_LF eps = scale_float(b,-(sintC)d-10);
                var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
                loop {
                        var cl_LF new_sum = sum + LF_to_LF(b,actuallen);
@@ -182,7 +181,7 @@ const cl_LF sinhx_naive (const cl_LF& x)
        var cl_LF z = square(powser_value); // sinh^2 als Ergebnis
        while (e > e_limit) {
                z = square(cl_float(1,x) + scale_float(z,1)) - cl_float(1,x); // z := (1+2*z)^2-1
-               e_limit++;
+               e--;
        }
        return z;
 }}