12 #include "cln/float.h"
15 #include "cln/lfloat.h"
17 #include "cln/integer.h"
20 #define MAYBE_INLINE inline
21 #include "cl_LF_zerop.cc"
22 #include "cl_LF_minusp.cc"
23 #include "cl_LF_exponent.cc"
27 // cl_F lnx_naive (const cl_F& x)
28 // cl_LF lnx_naive (const cl_LF& x)
31 // y:=x-1, e := Exponent aus (decode-float y), d := (float-digits y)
32 // Bei y=0.0 oder e<=-d liefere y
33 // (denn bei e<=-d ist y/2 < 2^(-d)/2 = 2^(-d-1), also
34 // 0 <= y - ln(x) < y^2/2 < 2^(-d-1)*y
35 // also ist ln(x)/y, auf d Bits gerundet, gleich y).
36 // Bei e<=-sqrt(d) verwende die Potenzreihe
37 // ln(x) = sum(j=0..inf,(-1)^j*y^(j+1)/(j+1)):
38 // a:=-y, b:=y, i:=1, sum:=0,
39 // while (/= sum (setq sum (+ sum (/ b i)))) do i:=i+1, b:=b*a.
41 // Sonst setze y := sqrt(x), berechne rekursiv z:=ln(y)
42 // und liefere 2*z = (scale-float z 1).
43 // Aufwand: asymptotisch d^0.5*M(d) = d^2.5 .
45 const cl_LF lnx_naive (const cl_LF& x)
47 var cl_LF y = x-cl_float(1,x);
48 if (zerop(y)) // y=0.0 -> y als Ergebnis
50 var uintL actuallen = TheLfloat(x)->len;
51 var uintC d = float_digits(x);
52 var sintE e = float_exponent(y);
53 if (e <= -(sintC)d) // e <= -d ?
54 return y; // ja -> y als Ergebnis
56 var uintL k = 0; // Rekursionszähler k:=0
57 // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
59 // Wähle für ln(1+y), naive1: limit_slope = 1.0,
60 // für ln(1+y), naive2: limit_slope = 11/16 = 0.7,
61 // für atanh(z), naive1: limit_slope = 0.6,
62 // für atanh(z), naive1: limit_slope = 0.5.
63 var sintL e_limit = -1-floor(isqrtC(d),2); // -1-floor(sqrt(d))
65 // e > -1-floor(sqrt(d)) -> muß |y| verkleinern.
66 x = sqrt(x); // x := (sqrt x)
67 y = x-cl_float(1,x); // y := (- x 1) und
68 e = float_exponent(y); // e neu berechnen
72 // Potenzreihe ln(1+y) anwenden:
74 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
79 // floating-point representation
81 var cl_LF new_sum = sum + b/(cl_I)i; // (+ sum (/ b i))
82 if (new_sum == sum) // = sum ?
83 break; // ja -> Potenzreihe abbrechen
90 // floating-point representation with smooth precision reduction
91 var cl_LF eps = scale_float(b,-(sintC)d-10);
93 var cl_LF new_sum = sum + LF_to_LF(b/(cl_I)i,actuallen); // (+ sum (/ b i))
94 if (new_sum == sum) // = sum ?
95 break; // ja -> Potenzreihe abbrechen
97 b = cl_LF_shortenwith(b,eps);
102 return scale_float(sum,k); // sum als Ergebnis, wegen Rekursion noch mal 2^k
104 var cl_LF z = y / (x+cl_float(1,x));
105 // Potenzreihe atanh(z) anwenden:
107 var cl_LF a = square(z); // a = x^2
108 var cl_LF b = cl_float(1,x); // b := (float 1 x)
109 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
112 // floating-point representation
114 var cl_LF new_sum = sum + b / (cl_I)i; // (+ sum (/ b i))
115 if (new_sum == sum) // = sum ?
116 break; // ja -> Potenzreihe abbrechen
123 // floating-point representation with smooth precision reduction
124 var cl_LF eps = scale_float(b,-(sintC)d-10);
126 var cl_LF new_sum = sum + LF_to_LF(b/(cl_I)i,actuallen); // (+ sum (/ b i))
127 if (new_sum == sum) // = sum ?
128 break; // ja -> Potenzreihe abbrechen
130 b = cl_LF_shortenwith(b,eps);
135 return scale_float(sum*z,k+1); // 2*sum*z als Ergebnis, wegen Rekursion noch mal 2^k
138 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
140 const cl_F lnx_naive (const cl_F& x)
143 DeclareType(cl_LF,x);
146 var cl_F y = x-cl_float(1,x);
147 if (zerop(y)) // y=0.0 -> y als Ergebnis
149 var uintC d = float_digits(x);
150 var sintE e = float_exponent(y);
151 if (e <= -(sintC)d) // e <= -d ?
152 return y; // ja -> y als Ergebnis
154 var uintL k = 0; // Rekursionszähler k:=0
155 // Bei e <= -1-floor(sqrt(d)) kann die Potenzreihe angewandt werden.
156 var sintL e_limit = -1-isqrtC(d); // -1-floor(sqrt(d))
157 while (e > e_limit) {
158 // e > -1-floor(sqrt(d)) -> muß |y| verkleinern.
159 x = sqrt(x); // x := (sqrt x)
160 y = x-cl_float(1,x); // y := (- x 1) und
161 e = float_exponent(y); // e neu berechnen
164 // Potenzreihe anwenden:
166 var cl_F sum = cl_float(0,x); // sum := (float 0 x)
170 var cl_F new_sum = sum + b/(cl_I)i; // (+ sum (/ b i))
171 if (new_sum == sum) // = sum ?
172 break; // ja -> Potenzreihe abbrechen
177 return scale_float(sum,k); // sum als Ergebnis, wegen Rekursion noch mal 2^k
179 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
181 const cl_LF lnx_ratseries (const cl_LF& x)
184 // Based on the same ideas as expx_ratseries.
187 // [x*exp(y) is invariant]
188 // x' := x-1. If x' = 0, terminate the loop.
189 // Choose approximation y' of log(x) = log(1+x'):
190 // If |x'| >= 1/2, set y' = 1/2 * sign(x').
191 // If |x'| < 2^-n with n maximal, set
192 // y' = truncate(x'*2^(2n))/2^(2n).
193 // Set y := y + y' and x := x*exp(-y').
194 var uintC len = TheLfloat(x)->len;
196 var cl_LF y = cl_I_to_LF(0,len);
198 var cl_LF x1 = x + cl_I_to_LF(-1,len);
199 var cl_idecoded_float x1_ = integer_decode_float(x1);
200 // x1 = (-1)^sign * 2^exponent * mantissa
201 if (zerop(x1_.mantissa))
203 var uintC lm = integer_length(x1_.mantissa);
204 var uintE me = cl_I_to_UE(- x1_.exponent);
207 var bool last_step = false;
208 if (lm >= me) { // |x'| >= 1/2 ?
209 p = x1_.sign; // 1 or -1
212 var uintE n = me - lm; // |x'| < 2^-n with n maximal
213 // Set p to the first n bits of |x'|:
215 p = x1_.mantissa >> (lm - n);
221 if (minusp(x1_.sign)) { p = -p; }
222 // If 2*n >= lm = intDsize*len, then within our
223 // precision exp(-y') = 1-y', (because |y'^2| < 2^-lm),
224 // and we know a priori that the iteration will stop
225 // after the next big multiplication. This saves one
226 // big multiplication at the end.
230 y = y + scale_float(cl_I_to_LF(p,len),-(sintE)lq);
233 x = x * cl_exp_aux(-p,lq,len);
237 // Bit complexity (N = length(x)): O(log(N)^2*M(N)).
239 // Timings of the above algorithms, on an i486 33 MHz, running Linux,
240 // applied to x = sqrt(sqrt(2)) = 1.189...
241 // N ln(1+y) ln(1+y) atanh z atanh z exp
242 // naive1 naive2 naive1 naive2 ratseries
243 // 10 0.019 0.016 0.013 0.012 0.036
244 // 25 0.077 0.056 0.057 0.040 0.087
245 // 50 0.30 0.21 0.23 0.15 0.21
246 // 100 1.24 0.81 0.92 0.59 0.61
247 // 250 8.8 5.8 6.3 4.3 2.77
248 // 500 43.9 28.8 29.7 21.0 9.8
249 // 1000 223 149 144 107 30
250 // ==> ratseries faster for N >= 110. (N = length before extended by the caller.)