src/float/transcendental/cl_LF_ln2.cc

   1 // cl_ln2().
   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_ln2_old (uintC len)
  18 {
  19         // Here, it is tricky to avoid a recursive loop. We assume that ln()
  20         // will not invoke cl_ln2() if its argument is between 2/3 and 4/3.
  21         // So we compute -2*ln(1/sqrt(2)).
  22         return -scale_float(ln(sqrt(scale_float(cl_I_to_LF(1,len),-1))),1);
  23 }
  24
  25 // ln 2 =
  26 //  = 2 atanh(1/3)
  27 //  = 4 atanh(1/7) + 2 atanh(1/17)
  28 //  = 14 atanh(1/31) + 10 atanh(1/49) + 6 atanh(1/161)
  29 //  = 144 atanh(1/251) + 54 atanh(1/449) - 38 atanh(1/4801) + 62 atanh(1/8749)
  30
  31 static inline const cl_LF compute_ln2_p2 (uintC len)
  32 {
  33         return scale_float(cl_atanh_recip(3,len),1);
  34 }
  35
  36 static inline const cl_LF compute_ln2_p23 (uintC len)
  37 {
  38         var uintC actuallen = len+1;
  39         return shorten(scale_float(cl_atanh_recip(7,actuallen),2)
  40                        + scale_float(cl_atanh_recip(17,actuallen),1),
  41                        len
  42                       );
  43 }
  44
  45 static inline const cl_LF compute_ln2_p235 (uintC len)
  46 {
  47         var uintC actuallen = len+1;
  48         return shorten(  The(cl_LF)(14 * cl_atanh_recip(31,actuallen))
  49                        + The(cl_LF)(10 * cl_atanh_recip(49,actuallen))
  50                        + The(cl_LF)( 6 * cl_atanh_recip(161,actuallen)),
  51                        len
  52                       );
  53 }
  54
  55 static inline const cl_LF compute_ln2_p2357 (uintC len)
  56 {
  57         var uintC actuallen = len+1;
  58         return shorten(  The(cl_LF)(144 * cl_atanh_recip(251,actuallen))
  59                        + The(cl_LF)( 54 * cl_atanh_recip(449,actuallen))
  60                        - The(cl_LF)( 38 * cl_atanh_recip(4801,actuallen))
  61                        + The(cl_LF)( 62 * cl_atanh_recip(8749,actuallen)),
  62                        len
  63                       );
  64 }
  65
  66 // Timings of the above algorithms, on an i486 33 MHz, running Linux.
  67 // N = 250
  68 // p2       1.71 s                        total 1.71 s
  69 // p23      0.88 + 0.60 s                 total 1.48 s
  70 // p235     0.51 + 0.48 + 0.40 s          total 1.39 s
  71 // p2357    0.38 + 0.37 + 0.31 + 0.29 s   total 1.35 s
  72 // In general, for fixed precision N, the computation time of atanh(1/m)
  73 // seems to be proportional to 1/log(m).
  74
  75 #define compute_ln2 compute_ln2_p2357
  76
  77 const cl_LF cl_ln2 (uintC len)
  78 {
  79         var uintC oldlen = TheLfloat(cl_LF_ln2)->len; // vorhandene Länge
  80         if (len < oldlen)
  81                 return shorten(cl_LF_ln2,len);
  82         if (len == oldlen)
  83                 return cl_LF_ln2;
  84
  85         // TheLfloat(cl_LF_ln2)->len um mindestens einen konstanten Faktor
  86         // > 1 wachsen lassen, damit es nicht zu häufig nachberechnet wird:
  87         var uintC newlen = len;
  88         oldlen += floor(oldlen,2); // oldlen * 3/2
  89         if (newlen < oldlen)
  90                 newlen = oldlen;
  91
  92         // gewünschte > vorhandene Länge -> muß nachberechnen:
  93         cl_LF_ln2 = compute_ln2(newlen);
  94         return (len < newlen ? shorten(cl_LF_ln2,len) : cl_LF_ln2);
  95 }
  96
  97 }  // namespace cln