src/float/transcendental/cl_F_exp.cc

   1 // exp().
   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_float.h"
  14 #include "cl_F.h"
  15 #include "cl_lfloat.h"
  16 #include "cl_LF.h"
  17
  18 #undef MAYBE_INLINE
  19 #define MAYBE_INLINE inline
  20 #include "cl_LF_minusp.cc"
  21 #include "cl_LF_exponent.cc"
  22
  23 // Division durch ln(2).
  24 inline const cl_F_div_t cl_floor_ln2 (const cl_F& x)
  25 {
  26         // Bei 0<=x<1/2 kann man sofort q:=0 setzen.
  27         if (!minusp(x) && (float_exponent(x) < 0))
  28                 return cl_F_div_t(0,x);
  29         else
  30                 return floor2(x,cl_ln2(x));
  31 }
  32 inline const cl_LF_div_t cl_floor_ln2 (const cl_LF& x)
  33 {
  34         // Bei 0<=x<1/2 kann man sofort q:=0 setzen.
  35         if (!minusp(x) && (float_exponent(x) < 0))
  36                 return cl_LF_div_t(0,x);
  37         else
  38                 return floor2(x,The(cl_LF)(cl_ln2(x)));
  39 }
  40
  41 const cl_F exp (const cl_F& x)
  42 {
  43 // Methode:
  44 // d := (float-digits x),
  45 // Genauigkeit um sqrt(d)+max(integer-length(e)) Bits erhöhen,
  46 // (q,r) := (floor x ln(2))
  47 // Ergebnis ist exp(q*ln(2)+r) = (scale-float exp(r) q).
  48
  49         // Rechengenauigkeit erhöhen und durch ln(2) dividieren:
  50         if (longfloatp(x) && (TheLfloat(x)->len >= 84)) {
  51                 DeclareType(cl_LF,x);
  52                 var cl_LF_div_t q_r = cl_floor_ln2(extend(x,TheLfloat(x)->len+1));
  53                 var cl_I& q = q_r.quotient;
  54                 var cl_LF& r = q_r.remainder;
  55                 return cl_float(scale_float(expx_ratseries(r),q),x);
  56         } else {
  57                 var cl_F_div_t q_r = cl_floor_ln2(cl_F_extendsqrtx(x));
  58                 var cl_I& q = q_r.quotient;
  59                 var cl_F& r = q_r.remainder;
  60                 return cl_float(scale_float(expx_naive(r),q),x);
  61         }
  62 }