src/float/transcendental/cl_F_cosh.cc

   1 // cosh().
   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 cosh (const cl_F& x)
  20 {
  21 // Methode:
  22 // Genauigkeit erhöhen,
  23 // e := Exponent aus (decode-float x), d := (float-digits x)
  24 // falls x=0.0 oder e<=(1-d)/2 liefere 1.0
  25 //   (denn bei e<=(1-d)/2 ist 1 <= cosh(x) = 1+x^2/2+... < 1+2^(-d),
  26 //    also ist cosh(x), auf d Bits gerundet, gleich 1.0).
  27 // falls e<0:
  28 //   y := x/2 = (scale-float x -1), (sinh(y)/y)^2 errechnen,
  29 //   cosh(x) = 1+x*y*(sinh(y)/y)^2 errechnen.
  30 // falls e>=0: y:=exp(x) errechnen, (scale-float (+ y (/ y)) -1) bilden.
  31
  32         var sintL e = float_exponent(x);
  33         if (e < 0) { // Exponent e abtesten
  34                 // e<0
  35                 if (zerop(x))
  36                         return cl_float(1,x);
  37                 var uintL d = float_digits(x);
  38                 if (e <= (1-(sintL)d)>>1) // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2) ?
  39                         return cl_float(1,x); // ja -> 1.0 als Ergebnis
  40                 // Rechengenauigkeit erhöhen
  41                 if (longfloatp(x)) {
  42                         DeclareType(cl_LF,x);
  43                         #if 0
  44                         if (TheLfloat(x)->len >= infty) {
  45                                 var cl_LF xx = extend(x,TheLfloat(x)->len+1);
  46                                 var cl_LF_cosh_sinh_t hyp = cl_coshsinh_ratseries(xx);
  47                                 return cln/float.hyp.cosh,x);
  48                         } else
  49                         #endif
  50                         if (TheLfloat(x)->len >= 600) {
  51                                 // verwende exp(x), schneller als cl_coshsinh_ratseries
  52                                 var cl_LF xx = extend(x,TheLfloat(x)->len+1);
  53                                 var cl_F y = exp(xx);
  54                                 var cl_F z = scale_float(y + recip(y), -1); // (/ (+ y (/ y)) 2)
  55                                 return cl_float(z,x);
  56                         } else {
  57                                 var cl_LF xx = The(cl_LF)(cl_F_extendsqrt(x));
  58                                 var cl_LF y = scale_float(xx,-1);
  59                                 // 1 + 2*sinh(y)^2, und wieder runden
  60                                 return cl_float(1 + scale_float(sinhx_naive(y),1), x);
  61                         }
  62                 } else {
  63                         var cl_F xx = cl_F_extendsqrt(x);
  64                         var cl_F y = scale_float(xx,-1);
  65                         // 1 + 2*y^2*(sinh(y)/y)^2, und wieder runden
  66                         return cl_float(1 + scale_float(square(y) * sinhxbyx_naive(y),1), x);
  67                 }
  68         } else {
  69                 // e>=0 -> verwende exp(x)
  70                 var cl_F y = exp(x);
  71                 return scale_float(y + recip(y), -1); // (/ (+ y (/ y)) 2)
  72         }
  73 }
  74
  75 // Timings of the three algorithms, on an i486 33 MHz, running Linux,
  76 // applied to x = sqrt(2)-1 = 0.414...
  77 //   N      naive  ratseries exp&recip
  78 //   10     0.008   0.037     0.012
  79 //   25     0.032   0.117     0.047
  80 //   50     0.11    0.33      0.017
  81 //  100     0.40    1.06      0.63
  82 //  250     2.65    5.2       3.3
  83 //  500    11.1    18.7      11.5
  84 // 1000    46      61        35
  85 // 2500   238     250       143
  86 // ==> exp&recip fastest for N >= 600.
  87
  88 }  // namespace cln