/cln-1.3.2/src/float/transcendental/cl_F_lnx.cc
C++ | 251 lines | 164 code | 13 blank | 74 comment | 31 complexity | 11fe2f6ee0eb006c847fcee58b706894 MD5 | raw file
Possible License(s): GPL-2.0
1// lnx().
2
3// General includes.
4#include "base/cl_sysdep.h"
5
6// Specification.
7#include "float/transcendental/cl_F_tran.h"
8
9
10// Implementation.
11
12#include "cln/float.h"
13#include "base/cl_low.h"
14#include "float/cl_F.h"
15#include "cln/lfloat.h"
16#include "float/lfloat/cl_LF.h"
17#include "cln/integer.h"
18
19#include "base/cl_inline.h"
20#include "float/lfloat/elem/cl_LF_zerop.cc"
21#include "float/lfloat/elem/cl_LF_minusp.cc"
22#include "float/lfloat/misc/cl_LF_exponent.cc"
23
24namespace cln {
25
26// cl_F lnx_naive (const cl_F& x)
27// cl_LF lnx_naive (const cl_LF& x)
28//
29// Methode:
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.
39// Ergebnis sum.
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 .
43
44const cl_LF lnx_naive (const cl_LF& x)
45{
46 var cl_LF y = x-cl_float(1,x);
47 if (zerop_inline(y)) // y=0.0 -> y als Ergebnis
48 return y;
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
54 { Mutable(cl_LF,x);
55 var uintL k = 0; // Rekursionszähler k:=0
56 // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
57 // angewandt werden.
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))
63 while (e > e_limit) {
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
68 k = k+1; // k:=k+1
69 }
70 if (0) {
71 // Potenzreihe ln(1+y) anwenden:
72 var int i = 1;
73 var cl_LF sum = cl_float(0,x); // sum := (float 0 x)
74 var cl_LF a = -y;
75 var cl_LF b = y;
76 if (0) {
77 // naive1:
78 // floating-point representation
79 loop {
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
83 sum = new_sum;
84 b = b*a;
85 i = i+1;
86 }
87 } else {
88 // naive2:
89 // floating-point representation with smooth precision reduction
90 var cl_LF eps = scale_float(b,-(sintC)d-10);
91 loop {
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
95 sum = new_sum;
96 b = cl_LF_shortenwith(b,eps);
97 b = b*a;
98 i = i+1;
99 }
100 }
101 return scale_float(sum,k); // sum als Ergebnis, wegen Rekursion noch mal 2^k
102 } else {
103 var cl_LF z = y / (x+cl_float(1,x));
104 // Potenzreihe atanh(z) anwenden:
105 var int i = 1;
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)
109 if (0) {
110 // naive1:
111 // floating-point representation
112 loop {
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
116 sum = new_sum;
117 b = b*a;
118 i = i+2;
119 }
120 } else {
121 // naive2:
122 // floating-point representation with smooth precision reduction
123 var cl_LF eps = scale_float(b,-(sintC)d-10);
124 loop {
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
128 sum = new_sum;
129 b = cl_LF_shortenwith(b,eps);
130 b = b*a;
131 i = i+2;
132 }
133 }
134 return scale_float(sum*z,k+1); // 2*sum*z als Ergebnis, wegen Rekursion noch mal 2^k
135 }
136}}
137// Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
138
139const cl_F lnx_naive (const cl_F& x)
140{
141 if (longfloatp(x)) {
142 DeclareType(cl_LF,x);
143 return lnx_naive(x);
144 }
145 var cl_F y = x-cl_float(1,x);
146 if (zerop(y)) // y=0.0 -> y als Ergebnis
147 return y;
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
152 { Mutable(cl_F,x);
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
161 k = k+1; // k:=k+1
162 }
163 // Potenzreihe anwenden:
164 var int i = 1;
165 var cl_F sum = cl_float(0,x); // sum := (float 0 x)
166 var cl_F a = -y;
167 var cl_F b = y;
168 loop {
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
172 sum = new_sum;
173 b = b*a;
174 i = i+1;
175 }
176 return scale_float(sum,k); // sum als Ergebnis, wegen Rekursion noch mal 2^k
177}}
178// Bit complexity (N = length(x)): O(N^(1/2)*M(N)).
179
180const cl_LF lnx_ratseries (const cl_LF& x)
181{
182 // Method:
183 // Based on the same ideas as expx_ratseries.
184 // y := 0.
185 // Loop
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;
194 { Mutable(cl_LF,x);
195 var cl_LF y = cl_I_to_LF(0,len);
196 loop {
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))
201 break;
202 var uintC lm = integer_length(x1_.mantissa);
203 var uintE me = cl_I_to_UE(- x1_.exponent);
204 var cl_I p;
205 var uintE lq;
206 var bool last_step = false;
207 if (lm >= me) { // |x'| >= 1/2 ?
208 p = x1_.sign; // 1 or -1
209 lq = 1;
210 } else {
211 var uintE n = me - lm; // |x'| < 2^-n with n maximal
212 // Set p to the first n bits of |x'|:
213 if (lm > n) {
214 p = x1_.mantissa >> (lm - n);
215 lq = 2*n;
216 } else {
217 p = x1_.mantissa;
218 lq = lm + n;
219 }
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.
226 if (2*n >= lm)
227 last_step = true;
228 }
229 y = y + scale_float(cl_I_to_LF(p,len),-(sintE)lq);
230 if (last_step)
231 break;
232 x = x * cl_exp_aux(-p,lq,len);
233 }
234 return y;
235}}
236// Bit complexity (N = length(x)): O(log(N)^2*M(N)).
237
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.)
250
251} // namespace cln