src/float/transcendental/cl_F_ln.cc

   1 // ln().
   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 "cl_SF.h"
  15 #include "cln/integer.h"
  16 #include "cln/lfloat.h"
  17 #include "cl_LF.h"
  18
  19 namespace cln {
  20
  21 const cl_F ln (const cl_F& x)
  22 {
  23 // Methode:
  24 // d := (float-digits x),
  25 // Genauigkeit um sqrt(d)+max(integer-length(e)) Bits erhöhen,
  26 // (m,e) := (decode-float x), so daß 1/2 <= m < 1.
  27 // m<2/3 -> m:=2m, e:=e-1, so daß 2/3 <= m <= 4/3.
  28 // ln(m) errechnen, ln(x)=ln(m)+e*ln(2) als Ergebnis.
  29
  30         // Rechengenauigkeit erhöhen und m,e,s bestimmen:
  31         if (longfloatp(x) && (TheLfloat(x)->len >= 110)) {
  32                 DeclareType(cl_LF,x);
  33                 var decoded_lfloat m_e_s = decode_float(extend(x,TheLfloat(x)->len+1));
  34                 var cl_LF& m = m_e_s.mantissa;
  35                 var cl_I& e = m_e_s.exponent;
  36                 if (m < make_SF(0,0+SF_exp_mid,floor(bit(SF_mant_len+2),3))) { // Short-Float 2/3
  37                         m = scale_float(m,1); // m verdoppeln
  38                         e = minus1(e); // e decrementieren
  39                 }
  40                 var cl_F res = lnx_ratseries(m);
  41                 if (!zerop(e))
  42                         res = res + cl_float(e,m)*cl_ln2(m); // ln(m)+e*ln(2)
  43                 return cl_float(res,x);
  44         } else {
  45                 var decoded_float m_e_s = decode_float(cl_F_extendsqrtx(x));
  46                 var cl_F& m = m_e_s.mantissa;
  47                 var cl_I& e = m_e_s.exponent;
  48                 if (m < make_SF(0,0+SF_exp_mid,floor(bit(SF_mant_len+2),3))) { // Short-Float 2/3
  49                         m = scale_float(m,1); // m verdoppeln
  50                         e = minus1(e); // e decrementieren
  51                 }
  52                 var cl_F res = lnx_naive(m);
  53                 if (!zerop(e))
  54                         res = res + cl_float(e,m)*cl_ln2(m); // ln(m)+e*ln(2)
  55                 return cl_float(res,x);
  56         }
  57 }
  58
  59 }  // namespace cln