12 #include "cln/lfloat.h"
17 static inline const cl_LF compute_ln2_old (uintC len)
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);
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)
31 static inline const cl_LF compute_ln2_p2 (uintC len)
33 return scale_float(cl_atanh_recip(3,len),1);
36 static inline const cl_LF compute_ln2_p23 (uintC len)
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),
45 static inline const cl_LF compute_ln2_p235 (uintC len)
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)),
55 static inline const cl_LF compute_ln2_p2357 (uintC len)
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)),
66 // Timings of the above algorithms, on an i486 33 MHz, running Linux.
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).
75 #define compute_ln2 compute_ln2_p2357
77 const cl_LF cl_ln2 (uintC len)
79 var uintC oldlen = TheLfloat(cl_LF_ln2)->len; // vorhandene Länge
81 return shorten(cl_LF_ln2,len);
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
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);