src/float/transcendental/cl_F_coshsinh.cc

   1 // cosh_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 cosh_sinh_t cosh_sinh (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,x)
  25 //   (denn bei e<=(1-d)/2 ist
  26 //    1 <= sinh(x)/x < cosh(x) = 1+x^2/2+... < 1+2^(-d),
  27 //    also ist cosh(x), auf d Bits gerundet, gleich 1.0
  28 //    und sinh(x), auf d Bits gerundet, gleich x).
  29 // falls e<0:
  30 //   y:=(sinh(x)/x)^2 errechnen,
  31 //   cosh(x) = sqrt(1+x^2*y) und sinh(x) = x*sqrt(y) errechnen.
  32 // falls e>=0: y:=exp(x) errechnen,
  33 //   (scale-float (+ y (/ y)) -1) und (scale-float (- y (/ y)) -1) bilden.
  34 // Genauigkeit wieder verringern.
  35
  36         var sintE e = float_exponent(x);
  37         if (e < 0) { // Exponent e abtesten
  38                 // e<0
  39                 if (zerop(x) || (e <= (1-(sintC)float_digits(x))>>1))
  40                         // e <= (1-d)/2 <==> e <= -ceiling((d-1)/2)
  41                         return cosh_sinh_t(cl_float(1,x),x);
  42                 // Rechengenauigkeit erhöhen
  43                 if (longfloatp(x)) {
  44                         DeclareType(cl_LF,x);
  45                         #if 0
  46                         if (TheLfloat(x)->len >= infty) {
  47                                 var cl_LF xx = extend(x,TheLfloat(x)->len+1);
  48                                 var cl_LF_cosh_sinh_t hyp = cl_coshsinh_ratseries(xx);
  49                                 return cosh_sinh_t(
  50                                         cl_float(hyp.cosh,x),
  51                                         cl_float(hyp.sinh,x)
  52                                        );
  53                         } else
  54                         #endif
  55                         if (TheLfloat(x)->len >= 585) {
  56                                 // verwende exp(x), schneller als cl_coshsinh_ratseries
  57                                 var cl_LF xx = extend(x,TheLfloat(x)->len+ceiling((uintE)(-e),intDsize));
  58                                 var cl_F y = exp(xx);
  59                                 var cl_F y_inv = recip(y);
  60                                 return cosh_sinh_t(
  61                                         cl_float(scale_float(y + y_inv, -1), x),
  62                                         cl_float(scale_float(y - y_inv, -1), x)
  63                                        );
  64                         } else {
  65                                 var cl_LF xx = The(cl_LF)(cl_F_extendsqrt(x));
  66                                 var cl_LF y = sinhx_naive(xx);
  67                                 var cl_LF z = sqrt(y);
  68                                 if (minusp(xx))
  69                                         z = -z;
  70                                 return cosh_sinh_t(
  71                                         cl_float(sqrt(1+y),x), // sqrt(1+y)
  72                                         cl_float(z,x)
  73                                        );
  74                         }
  75                 } else {
  76                         var cl_F xx = cl_F_extendsqrt(x);
  77                         var cl_F y = sinhxbyx_naive(xx);
  78                         return cosh_sinh_t(
  79                                 cl_float(sqrt(1+square(xx)*y),x), // sqrt(1+x^2*y)
  80                                 cl_float(xx*sqrt(y),x)
  81                                );
  82                 }
  83         } else {
  84                 // e>=0 -> verwende exp(x)
  85                 var cl_F y = exp(x);
  86                 var cl_F y_inv = recip(y);
  87                 return cosh_sinh_t(
  88                         scale_float(y+y_inv,-1),
  89                         scale_float(y-y_inv,-1)
  90                        );
  91         }
  92 }
  93
  94 }  // namespace cln