src/float/transcendental/cl_F_exp.cc

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