9 #include "cln/lfloat.h"
15 #include "cln/float.h"
18 #include "cl_inline.h"
19 #include "cl_LF_zerop.cc"
20 #include "cl_LF_minusp.cc"
21 #include "cl_LF_exponent.cc"
25 // cl_F atanhx (const cl_F& x)
26 // cl_LF atanhx (const cl_LF& x)
29 // e := Exponent aus (decode-float x), d := (float-digits x)
30 // Bei x=0.0 oder e<=-d/2 liefere x
31 // (denn bei e<=-d/2 ist x^2 < 2^(-d), also
32 // 1 <= atanh(x)/x = 1+x^2/3+x^4/5+... < 1+x^2/2 < 1+2^(-d-1) < 1+2^(-d),
33 // also ist atanh(x)/x, auf d Bits gerundet, gleich 1.0).
34 // Bei großem d verwende die Formel ln((1+x)/(1-x))/2 (asymptotisch schneller),
35 // aber erhöhe die Genauigkeit, so daß beim Bilden von 1+x keine Bits verloren
37 // Bei e<=-sqrt(d) verwende die Potenzreihe
38 // atanh(x)/x = sum(j=0..inf,(x^2)^j/(2j+1)):
39 // a:=x^2, b:=1, i:=1, sum:=0,
40 // while (/= sum (setq sum (+ sum (/ b i)))) do i:=i+2, b:=b*a.
42 // Sonst setze y := x/(1+sqrt(1-x^2)), berechne rekursiv z:=atanh(y)
43 // und liefere 2*z = (scale-float z 1).
44 // Diese Rekursion wird entrekursiviert. Statt k mal hintereinander
45 // x := x/(1+sqrt(1-x^2)) zu bilden, arbeitet man lieber mit den Kehrwerten,
46 // setzt also x := 1/|x|, dann k mal x := x+sqrt(x^2-1), dann x := +- 1/x.
47 // Aufwand: asymptotisch d^2.5 .
49 const cl_LF atanhx (const cl_LF& x)
53 var uintC actuallen = TheLfloat(x)->len;
54 var uintC d = float_digits(x);
55 var sintE e = float_exponent_inline(x);
56 if (e <= (sintC)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2)
57 return x; // ja -> x als Ergebnis
58 if (actuallen >= 34) {
60 var cl_LF xx = extend(x,TheLfloat(x)->len+ceiling((uintE)(-e),intDsize));
61 return cl_float(scale_float(ln((1+xx)/(1-xx)),-1),x);
63 var uintL k = 0; // Rekursionszähler k:=0
64 // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
65 // angewandt werden. limit_slope = 1.0 ist schlecht (ca. 15% zu
66 // schlecht). Ein guter Wert ist:
67 // für naive1: limit_scope = 0.625 = 5/8,
68 // für naive2: limit_scope = 0.4 = 13/32.
69 var uintL sqrt_d = floor(isqrtC(d)*13,32); // limit_slope*floor(sqrt(d))
71 if (e >= (sintL)(-sqrt_d)) {
72 // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
73 var sintL e_limit = 1+sqrt_d; // 1+limit_slope*floor(sqrt(d))
74 xx = recip(abs(xx)); // 1/|x|
76 // nächstes x nach der Formel x := x+sqrt(x^2 - 1) berechnen:
77 xx = sqrt(square(xx) + cl_float(-1,xx)) + xx;
79 } until (float_exponent_inline(xx) > e_limit);
80 // Schleifenende mit Exponent(x) > 1+limit_slope*floor(sqrt(d)),
81 // also x >= 2^(1+limit_slope*floor(sqrt(d))),
82 // also 1/x <= 2^(-1-limit_slope*floor(sqrt(d))).
83 // Nun kann die Potenzreihe auf 1/x angewandt werden.
86 xx = - xx; // Vorzeichen wieder rein
88 // Potenzreihe anwenden:
90 var cl_LF a = square(xx); // a = x^2
91 var cl_LF b = cl_float(1,xx); // b := (float 1 x)
92 var cl_LF sum = cl_float(0,xx); // sum := (float 0 x)
95 // floating-point representation
97 var cl_LF new_sum = sum + b/(cl_I)i; // (+ sum (/ b i))
98 if (new_sum == sum) // = sum ?
99 break; // ja -> Potenzreihe abbrechen
106 // floating-point representation with smooth precision reduction
107 var cl_LF eps = scale_float(b,-(sintC)d-10);
109 var cl_LF new_sum = sum + LF_to_LF(b/(cl_I)i,actuallen); // (+ sum (/ b i))
110 if (new_sum == sum) // = sum ?
111 break; // ja -> Potenzreihe abbrechen
113 b = cl_LF_shortenwith(b,eps);
118 var cl_LF erg = sum*xx; // sum*x als Ergebnis
119 return scale_float(erg,k); // wegen Rekursion noch mal 2^k
121 // Bit complexity (N = length(x)): O(log(N)^2*M(N)).
123 const cl_F atanhx (const cl_F& x)
126 DeclareType(cl_LF,x);
131 var uintC d = float_digits(x);
132 var sintE e = float_exponent(x);
133 if (e <= (sintC)(-d)>>1) // e <= -d/2 <==> e <= -ceiling(d/2)
134 return x; // ja -> x als Ergebnis
135 var uintL k = 0; // Rekursionszähler k:=0
136 // Bei e <= -1-limit_slope*floor(sqrt(d)) kann die Potenzreihe
137 // angewandt werden. limit_slope = 1.0 ist schlecht (ca. 15% zu
138 // schlecht). Ein guter Wert ist limit_scope = 0.625 = 5/8.
139 var uintL sqrt_d = floor(isqrtC(d)*5,8); // limit_slope*floor(sqrt(d))
141 if (e >= (sintL)(-sqrt_d)) {
142 // e > -1-limit_slope*floor(sqrt(d)) -> muß |x| verkleinern.
143 var sintL e_limit = 1+sqrt_d; // 1+limit_slope*floor(sqrt(d))
144 xx = recip(abs(xx)); // 1/|x|
146 // nächstes x nach der Formel x := x+sqrt(x^2 - 1) berechnen:
147 xx = sqrt(square(xx) + cl_float(-1,xx)) + xx;
149 } until (float_exponent(xx) > e_limit);
150 // Schleifenende mit Exponent(x) > 1+limit_slope*floor(sqrt(d)),
151 // also x >= 2^(1+limit_slope*floor(sqrt(d))),
152 // also 1/x <= 2^(-1-limit_slope*floor(sqrt(d))).
153 // Nun kann die Potenzreihe auf 1/x angewandt werden.
156 xx = - xx; // Vorzeichen wieder rein
158 // Potenzreihe anwenden:
160 var cl_F a = square(xx); // a = x^2
161 var cl_F b = cl_float(1,xx); // b := (float 1 x)
162 var cl_F sum = cl_float(0,xx); // sum := (float 0 x)
164 var cl_F new_sum = sum + b / (cl_I)i; // (+ sum (/ b i))
165 if (new_sum == sum) // = sum ?
166 break; // ja -> Potenzreihe abbrechen
171 var cl_F erg = sum*xx; // sum*x als Ergebnis
172 return scale_float(erg,k); // wegen Rekursion noch mal 2^k
174 // Bit complexity (N = length(x)): O(log(N)^2*M(N)).
176 // Timings of the above algorithms, on an i486 33 MHz, running Linux,
177 // applied to x = sqrt(2)-1 = 0.414...
178 // N naive1 naive2 use ln
179 // 10 0.013 0.013 0.015
180 // 25 0.064 0.050 0.049
181 // 50 0.25 0.018 0.17
182 // 100 1.07 0.75 0.64
186 // ==> using ln faster for N >= 34.