src/float/transcendental/cl_LF_ln10.cc

   1 // cl_ln10().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_F_tran.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/lfloat.h"
  13 #include "cl_LF.h"
  14
  15 namespace cln {
  16
  17 static inline const cl_LF compute_ln10_old (uintC len)
  18 {
  19         return ln(cl_I_to_LF(10,len));
  20 }
  21
  22 // ln 10 =
  23 // = 46 atanh(1/31) + 34 atanh(1/49) + 20 atanh(1/161)
  24 // = 478 atanh(1/251) + 180 atanh(1/449) - 126 atanh(1/4801) + 206 atanh(1/8749)
  25
  26 static inline const cl_LF compute_ln10_p235 (uintC len)
  27 {
  28         var uintC actuallen = len+1;
  29         return shorten(  The(cl_LF)(46 * cl_atanh_recip(31,actuallen))
  30                        + The(cl_LF)(34 * cl_atanh_recip(49,actuallen))
  31                        + The(cl_LF)(20 * cl_atanh_recip(161,actuallen)),
  32                        len
  33                       );
  34 }
  35
  36 static inline const cl_LF compute_ln10_p2357 (uintC len)
  37 {
  38         var uintC actuallen = len+1;
  39         return shorten(  The(cl_LF)(478 * cl_atanh_recip(251,actuallen))
  40                        + The(cl_LF)(180 * cl_atanh_recip(449,actuallen))
  41                        - The(cl_LF)(126 * cl_atanh_recip(4801,actuallen))
  42                        + The(cl_LF)(206 * cl_atanh_recip(8749,actuallen)),
  43                        len
  44                       );
  45 }
  46
  47 #define compute_ln10 compute_ln10_p2357
  48
  49 const cl_LF cl_ln10 (uintC len)
  50 {
  51         var uintC oldlen = TheLfloat(cl_LF_ln10)->len; // vorhandene Länge
  52         if (len < oldlen)
  53                 return shorten(cl_LF_ln10,len);
  54         if (len == oldlen)
  55                 return cl_LF_ln10;
  56
  57         // TheLfloat(cl_LF_ln10)->len um mindestens einen konstanten Faktor
  58         // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
  59         var uintC newlen = len;
  60         oldlen += floor(oldlen,2); // oldlen * 3/2
  61         if (newlen < oldlen)
  62                 newlen = oldlen;
  63
  64         // gewünschte > vorhandene Länge -> muß nachberechnen:
  65         cl_LF_ln10 = compute_ln10(newlen);
  66         return (len < newlen ? shorten(cl_LF_ln10,len) : cl_LF_ln10);
  67 }
  68
  69 }  // namespace cln