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