src/float/transcendental/cl_F_sinhx.cc

   1 // sinhxbyx(), sinhx().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_F_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/float.h"
  13 #include "cl_low.h"
  14 #include "cl_F.h"
  15 #include "cln/lfloat.h"
  16 #include "cl_LF.h"
  17 #include "cln/integer.h"
  18
  19 #undef MAYBE_INLINE
  20 #define MAYBE_INLINE inline
  21 #include "cl_LF_zerop.cc"
  22 #include "cl_LF_exponent.cc"
  23
  24 namespace cln {
  25
  26 // sinhxbyx is mainly for cl_SF, cl_FF, cl_DF, where we want to avoid underflow.
  27
  28 const cl_F sinhxbyx_naive (const cl_F& x)
  29 {
  30 // Methode:
  31 // e := Exponent aus (decode-float x), d := (float-digits x)
  32 // Bei x=0.0 oder e<=(1-d)/2 liefere 1.0
  33 //   (denn bei e<=(1-d)/2 ist x^2/6 < x^2/4 < 2^(1-d)/4 = 2^(-d-1), also
  34 //   1 <= sinh(x)/x = 1+x^2/6+... < 1+2^(-d-1), also 1 <= (sinh(x)/x)^2 < 1+2^(-d),
  35 //   also ist (sinh(x)/x)^2, auf d Bits gerundet, gleich 1.0).
  36 // Bei e<=-sqrt(d) verwende die Potenzreihe
  37 //   sinh(x)/x = sum(j=0..inf,(x^2)^j/(2j+1)!):
  38 //   a:=x^2, b:=1, i:=1, sum:=0,
  39 //   while (/= sum (setq sum (+ sum b))) do b:=b*a/((i+1)*(i+2)), i:=i+2.
  40 //   Ergebnis sum^2.
  41 // Sonst setze y := x/2 = (scale-float x -1),
  42 //   berechne rekursiv z:=(sinh(y)/y)^2 und liefere z*(1+y^2*z).
  43 // [Die Grenze sqrt(d) ergibt sich so:
  44 //  Man braucht bei der Potenzreihe mit x=2^-k etwa j Glieder, mit
  45 //  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand beträgt etwa 2.8*(j/2)
  46 //  Multiplikationen von d-Bit-Zahlen. Bei Halbierungen bis x=2^-k ist der
  47 //  Gesamtaufwand etwa 2*(k+e)+1.4*j(k). Dieses minimieren nach k: Soll sein
  48 //  -1.4 = d/dk j(k) = (d/dj k(j))^-1 = - j^2/(d+j)*ln 2, also j^2=2(d+j),
  49 //  grob j=sqrt(2d) und damit k=sqrt(d).]
  50 // Aufwand: asymptotisch d^2.5 .
  51
  52         if (zerop(x))
  53                 return cl_float(1,x);
  54         var uintL d = float_digits(x);
  55         var sintL e = float_exponent(x);
  56         if (e <= (1-(sintL)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
  57                 return cl_float(1,x); // ja -> 1.0 als Ergebnis
  58  {      Mutable(cl_F,x);
  59         // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
  60         // angewandt werden. Wähle limit_slope = 13/32 = 0.4.
  61         var sintL e_limit = -1-floor(isqrt(d)*13,32); // -1-floor(sqrt(d))
  62         if (e > e_limit) {
  63                 // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
  64                 x = scale_float(x,e_limit-e);
  65                 // Neuer Exponent = e_limit.
  66         }
  67         var cl_F x2 = square(x);        // x^2
  68         // Potenzreihe anwenden:
  69         var cl_F a = x2; // a := x^2
  70         var int i = 1;
  71         var cl_F b = cl_float(1,x); // b := (float 1 x)
  72         var cl_F sum = cl_float(0,x); // sum := (float 0 x)
  73         loop {
  74                 var cl_F new_sum = sum + b;
  75                 if (new_sum == sum) // = sum ?
  76                         break; // ja -> Potenzreihe abbrechen
  77                 sum = new_sum;
  78                 b = (b*a)/(cl_I)((i+1)*(i+2));
  79                 i = i+2;
  80         }
  81         var cl_F z = square(sum); // sum^2 als Ergebnis
  82         while (e > e_limit) {
  83                 z = z + x2 * square(z);
  84                 x2 = scale_float(x2,2); // x^2 := x^2*4
  85                 e_limit++;
  86         }
  87         return z;
  88 }}
  89 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
  90
  91 const cl_LF sinhx_naive (const cl_LF& x)
  92 {
  93 // Methode:
  94 // e := Exponent aus (decode-float x), d := (float-digits x)
  95 // Bei x=0.0 oder e<=(1-d)/2 liefere x
  96 //   (denn bei e<=(1-d)/2 ist x^2/6 < x^2/4 < 2^(1-d)/4 = 2^(-d-1), also
  97 //   1 <= sinh(x)/x = 1+x^2/6+... < 1+2^(-d-1), also ist sinh(x)^2, auf d Bits
  98 //   gerundet, gleich x).
  99 // Bei e<=-sqrt(d) verwende die Potenzreihe
 100 //   sinh(x) = sum(j=0..inf,x*(x^2)^j/(2j+1)!):
 101 //   a:=x^2, b:=x, i:=1, sum:=0,
 102 //   while (/= sum (setq sum (+ sum b))) do b:=b*a/((i+1)*(i+2)), i:=i+2.
 103 //   Ergebnis sum^2.
 104 // Sonst setze y := x/2 = (scale-float x -1),
 105 //   berechne rekursiv z:=sinh(y)^2 und liefere 4*z*(1+z) = (1+2*z)^2-1.
 106 // [Die Grenze sqrt(d) ergibt sich so:
 107 //  Man braucht bei der Potenzreihe mit x=2^-k etwa j Glieder, mit
 108 //  k*j*ln 2 + j*(ln j - 1) = d, und der Aufwand beträgt etwa 2.8*(j/2)
 109 //  Multiplikationen von d-Bit-Zahlen. Bei Halbierungen bis x=2^-k ist der
 110 //  Gesamtaufwand etwa 2*(k+e)+1.4*j(k). Dieses minimieren nach k: Soll sein
 111 //  -1.4 = d/dk j(k) = (d/dj k(j))^-1 = - j^2/(d+j)*ln 2, also j^2=2(d+j),
 112 //  grob j=sqrt(2d) und damit k=sqrt(d).]
 113 // Aufwand: asymptotisch d^2.5 .
 114
 115         if (zerop(x))
 116                 return x;
 117         var uintL actuallen = TheLfloat(x)->len;
 118         var uintL d = float_digits(x);
 119         var sintL e = float_exponent(x);
 120         if (e <= (1-(sintL)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
 121                 return x; // ja -> x als Ergebnis
 122  {      Mutable(cl_LF,x);
 123         var sintL ee = e;
 124         // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
 125         // angewandt werden. Ein guter Wert für naive1 ist limit_slope = 0.6,
 126         // für naive3 aber limit_slope = 0.5.
 127         var sintL e_limit = -1-floor(isqrt(d),2); // -1-floor(sqrt(d))
 128         if (e > e_limit) {
 129                 // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
 130                 x = scale_float(x,e_limit-e);
 131                 ee = e_limit; // Neuer Exponent = e_limit.
 132         }
 133         var cl_LF x2 = square(x); // x^2
 134         // Potenzreihe anwenden:
 135         var cl_LF powser_value;
 136         var cl_LF a = x2; // a := x^2
 137         var int i = 1;
 138         if (0) {
 139                 // naive1:
 140                 // fixed-point representation
 141                 d = d-ee; // fixed-point representation with d mantissa bits
 142                 var cl_I b = round1(scale_float(x,d)); // b := x
 143                 var cl_I sum = 0; // sum := (float 0 x)
 144                 loop {
 145                         if (b == 0) break;
 146                         sum = sum + b;
 147                         b = round1(round1(The(cl_LF)(b*a)),(cl_I)((i+1)*(i+2)));
 148                         i = i+2;
 149                 }
 150                 powser_value = scale_float(cl_float(sum,x),-d);
 151         } else if (actuallen <= 7) { // Break-even-Point before extendsqrt: N<=6
 152                 // naive2:
 153                 // floating-point representation
 154                 var cl_LF b = x; // b := x
 155                 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
 156                 loop {
 157                         var cl_LF new_sum = sum + b;
 158                         if (new_sum == sum) // = sum ?
 159                                 break; // ja -> Potenzreihe abbrechen
 160                         sum = new_sum;
 161                         b = (b*a)/(cl_I)((i+1)*(i+2));
 162                         i = i+2;
 163                 }
 164                 powser_value = sum;
 165         } else {
 166                 // naive3:
 167                 // floating-point representation with smooth precision reduction
 168                 var cl_LF b = x; // b := x
 169                 var cl_LF eps = scale_float(b,-(sintL)d-10);
 170                 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
 171                 loop {
 172                         var cl_LF new_sum = sum + LF_to_LF(b,actuallen);
 173                         if (new_sum == sum) // = sum ?
 174                                 break; // ja -> Potenzreihe abbrechen
 175                         sum = new_sum;
 176                         b = cl_LF_shortenwith(b,eps);
 177                         b = (b*a)/(cl_I)((i+1)*(i+2));
 178                         i = i+2;
 179                 }
 180                 powser_value = sum;
 181         }
 182         var cl_LF z = square(powser_value); // sinh^2 als Ergebnis
 183         while (e > e_limit) {
 184                 z = square(cl_float(1,x) + scale_float(z,1)) - cl_float(1,x); // z := (1+2*z)^2-1
 185                 e_limit++;
 186         }
 187         return z;
 188 }}
 189 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
 190
 191 // Timings of the three variants, on an i486 33 MHz, running Linux,
 192 // applied to x = sqrt(2)-1 = 0.414...
 193 //   N     naive1  naive2  naive3  ratseries exp&recip
 194 //    4     0.0055  0.0039  0.0041  0.021     0.0046
 195 //    6     0.0073  0.0054  0.0054  0.029     0.0062
 196 //    8     0.0093  0.0075  0.0070  0.036     0.0081
 197 //   10     0.011   0.010   0.009   0.046     0.0011
 198 //   25     0.041   0.046   0.033   0.133     0.043
 199 //   50     0.14    0.18    0.12    0.36      0.16
 200 //  100     0.56    0.70    0.43    1.12      0.61
 201 //  250     3.5     4.5     2.7     5.3       3.3
 202 //  500    14.9    19.4    11.4    19.0      11.4
 203 // 1000    63      82      47      63        35
 204 // 2500   328     381     243     261       143
 205 // ==> naive2 fastest for N <= 6,
 206 //     naive3 fastest for 6 <= N <= 500,
 207 //     exp&recip (which uses exp's own ratseries) fastest for N >= 500.
 208
 209 }  // namespace cln