src/float/transcendental/cl_F_cosh.cc

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