12 #include "cln/float.h"
15 #include "cln/lfloat.h"
17 #include "cln/integer.h"
19 #include "cl_inline.h"
20 #include "cl_LF_zerop.cc"
21 #include "cl_LF_minusp.cc"
22 #include "cl_LF_exponent.cc"
26 // cl_F lnx_naive (const cl_F& x)
27 // cl_LF lnx_naive (const cl_LF& x)
30 // y:=x-1, e := Exponent aus (decode-float y), d := (float-digits y)
31 // Bei y=0.0 oder e<=-d liefere y
32 // (denn bei e<=-d ist y/2 < 2^(-d)/2 = 2^(-d-1), also
33 // 0 <= y - ln(x) < y^2/2 < 2^(-d-1)*y
34 // also ist ln(x)/y, auf d Bits gerundet, gleich y).
35 // Bei e<=-sqrt(d) verwende die Potenzreihe
36 // ln(x) = sum(j=0..inf,(-1)^j*y^(j+1)/(j+1)):
37 // a:=-y, b:=y, i:=1, sum:=0,
38 // while (/= sum (setq sum (+ sum (/ b i)))) do i:=i+1, b:=b*a.
40 // Sonst setze y := sqrt(x), berechne rekursiv z:=ln(y)
41 // und liefere 2*z = (scale-float z 1).
42 // Aufwand: asymptotisch d^0.5*M(d) = d^2.5 .
44 const cl_LF lnx_naive (const cl_LF& x)
46 var cl_LF y = x-cl_float(1,x);
47 if (zerop_inline(y)) // y=0.0 -> y als Ergebnis
49 var uintC actuallen = TheLfloat(x)->len;
50 var uintC d = float_digits(x);
51 var sintE e = float_exponent_inline(y);
52 if (e <= -(sintC)d) // e <= -d ?
53 return y; // ja -> y als Ergebnis
55 var uintL k = 0; // Rekursionszähler k:=0
56 // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
58 // Wähle für ln(1+y), naive1: limit_slope = 1.0,
59 // für ln(1+y), naive2: limit_slope = 11/16 = 0.7,
60 // für atanh(z), naive1: limit_slope = 0.6,
61 // für atanh(z), naive1: limit_slope = 0.5.
62 var sintL e_limit = -1-floor(isqrtC(d),2); // -1-floor(sqrt(d))
64 // e > -1-floor(sqrt(d)) -> muß |y| verkleinern.
65 x = sqrt(x); // x := (sqrt x)
66 y = x-cl_float(1,x); // y := (- x 1) und
67 e = float_exponent_inline(y); // e neu berechnen
71 // Potenzreihe ln(1+y) anwenden:
73 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
78 // floating-point representation
80 var cl_LF new_sum = sum + b/(cl_I)i; // (+ sum (/ b i))
81 if (new_sum == sum) // = sum ?
82 break; // ja -> Potenzreihe abbrechen
89 // floating-point representation with smooth precision reduction
90 var cl_LF eps = scale_float(b,-(sintC)d-10);
92 var cl_LF new_sum = sum + LF_to_LF(b/(cl_I)i,actuallen); // (+ sum (/ b i))
93 if (new_sum == sum) // = sum ?
94 break; // ja -> Potenzreihe abbrechen
96 b = cl_LF_shortenwith(b,eps);
101 return scale_float(sum,k); // sum als Ergebnis, wegen Rekursion noch mal 2^k
103 var cl_LF z = y / (x+cl_float(1,x));
104 // Potenzreihe atanh(z) anwenden:
106 var cl_LF a = square(z); // a = x^2
107 var cl_LF b = cl_float(1,x); // b := (float 1 x)
108 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
111 // floating-point representation
113 var cl_LF new_sum = sum + b / (cl_I)i; // (+ sum (/ b i))
114 if (new_sum == sum) // = sum ?
115 break; // ja -> Potenzreihe abbrechen
122 // floating-point representation with smooth precision reduction
123 var cl_LF eps = scale_float(b,-(sintC)d-10);
125 var cl_LF new_sum = sum + LF_to_LF(b/(cl_I)i,actuallen); // (+ sum (/ b i))
126 if (new_sum == sum) // = sum ?
127 break; // ja -> Potenzreihe abbrechen
129 b = cl_LF_shortenwith(b,eps);
134 return scale_float(sum*z,k+1); // 2*sum*z als Ergebnis, wegen Rekursion noch mal 2^k
137 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
139 const cl_F lnx_naive (const cl_F& x)
142 DeclareType(cl_LF,x);
145 var cl_F y = x-cl_float(1,x);
146 if (zerop(y)) // y=0.0 -> y als Ergebnis
148 var uintC d = float_digits(x);
149 var sintE e = float_exponent(y);
150 if (e <= -(sintC)d) // e <= -d ?
151 return y; // ja -> y als Ergebnis
153 var uintL k = 0; // Rekursionszähler k:=0
154 // Bei e <= -1-floor(sqrt(d)) kann die Potenzreihe angewandt werden.
155 var sintL e_limit = -1-isqrtC(d); // -1-floor(sqrt(d))
156 while (e > e_limit) {
157 // e > -1-floor(sqrt(d)) -> muß |y| verkleinern.
158 x = sqrt(x); // x := (sqrt x)
159 y = x-cl_float(1,x); // y := (- x 1) und
160 e = float_exponent(y); // e neu berechnen
163 // Potenzreihe anwenden:
165 var cl_F sum = cl_float(0,x); // sum := (float 0 x)
169 var cl_F new_sum = sum + b/(cl_I)i; // (+ sum (/ b i))
170 if (new_sum == sum) // = sum ?
171 break; // ja -> Potenzreihe abbrechen
176 return scale_float(sum,k); // sum als Ergebnis, wegen Rekursion noch mal 2^k
178 // Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
180 const cl_LF lnx_ratseries (const cl_LF& x)
183 // Based on the same ideas as expx_ratseries.
186 // [x*exp(y) is invariant]
187 // x' := x-1. If x' = 0, terminate the loop.
188 // Choose approximation y' of log(x) = log(1+x'):
189 // If |x'| >= 1/2, set y' = 1/2 * sign(x').
190 // If |x'| < 2^-n with n maximal, set
191 // y' = truncate(x'*2^(2n))/2^(2n).
192 // Set y := y + y' and x := x*exp(-y').
193 var uintC len = TheLfloat(x)->len;
195 var cl_LF y = cl_I_to_LF(0,len);
197 var cl_LF x1 = x + cl_I_to_LF(-1,len);
198 var cl_idecoded_float x1_ = integer_decode_float(x1);
199 // x1 = (-1)^sign * 2^exponent * mantissa
200 if (zerop(x1_.mantissa))
202 var uintC lm = integer_length(x1_.mantissa);
203 var uintE me = cl_I_to_UE(- x1_.exponent);
206 var bool last_step = false;
207 if (lm >= me) { // |x'| >= 1/2 ?
208 p = x1_.sign; // 1 or -1
211 var uintE n = me - lm; // |x'| < 2^-n with n maximal
212 // Set p to the first n bits of |x'|:
214 p = x1_.mantissa >> (lm - n);
220 if (minusp(x1_.sign)) { p = -p; }
221 // If 2*n >= lm = intDsize*len, then within our
222 // precision exp(-y') = 1-y', (because |y'^2| < 2^-lm),
223 // and we know a priori that the iteration will stop
224 // after the next big multiplication. This saves one
225 // big multiplication at the end.
229 y = y + scale_float(cl_I_to_LF(p,len),-(sintE)lq);
232 x = x * cl_exp_aux(-p,lq,len);
236 // Bit complexity (N = length(x)): O(log(N)^2*M(N)).
238 // Timings of the above algorithms, on an i486 33 MHz, running Linux,
239 // applied to x = sqrt(sqrt(2)) = 1.189...
240 // N ln(1+y) ln(1+y) atanh z atanh z exp
241 // naive1 naive2 naive1 naive2 ratseries
242 // 10 0.019 0.016 0.013 0.012 0.036
243 // 25 0.077 0.056 0.057 0.040 0.087
244 // 50 0.30 0.21 0.23 0.15 0.21
245 // 100 1.24 0.81 0.92 0.59 0.61
246 // 250 8.8 5.8 6.3 4.3 2.77
247 // 500 43.9 28.8 29.7 21.0 9.8
248 // 1000 223 149 144 107 30
249 // ==> ratseries faster for N >= 110. (N = length before extended by the caller.)