// sinhxbyx(), sinhx().
// 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/float.h"
-#include "cl_low.h"
-#include "cl_F.h"
+#include "base/cl_low.h"
+#include "float/cl_F.h"
#include "cln/lfloat.h"
-#include "cl_LF.h"
+#include "float/lfloat/cl_LF.h"
#include "cln/integer.h"
-#undef MAYBE_INLINE
-#define MAYBE_INLINE inline
-#include "cl_LF_zerop.cc"
-#include "cl_LF_exponent.cc"
+#include "base/cl_inline.h"
+#include "float/lfloat/elem/cl_LF_zerop.cc"
+#include "float/lfloat/misc/cl_LF_exponent.cc"
namespace cln {
// 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),
if (zerop(x))
return cl_float(1,x);
var uintC d = float_digits(x);
- var sintL e = float_exponent(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.
}
while (e > e_limit) {
z = z + x2 * square(z);
x2 = scale_float(x2,2); // x^2 := x^2*4
- e_limit++;
+ e--;
}
return z;
}}
// 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 uintC actuallen = TheLfloat(x)->len;
var uintC d = float_digits(x);
- var sintL e = float_exponent(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.
}
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;
}}