1 // Computation of artanh(1/m) (m integer) to high precision.
3 #include "cl_integer.h"
4 #include "cl_rational.h"
6 #include "cl_complex.h"
9 #include "cl_LF_tran.h"
10 #include "cl_alloca.h"
13 #include "cl_timing.h"
17 #define floor cln_floor
20 // Method 1: atanh(1/m) = sum(n=0..infty, 1/(2n+1) * 1/m^(2n+1))
21 // Method 2: atanh(1/m) = sum(n=0..infty, (-4)^n*n!^2/(2n+1)! * m/(m^2-1)^(n+1))
22 // a. Using long floats. [N^2]
23 // b. Simulating long floats using integers. [N^2]
24 // c. Using integers, no binary splitting. [N^2]
25 // d. Using integers, with binary splitting. [FAST]
26 // Method 3: general built-in algorithm. [FAST]
27 // Method 4: atanh(x) = 1/2 ln((1+x)/(1-x)),
28 // using the general built-in algorithm [FAST]
31 // Method 1: atanh(1/m) = sum(n=0..infty, 1/(2n+1) * 1/m^(2n+1))
33 const cl_LF atanh_recip_1a (cl_I m, uintC len)
35 var uintC actuallen = len + 1;
36 var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintL)actuallen);
38 var cl_LF fterm = cl_I_to_LF(1,actuallen)/m;
39 var cl_LF fsum = fterm;
40 for (var uintL n = 1; fterm >= eps; n++) {
42 fterm = cl_LF_shortenwith(fterm,eps);
43 fsum = fsum + LF_to_LF(fterm/(2*n+1),actuallen);
45 return shorten(fsum,len);
48 const cl_LF atanh_recip_1b (cl_I m, uintC len)
50 var uintC actuallen = len + 1;
52 var cl_I fterm = floor1((cl_I)1 << (intDsize*actuallen), m);
53 var cl_I fsum = fterm;
54 for (var uintL n = 1; fterm > 0; n++) {
55 fterm = floor1(fterm,m2);
56 fsum = fsum + floor1(fterm,2*n+1);
58 return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintL)actuallen);
61 const cl_LF atanh_recip_1c (cl_I m, uintC len)
63 var uintC actuallen = len + 1;
65 var sintL N = (sintL)(0.69314718*intDsize/2*actuallen/log(cl_double_approx(m))) + 1;
66 var cl_I num = 0, den = 1; // "lazy rational number"
67 for (sintL n = N-1; n>=0; n--) {
68 // Multiply sum with 1/m^2:
71 num = num*(2*n+1) + den;
75 var cl_LF result = cl_I_to_LF(num,actuallen)/cl_I_to_LF(den,actuallen);
76 return shorten(result,len);
79 const cl_LF atanh_recip_1d (cl_I m, uintC len)
81 var uintC actuallen = len + 1;
83 var uintL N = (uintL)(0.69314718*intDsize/2*actuallen/log(cl_double_approx(m))) + 1;
85 var cl_I* bv = (cl_I*) cl_alloca(N*sizeof(cl_I));
86 var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
88 for (n = 0; n < N; n++) {
89 new (&bv[n]) cl_I ((cl_I)(2*n+1));
90 new (&qv[n]) cl_I (n==0 ? m : m2);
92 var cl_rational_series series;
93 series.av = NULL; series.bv = bv;
94 series.pv = NULL; series.qv = qv; series.qsv = NULL;
95 var cl_LF result = eval_rational_series(N,series,actuallen);
96 for (n = 0; n < N; n++) {
100 return shorten(result,len);
104 // Method 2: atanh(1/m) = sum(n=0..infty, (-4)^n*n!^2/(2n+1)! * m/(m^2-1)^(n+1))
106 const cl_LF atanh_recip_2a (cl_I m, uintC len)
108 var uintC actuallen = len + 1;
109 var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintL)actuallen);
111 var cl_LF fterm = cl_I_to_LF(m,actuallen)/m2;
112 var cl_LF fsum = fterm;
113 for (var uintL n = 1; fterm >= eps; n++) {
114 fterm = The(cl_LF)((2*n)*fterm)/((2*n+1)*m2);
115 fterm = cl_LF_shortenwith(fterm,eps);
117 fsum = fsum + LF_to_LF(fterm,actuallen);
119 fsum = fsum - LF_to_LF(fterm,actuallen);
121 return shorten(fsum,len);
124 const cl_LF atanh_recip_2b (cl_I m, uintC len)
126 var uintC actuallen = len + 1;
128 var cl_I fterm = floor1((cl_I)m << (intDsize*actuallen), m2);
129 var cl_I fsum = fterm;
130 for (var uintL n = 1; fterm > 0; n++) {
131 fterm = floor1((2*n)*fterm,(2*n+1)*m2);
137 return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintL)actuallen);
140 const cl_LF atanh_recip_2c (cl_I m, uintC len)
142 var uintC actuallen = len + 1;
144 var uintL N = (uintL)(0.69314718*intDsize*actuallen/log(cl_double_approx(m2))) + 1;
145 var cl_I num = 0, den = 1; // "lazy rational number"
146 for (uintL n = N; n>0; n--) {
147 // Multiply sum with -(2n)/(2n+1)(m^2+1):
149 den = - den * ((2*n+1)*m2);
155 var cl_LF result = cl_I_to_LF(num,actuallen)/cl_I_to_LF(den,actuallen);
156 return shorten(result,len);
159 const cl_LF atanh_recip_2d (cl_I m, uintC len)
161 var uintC actuallen = len + 1;
163 var uintL N = (uintL)(0.69314718*intDsize*actuallen/log(cl_double_approx(m2))) + 1;
165 var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
166 var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
168 new (&pv[0]) cl_I (m);
169 new (&qv[0]) cl_I (m2);
170 for (n = 1; n < N; n++) {
171 new (&pv[n]) cl_I (-(cl_I)(2*n));
172 new (&qv[n]) cl_I ((2*n+1)*m2);
174 var cl_rational_series series;
175 series.av = NULL; series.bv = NULL;
176 series.pv = pv; series.qv = qv; series.qsv = NULL;
177 var cl_LF result = eval_rational_series(N,series,actuallen);
178 for (n = 0; n < N; n++) {
182 return shorten(result,len);
186 // Main program: Compute and display the timings.
188 int main (int argc, char * argv[])
191 if ((argc >= 3) && !strcmp(argv[1],"-r")) {
192 repetitions = atoi(argv[2]);
193 argc -= 2; argv += 2;
197 cl_I m = (cl_I)argv[1];
198 uintL len = atoi(argv[2]);
200 ln(cl_I_to_LF(1000,len+10)); // fill cache
203 for (int rep = repetitions; rep > 0; rep--)
204 { p = atanh_recip_1a(m,len); }
208 for (int rep = repetitions; rep > 0; rep--)
209 { p = atanh_recip_1b(m,len); }
213 for (int rep = repetitions; rep > 0; rep--)
214 { p = atanh_recip_1c(m,len); }
218 for (int rep = repetitions; rep > 0; rep--)
219 { p = atanh_recip_1d(m,len); }
224 for (int rep = repetitions; rep > 0; rep--)
225 { p = atanh_recip_2a(m,len); }
229 for (int rep = repetitions; rep > 0; rep--)
230 { p = atanh_recip_2b(m,len); }
234 for (int rep = repetitions; rep > 0; rep--)
235 { p = atanh_recip_2c(m,len); }
239 for (int rep = repetitions; rep > 0; rep--)
240 { p = atanh_recip_2d(m,len); }
245 for (int rep = repetitions; rep > 0; rep--)
246 { p = The(cl_LF)(atanh(cl_RA_to_LF(1/(cl_RA)m,len))); }
251 for (int rep = repetitions; rep > 0; rep--)
252 { p = The(cl_LF)(scale_float(ln(cl_RA_to_LF((cl_RA)(m+1)/(cl_RA)(m-1),len)),-1)); }
258 // Timings of the above algorithms, on an i486 33 MHz, running Linux.
259 // m = 3 -> 1/2 ln(2)
260 // N 1a 1b 1c 1d 2a 2b 2c 2d 3
261 // 10 0.021 0.014 0.019 0.012 0.029 0.015 0.023 0.015 0.015
262 // 25 0.060 0.041 0.073 0.041 0.082 0.046 0.086 0.051 0.066
263 // 50 0.164 0.110 0.258 0.120 0.203 0.124 0.295 0.142
264 // 100 0.49 0.35 1.05 0.37 0.60 0.35 1.19 0.42
265 // 250 2.5 1.9 7.2 1.7 2.9 1.9 8.0 1.8
266 // 500 10.1 7.2 33.4 5.5 10.7 7.3 36.5 5.9
267 // 1000 38 30 145 16.1 39 29 158 16.8
268 // 2500 231 188 976 53 237 186 1081 58
269 // asymp. N^2 N^2 N^2 FAST N^2 N^2 N^2 FAST
271 // m = 9 -> 1/2 ln(5/4)
272 // N 1a 1b 1c 1d 2a 2b 2c 2d 3 4
273 // 10 0.0106 0.0072 0.0084 0.0061 0.0139 0.0073 0.0098 0.0073 0.0140 0.0211
274 // 25 0.031 0.021 0.029 0.019 0.039 0.022 0.031 0.022 0.063 0.081
275 // 50 0.083 0.057 0.091 0.056 0.098 0.058 0.098 0.060 0.232 0.212
276 // 100 0.25 0.17 0.32 0.16 0.28 0.17 0.34 0.17 0.60 0.59
277 // 250 1.28 0.94 2.11 0.77 1.40 0.91 2.18 0.76 2.76 2.76
278 // 500 5.1 3.6 9.4 2.5 5.2 3.4 9.3 2.4 10.4 9.7
279 // 1000 19.1 14.7 42 7.8 18.5 13.6 42 7.4 31 30
280 // 2500 116 93 279 29.6 113 86 278 30.0 129 125
281 // asymp. N^2 N^2 N^2 FAST N^2 N^2 N^2 FAST FAST FAST