]> www.ginac.de Git - cln.git/blobdiff - src/float/transcendental/cl_F_atanhx.cc
Fix linking problems on some platforms caused by inline/non-inline versions
[cln.git] / src / float / transcendental / cl_F_atanhx.cc
index fc732eb1badbbff78adad4d33f06cfd83bd55ca2..c2f0475ebd8e4b77d4e0ec687aa00036f41c7d72 100644 (file)
@@ -15,8 +15,7 @@
 #include "cln/float.h"
 #include "cl_low.h"
 
-#undef MAYBE_INLINE
-#define MAYBE_INLINE inline
+#include "cl_inline.h"
 #include "cl_LF_zerop.cc"
 #include "cl_LF_minusp.cc"
 #include "cl_LF_exponent.cc"
@@ -32,8 +31,8 @@ namespace cln {
 //   (denn bei e<=-d/2 ist x^2 < 2^(-d), also
 //   1 <= atanh(x)/x = 1+x^2/3+x^4/5+... < 1+x^2/2 < 1+2^(-d-1) < 1+2^(-d),
 //   also ist atanh(x)/x, auf d Bits gerundet, gleich 1.0).
-// Bei großem d verwende die Formel ln((1+x)/(1-x))/2 (asymptotisch schneller),
-//   aber erhöhe die Genauigkeit, so daß beim Bilden von 1+x keine Bits verloren
+// Bei großem d verwende die Formel ln((1+x)/(1-x))/2 (asymptotisch schneller),
+//   aber erhöhe die Genauigkeit, so daß beim Bilden von 1+x keine Bits verloren
 //   gehen.
 // Bei e<=-sqrt(d) verwende die Potenzreihe
 //   atanh(x)/x = sum(j=0..inf,(x^2)^j/(2j+1)):
@@ -49,41 +48,41 @@ namespace cln {
 
 const cl_LF atanhx (const cl_LF& x)
 {
-       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 <= (sintL)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2)
+       var uintC actuallen = TheLfloat(x)->len;
+       var uintC d = float_digits(x);
+       var sintE e = float_exponent_inline(x);
+       if (e <= (sintC)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2)
                return x; // ja -> x als Ergebnis
        if (actuallen >= 34) {
                DeclareType(cl_LF,x);
-               var cl_LF xx = extend(x,TheLfloat(x)->len+ceiling((uintL)(-e),intDsize));
+               var cl_LF xx = extend(x,TheLfloat(x)->len+ceiling((uintE)(-e),intDsize));
                return cl_float(scale_float(ln((1+xx)/(1-xx)),-1),x);
        }
-       var uintL k = 0; // Rekursionszähler k:=0
+       var uintL k = 0; // Rekursionszähler k:=0
        // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
        // angewandt werden. limit_slope = 1.0 ist schlecht (ca. 15% zu
        // schlecht). Ein guter Wert ist:
-       // für naive1: limit_scope = 0.625 = 5/8,
-       // für naive2: limit_scope = 0.4 = 13/32.
-       var uintL sqrt_d = floor(isqrt(d)*13,32); // limit_slope*floor(sqrt(d))
+       // für naive1: limit_scope = 0.625 = 5/8,
+       // für naive2: limit_scope = 0.4 = 13/32.
+       var uintL sqrt_d = floor(isqrtC(d)*13,32); // limit_slope*floor(sqrt(d))
        var cl_LF xx = x;
        if (e >= (sintL)(-sqrt_d)) {
-               // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
+               // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
                var sintL e_limit = 1+sqrt_d; // 1+limit_slope*floor(sqrt(d))
                xx = recip(abs(xx)); // 1/|x|
                do {
-                 // nächstes x nach der Formel x := x+sqrt(x^2 - 1) berechnen:
+                 // nächstes x nach der Formel x := x+sqrt(x^2 - 1) berechnen:
                  xx = sqrt(square(xx) + cl_float(-1,xx)) + xx;
                  k = k+1;
-               } until (float_exponent(xx) > e_limit);
+               } until (float_exponent_inline(xx) > e_limit);
                // Schleifenende mit Exponent(x) > 1+limit_slope*floor(sqrt(d)),
                // also x >= 2^(1+limit_slope*floor(sqrt(d))),
                // also 1/x <= 2^(-1-limit_slope*floor(sqrt(d))).
                // Nun kann die Potenzreihe auf 1/x angewandt werden.
                xx = recip(xx);
-               if (minusp(x))
+               if (minusp_inline(x))
                        xx = - xx; // Vorzeichen wieder rein
        }
        // Potenzreihe anwenden:
@@ -105,7 +104,7 @@ const cl_LF atanhx (const cl_LF& x)
        } else {
                // naive2:
                // floating-point representation with smooth precision reduction
-               var cl_LF eps = scale_float(b,-(sintL)d-10);
+               var cl_LF eps = scale_float(b,-(sintC)d-10);
                loop {
                        var cl_LF new_sum = sum + LF_to_LF(b/(cl_I)i,actuallen); // (+ sum (/ b i))
                        if (new_sum == sum) // = sum ?
@@ -129,22 +128,22 @@ const cl_F atanhx (const cl_F& x)
        }
        if (zerop(x))
                return x;
-       var uintL d = float_digits(x);
-       var sintL e = float_exponent(x);
-       if (e <= (sintL)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2)
+       var uintC d = float_digits(x);
+       var sintE e = float_exponent(x);
+       if (e <= (sintC)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2)
                return x; // ja -> x als Ergebnis
-       var uintL k = 0; // Rekursionszähler k:=0
+       var uintL k = 0; // Rekursionszähler k:=0
        // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
        // angewandt werden. limit_slope = 1.0 ist schlecht (ca. 15% zu
        // schlecht). Ein guter Wert ist limit_scope = 0.625 = 5/8.
-       var uintL sqrt_d = floor(isqrt(d)*5,8); // limit_slope*floor(sqrt(d))
+       var uintL sqrt_d = floor(isqrtC(d)*5,8); // limit_slope*floor(sqrt(d))
        var cl_F xx = x;
        if (e >= (sintL)(-sqrt_d)) {
-               // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
+               // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
                var sintL e_limit = 1+sqrt_d; // 1+limit_slope*floor(sqrt(d))
                xx = recip(abs(xx)); // 1/|x|
                do {
-                 // nächstes x nach der Formel x := x+sqrt(x^2 - 1) berechnen:
+                 // nächstes x nach der Formel x := x+sqrt(x^2 - 1) berechnen:
                  xx = sqrt(square(xx) + cl_float(-1,xx)) + xx;
                  k = k+1;
                } until (float_exponent(xx) > e_limit);