src/float/transcendental/cl_F_sinh.cc

   1 // sinh().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/float.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_F_tran.h"
  13 #include "cl_F.h"
  14 #include "cln/lfloat.h"
  15 #include "cl_LF.h"
  16
  17 namespace cln {
  18
  19 const cl_F sinh (const cl_F& x)
  20 {
  21 // Methode:
  22 // Genauigkeit erhöhen,
  23 // e := Exponent aus (decode-float x)
  24 // falls e<0: (sinh(x)/x)^2 errechnen, Wurzel ziehen, mit x multiplizieren.
  25 // falls e>=0: y:=exp(x) errechnen, (scale-float (- y (/ y)) -1) bilden.
  26
  27         if (float_exponent(x) < 0) { // Exponent e abtesten
  28                 // e<0
  29                 // Rechengenauigkeit erhöhen
  30                 if (longfloatp(x)) {
  31                         DeclareType(cl_LF,x);
  32                         #if 0
  33                         if (TheLfloat(x)->len >= infty) {
  34                                 var cl_LF xx = extend(x,TheLfloat(x)->len+1);
  35                                 var cl_LF_cosh_sinh_t hyp = cl_coshsinh_ratseries(xx);
  36                                 return cln/float.hyp.sinh,x);
  37                         } else
  38                         #endif
  39                         if ((TheLfloat(x)->len >= 500)
  40                             && (float_exponent(x) > (-(sintC)float_digits(x))>>1)) {
  41                                 // verwende exp(x), schneller als cl_coshsinh_ratseries
  42                                 // (aber nur bei 0 > e > -d/2, denn wir müssen, um
  43                                 // Auslöschung zu verhindern, |e| Bits dazunehmen)
  44                                 var cl_LF xx = extend(x,TheLfloat(x)->len+ceiling((uintE)(-float_exponent(x)),intDsize));
  45                                 var cl_F y = exp(xx);
  46                                 var cl_F z = scale_float(y - recip(y), -1); // (/ (- y (/ y)) 2)
  47                                 return cl_float(z,x);
  48                         } else {
  49                                 var cl_LF xx = The(cl_LF)(cl_F_extendsqrt(x));
  50                                 // Wurzel aus sinh(x)^2 bilden
  51                                 var cl_LF z = sqrt(sinhx_naive(xx));
  52                                 if (minusp(xx))
  53                                         z = -z;
  54                                 return cl_float(z,x);
  55                         }
  56                 } else {
  57                         var cl_F xx = cl_F_extendsqrt(x);
  58                         // Wurzel aus (sinh(x)/x)^2 mit x multiplizieren und wieder runden
  59                         return cl_float(sqrt(sinhxbyx_naive(xx))*xx,x);
  60                 }
  61         } else {
  62                 // e>=0 -> verwende exp(x)
  63                 var cl_F y = exp(x);
  64                 return scale_float(y - recip(y), -1); // (/ (- y (/ y)) 2)
  65         }
  66 }
  67
  68 // Timings of the two algorithms, on an i486 33 MHz, running Linux,
  69 // applied to x = sqrt(2)-1 = 0.414...
  70 //   N      naive  ratseries
  71 //   10     0.008   0.037
  72 //   25     0.034   0.115
  73 //   50     0.13    0.33
  74 //  100     0.50    1.07
  75 //  250     3.3     5.2
  76 //  500    14.2    18.8
  77 // 1000    59      61
  78 // 2500   297     247
  79 // ==> ratseries faster for N >= 1300.
  80
  81 }  // namespace cln