src/float/transcendental/cl_F_coshsinh.cc

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