cln/examples/atan_recip.cc

// Computation of arctan(1/m) (m integer) to high precision.

#include "cln/integer.h"
#include "cln/rational.h"
#include "cln/real.h"
#include "cln/lfloat.h"
#include "cl_LF.h"
#include "cl_LF_tran.h"
#include "cl_alloca.h"
#include <cstdlib>
#include <cstring>
#include "cln/timing.h"

#undef floor
#include <cmath>
#define floor cln_floor

using namespace cln;

// Method 1: atan(1/m) = sum(n=0..infty, (-1)^n/(2n+1) * 1/m^(2n+1))
// Method 2: atan(1/m) = sum(n=0..infty, 4^n*n!^2/(2n+1)! * m/(m^2+1)^(n+1))
// a. Using long floats.                     [N^2]
// b. Simulating long floats using integers. [N^2]
// c. Using integers, no binary splitting.   [N^2]
// d. Using integers, with binary splitting. [FAST]
// Method 3: general built-in algorithm.     [FAST]


// Method 1: atan(1/m) = sum(n=0..infty, (-1)^n/(2n+1) * 1/m^(2n+1))

const cl_LF atan_recip_1a (cl_I m, uintC len)
{
        var uintC actuallen = len + 1;
        var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintC)actuallen);
        var cl_I m2 = m*m;
        var cl_LF fterm = cl_I_to_LF(1,actuallen)/m;
        var cl_LF fsum = fterm;
        for (var uintC n = 1; fterm >= eps; n++) {
                fterm = fterm/m2;
                fterm = cl_LF_shortenwith(fterm,eps);
                if ((n % 2) == 0)
                        fsum = fsum + LF_to_LF(fterm/(2*n+1),actuallen);
                else
                        fsum = fsum - LF_to_LF(fterm/(2*n+1),actuallen);
        }
        return shorten(fsum,len);
}

const cl_LF atan_recip_1b (cl_I m, uintC len)
{
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m;
        var cl_I fterm = floor1((cl_I)1 << (intDsize*actuallen), m);
        var cl_I fsum = fterm;
        for (var uintC n = 1; fterm > 0; n++) {
                fterm = floor1(fterm,m2);
                if ((n % 2) == 0)
                        fsum = fsum + floor1(fterm,2*n+1);
                else
                        fsum = fsum - floor1(fterm,2*n+1);
        }
        return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintC)actuallen);
}

const cl_LF atan_recip_1c (cl_I m, uintC len)
{
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m;
        var sintC N = (sintC)(0.69314718*intDsize/2*actuallen/log(double_approx(m))) + 1;
        var cl_I num = 0, den = 1; // "lazy rational number"
        for (sintC n = N-1; n>=0; n--) {
                // Multiply sum with 1/m^2:
                den = den * m2;
                // Add (-1)^n/(2n+1):
                if ((n % 2) == 0)
                        num = num*(2*n+1) + den;
                else
                        num = num*(2*n+1) - den;
                den = den*(2*n+1);
        }
        den = den*m;
        var cl_LF result = cl_I_to_LF(num,actuallen)/cl_I_to_LF(den,actuallen);
        return shorten(result,len);
}

const cl_LF atan_recip_1d (cl_I m, uintC len)
{
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m;
        var uintC N = (uintC)(0.69314718*intDsize/2*actuallen/log(double_approx(m))) + 1;
        CL_ALLOCA_STACK;
        var cl_I* bv = (cl_I*) cl_alloca(N*sizeof(cl_I));
        var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
        var uintC n;
        for (n = 0; n < N; n++) {
                new (&bv[n]) cl_I ((n % 2) == 0 ? (cl_I)(2*n+1) : -(cl_I)(2*n+1));
                new (&qv[n]) cl_I (n==0 ? m : m2);
        }
        var cl_rational_series series;
        series.av = NULL; series.bv = bv;
        series.pv = NULL; series.qv = qv; series.qsv = NULL;
        var cl_LF result = eval_rational_series(N,series,actuallen);
        for (n = 0; n < N; n++) {
                bv[n].~cl_I();
                qv[n].~cl_I();
        }
        return shorten(result,len);
}


// Method 2: atan(1/m) = sum(n=0..infty, 4^n*n!^2/(2n+1)! * m/(m^2+1)^(n+1))

const cl_LF atan_recip_2a (cl_I m, uintC len)
{
        var uintC actuallen = len + 1;
        var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintC)actuallen);
        var cl_I m2 = m*m+1;
        var cl_LF fterm = cl_I_to_LF(m,actuallen)/m2;
        var cl_LF fsum = fterm;
        for (var uintC n = 1; fterm >= eps; n++) {
                fterm = The(cl_LF)((2*n)*fterm)/((2*n+1)*m2);
                fterm = cl_LF_shortenwith(fterm,eps);
                fsum = fsum + LF_to_LF(fterm,actuallen);
        }
        return shorten(fsum,len);
}

const cl_LF atan_recip_2b (cl_I m, uintC len)
{
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m+1;
        var cl_I fterm = floor1((cl_I)m << (intDsize*actuallen), m2);
        var cl_I fsum = fterm;
        for (var uintC n = 1; fterm > 0; n++) {
                fterm = floor1((2*n)*fterm,(2*n+1)*m2);
                fsum = fsum + fterm;
        }
        return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintC)actuallen);
}

const cl_LF atan_recip_2c (cl_I m, uintC len)
{
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m+1;
        var uintC N = (uintC)(0.69314718*intDsize*actuallen/log(double_approx(m2))) + 1;
        var cl_I num = 0, den = 1; // "lazy rational number"
        for (uintC n = N; n>0; n--) {
                // Multiply sum with (2n)/(2n+1)(m^2+1):
                num = num * (2*n);
                den = den * ((2*n+1)*m2);
                // Add 1:
                num = num + den;
        }
        num = num*m;
        den = den*m2;
        var cl_LF result = cl_I_to_LF(num,actuallen)/cl_I_to_LF(den,actuallen);
        return shorten(result,len);
}

const cl_LF atan_recip_2d (cl_I m, uintC len)
{
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m+1;
        var uintC N = (uintC)(0.69314718*intDsize*actuallen/log(double_approx(m2))) + 1;
        CL_ALLOCA_STACK;
        var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
        var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
        var uintC n;
        new (&pv[0]) cl_I (m);
        new (&qv[0]) cl_I (m2);
        for (n = 1; n < N; n++) {
                new (&pv[n]) cl_I (2*n);
                new (&qv[n]) cl_I ((2*n+1)*m2);
        }
        var cl_rational_series series;
        series.av = NULL; series.bv = NULL;
        series.pv = pv; series.qv = qv; series.qsv = NULL;
        var cl_LF result = eval_rational_series(N,series,actuallen);
        for (n = 0; n < N; n++) {
                pv[n].~cl_I();
                qv[n].~cl_I();
        }
        return shorten(result,len);
}


// Main program: Compute and display the timings.

int main (int argc, char * argv[])
{
        int repetitions = 1;
        if ((argc >= 3) && !strcmp(argv[1],"-r")) {
                repetitions = atoi(argv[2]);
                argc -= 2; argv += 2;
        }
        if (argc < 2)
                exit(1);
        cl_I m = (cl_I)argv[1];
        uintC len = atol(argv[2]);
        cl_LF p;
        ln(cl_I_to_LF(1000,len+10)); // fill cache
        // Method 1.
        { CL_TIMING;
          for (int rep = repetitions; rep > 0; rep--)
            { p = atan_recip_1a(m,len); }
        }
        cout << p << endl;
        { CL_TIMING;
          for (int rep = repetitions; rep > 0; rep--)
            { p = atan_recip_1b(m,len); }
        }
        cout << p << endl;
        { CL_TIMING;
          for (int rep = repetitions; rep > 0; rep--)
            { p = atan_recip_1c(m,len); }
        }
        cout << p << endl;
        { CL_TIMING;
          for (int rep = repetitions; rep > 0; rep--)
            { p = atan_recip_1d(m,len); }
        }
        cout << p << endl;
        // Method 2.
        { CL_TIMING;
          for (int rep = repetitions; rep > 0; rep--)
            { p = atan_recip_2a(m,len); }
        }
        cout << p << endl;
        { CL_TIMING;
          for (int rep = repetitions; rep > 0; rep--)
            { p = atan_recip_2b(m,len); }
        }
        cout << p << endl;
        { CL_TIMING;
          for (int rep = repetitions; rep > 0; rep--)
            { p = atan_recip_2c(m,len); }
        }
        cout << p << endl;
        { CL_TIMING;
          for (int rep = repetitions; rep > 0; rep--)
            { p = atan_recip_2d(m,len); }
        }
        cout << p << endl;
        // Method 3.
        { CL_TIMING;
          for (int rep = repetitions; rep > 0; rep--)
            { p = The(cl_LF)(atan(cl_RA_to_LF(1/(cl_RA)m,len))); }
        }
        cout << p << endl;
}


// Timings of the above algorithms, on an i486 33 MHz, running Linux.
// m = 390112. (For Jörg Arndt's formula (4.15).)
//    N      1a     1b     1c     1d      2a     2b     2c     2d      3
//    10     0.0027 0.0018 0.0019 0.0019  0.0032 0.0022 0.0019 0.0019  0.0042
//    25     0.0085 0.0061 0.0058 0.0061  0.0095 0.0069 0.0056 0.0061  0.028
//    50     0.024  0.018  0.017  0.017   0.026  0.020  0.016  0.017   0.149
//   100     0.075  0.061  0.057  0.054   0.079  0.065  0.052  0.052   0.71
//   250     0.41   0.33   0.32   0.26    0.42   0.36   0.28   0.24    3.66
//   500     1.57   1.31   1.22   0.88    1.57   1.36   1.10   0.83   13.7
//  1000     6.08   5.14   4.56   2.76    6.12   5.36   4.06   2.58   45.5
//  2500    36.5   32.2   25.8   10.2    38.4   33.6   22.2    9.1   191
//  5000
// 10000
// asymp.    N^2    N^2    N^2    FAST    N^2    N^2    N^2    FAST    FAST
//
// m = 319. (For Jörg Arndt's formula (4.7).)
//    N      1a     1b     1c     1d      2a     2b     2c     2d      3
//  1000     6.06   4.40   9.17   3.82    5.29   3.90   7.50   3.53   50.3
//
// m = 18. (For Jörg Arndt's formula (4.4).)
//    N      1a     1b     1c     1d      2a     2b     2c     2d      3
//  1000    11.8    9.0   22.3    6.0    10.2    7.7   17.1    5.7    54.3
1	kreckel	1.1	// Computation of arctan(1/m) (m integer) to high precision.
2
3	kreckel	1.2	#include "cln/integer.h"
4			#include "cln/rational.h"
5			#include "cln/real.h"
6			#include "cln/lfloat.h"
7	kreckel	1.1	#include "cl_LF.h"
8			#include "cl_LF_tran.h"
9			#include "cl_alloca.h"
10	kreckel	1.3	#include <cstdlib>
11			#include <cstring>
12	kreckel	1.2	#include "cln/timing.h"
13	kreckel	1.1
14			#undef floor
15	kreckel	1.3	#include <cmath>
16	kreckel	1.1	#define floor cln_floor
17
18	kreckel	1.2	using namespace cln;
19	kreckel	1.1
20			// Method 1: atan(1/m) = sum(n=0..infty, (-1)^n/(2n+1) * 1/m^(2n+1))
21			// Method 2: atan(1/m) = sum(n=0..infty, 4^nn!^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
28
29			// Method 1: atan(1/m) = sum(n=0..infty, (-1)^n/(2n+1) * 1/m^(2n+1))
30
31			const cl_LF atan_recip_1a (cl_I m, uintC len)
32			{
33			var uintC actuallen = len + 1;
34	kreckel	1.4	var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintC)actuallen);
35	kreckel	1.1	var cl_I m2 = m*m;
36			var cl_LF fterm = cl_I_to_LF(1,actuallen)/m;
37			var cl_LF fsum = fterm;
38	kreckel	1.4	for (var uintC n = 1; fterm >= eps; n++) {
39	kreckel	1.1	fterm = fterm/m2;
40			fterm = cl_LF_shortenwith(fterm,eps);
41			if ((n % 2) == 0)
42			fsum = fsum + LF_to_LF(fterm/(2*n+1),actuallen);
43			else
44			fsum = fsum - LF_to_LF(fterm/(2*n+1),actuallen);
45			}
46			return shorten(fsum,len);
47			}
48
49			const cl_LF atan_recip_1b (cl_I m, uintC len)
50			{
51			var uintC actuallen = len + 1;
52			var cl_I m2 = m*m;
53			var cl_I fterm = floor1((cl_I)1 << (intDsize*actuallen), m);
54			var cl_I fsum = fterm;
55	kreckel	1.4	for (var uintC n = 1; fterm > 0; n++) {
56	kreckel	1.1	fterm = floor1(fterm,m2);
57			if ((n % 2) == 0)
58			fsum = fsum + floor1(fterm,2*n+1);
59			else
60			fsum = fsum - floor1(fterm,2*n+1);
61			}
62	kreckel	1.4	return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintC)actuallen);
63	kreckel	1.1	}
64
65			const cl_LF atan_recip_1c (cl_I m, uintC len)
66			{
67			var uintC actuallen = len + 1;
68			var cl_I m2 = m*m;
69	kreckel	1.4	var sintC N = (sintC)(0.69314718intDsize/2actuallen/log(double_approx(m))) + 1;
70	kreckel	1.1	var cl_I num = 0, den = 1; // "lazy rational number"
71	kreckel	1.4	for (sintC n = N-1; n>=0; n--) {
72	kreckel	1.1	// Multiply sum with 1/m^2:
73			den = den * m2;
74			// Add (-1)^n/(2n+1):
75			if ((n % 2) == 0)
76			num = num(2n+1) + den;
77			else
78			num = num(2n+1) - den;
79			den = den(2n+1);
80			}
81			den = den*m;
82			var cl_LF result = cl_I_to_LF(num,actuallen)/cl_I_to_LF(den,actuallen);
83			return shorten(result,len);
84			}
85
86			const cl_LF atan_recip_1d (cl_I m, uintC len)
87			{
88			var uintC actuallen = len + 1;
89			var cl_I m2 = m*m;
90	kreckel	1.4	var uintC N = (uintC)(0.69314718intDsize/2actuallen/log(double_approx(m))) + 1;
91	kreckel	1.1	CL_ALLOCA_STACK;
92			var cl_I* bv = (cl_I) cl_alloca(Nsizeof(cl_I));
93			var cl_I* qv = (cl_I) cl_alloca(Nsizeof(cl_I));
94	kreckel	1.4	var uintC n;
95	kreckel	1.1	for (n = 0; n < N; n++) {
96			new (&bv[n]) cl_I ((n % 2) == 0 ? (cl_I)(2n+1) : -(cl_I)(2n+1));
97			new (&qv[n]) cl_I (n==0 ? m : m2);
98			}
99			var cl_rational_series series;
100			series.av = NULL; series.bv = bv;
101			series.pv = NULL; series.qv = qv; series.qsv = NULL;
102			var cl_LF result = eval_rational_series(N,series,actuallen);
103			for (n = 0; n < N; n++) {
104			bv[n].~cl_I();
105			qv[n].~cl_I();
106			}
107			return shorten(result,len);
108			}
109
110
111			// Method 2: atan(1/m) = sum(n=0..infty, 4^nn!^2/(2n+1)! m/(m^2+1)^(n+1))
112
113			const cl_LF atan_recip_2a (cl_I m, uintC len)
114			{
115			var uintC actuallen = len + 1;
116	kreckel	1.4	var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintC)actuallen);
117	kreckel	1.1	var cl_I m2 = m*m+1;
118			var cl_LF fterm = cl_I_to_LF(m,actuallen)/m2;
119			var cl_LF fsum = fterm;
120	kreckel	1.4	for (var uintC n = 1; fterm >= eps; n++) {
121	kreckel	1.1	fterm = The(cl_LF)((2n)fterm)/((2n+1)m2);
122			fterm = cl_LF_shortenwith(fterm,eps);
123			fsum = fsum + LF_to_LF(fterm,actuallen);
124			}
125			return shorten(fsum,len);
126			}
127
128			const cl_LF atan_recip_2b (cl_I m, uintC len)
129			{
130			var uintC actuallen = len + 1;
131			var cl_I m2 = m*m+1;
132			var cl_I fterm = floor1((cl_I)m << (intDsize*actuallen), m2);
133			var cl_I fsum = fterm;
134	kreckel	1.4	for (var uintC n = 1; fterm > 0; n++) {
135	kreckel	1.1	fterm = floor1((2n)fterm,(2n+1)m2);
136			fsum = fsum + fterm;
137			}
138	kreckel	1.4	return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintC)actuallen);
139	kreckel	1.1	}
140
141			const cl_LF atan_recip_2c (cl_I m, uintC len)
142			{
143			var uintC actuallen = len + 1;
144			var cl_I m2 = m*m+1;
145	kreckel	1.4	var uintC N = (uintC)(0.69314718intDsizeactuallen/log(double_approx(m2))) + 1;
146	kreckel	1.1	var cl_I num = 0, den = 1; // "lazy rational number"
147	kreckel	1.4	for (uintC n = N; n>0; n--) {
148	kreckel	1.1	// Multiply sum with (2n)/(2n+1)(m^2+1):
149			num = num * (2*n);
150			den = den * ((2n+1)m2);
151			// Add 1:
152			num = num + den;
153			}
154			num = num*m;
155			den = den*m2;
156			var cl_LF result = cl_I_to_LF(num,actuallen)/cl_I_to_LF(den,actuallen);
157			return shorten(result,len);
158			}
159
160			const cl_LF atan_recip_2d (cl_I m, uintC len)
161			{
162			var uintC actuallen = len + 1;
163			var cl_I m2 = m*m+1;
164	kreckel	1.4	var uintC N = (uintC)(0.69314718intDsizeactuallen/log(double_approx(m2))) + 1;
165	kreckel	1.1	CL_ALLOCA_STACK;
166			var cl_I* pv = (cl_I) cl_alloca(Nsizeof(cl_I));
167			var cl_I* qv = (cl_I) cl_alloca(Nsizeof(cl_I));
168	kreckel	1.4	var uintC n;
169	kreckel	1.1	new (&pv[0]) cl_I (m);
170			new (&qv[0]) cl_I (m2);
171			for (n = 1; n < N; n++) {
172			new (&pv[n]) cl_I (2*n);
173			new (&qv[n]) cl_I ((2n+1)m2);
174			}
175			var cl_rational_series series;
176			series.av = NULL; series.bv = NULL;
177			series.pv = pv; series.qv = qv; series.qsv = NULL;
178			var cl_LF result = eval_rational_series(N,series,actuallen);
179			for (n = 0; n < N; n++) {
180			pv[n].~cl_I();
181			qv[n].~cl_I();
182			}
183			return shorten(result,len);
184			}
185
186
187			// Main program: Compute and display the timings.
188
189			int main (int argc, char * argv[])
190			{
191			int repetitions = 1;
192			if ((argc >= 3) && !strcmp(argv[1],"-r")) {
193			repetitions = atoi(argv[2]);
194			argc -= 2; argv += 2;
195			}
196			if (argc < 2)
197			exit(1);
198			cl_I m = (cl_I)argv[1];
199	kreckel	1.4	uintC len = atol(argv[2]);
200	kreckel	1.1	cl_LF p;
201			ln(cl_I_to_LF(1000,len+10)); // fill cache
202			// Method 1.
203			{ CL_TIMING;
204			for (int rep = repetitions; rep > 0; rep--)
205			{ p = atan_recip_1a(m,len); }
206			}
207			cout << p << endl;
208			{ CL_TIMING;
209			for (int rep = repetitions; rep > 0; rep--)
210			{ p = atan_recip_1b(m,len); }
211			}
212			cout << p << endl;
213			{ CL_TIMING;
214			for (int rep = repetitions; rep > 0; rep--)
215			{ p = atan_recip_1c(m,len); }
216			}
217			cout << p << endl;
218			{ CL_TIMING;
219			for (int rep = repetitions; rep > 0; rep--)
220			{ p = atan_recip_1d(m,len); }
221			}
222			cout << p << endl;
223			// Method 2.
224			{ CL_TIMING;
225			for (int rep = repetitions; rep > 0; rep--)
226			{ p = atan_recip_2a(m,len); }
227			}
228			cout << p << endl;
229			{ CL_TIMING;
230			for (int rep = repetitions; rep > 0; rep--)
231			{ p = atan_recip_2b(m,len); }
232			}
233			cout << p << endl;
234			{ CL_TIMING;
235			for (int rep = repetitions; rep > 0; rep--)
236			{ p = atan_recip_2c(m,len); }
237			}
238			cout << p << endl;
239			{ CL_TIMING;
240			for (int rep = repetitions; rep > 0; rep--)
241			{ p = atan_recip_2d(m,len); }
242			}
243			cout << p << endl;
244			// Method 3.
245			{ CL_TIMING;
246			for (int rep = repetitions; rep > 0; rep--)
247			{ p = The(cl_LF)(atan(cl_RA_to_LF(1/(cl_RA)m,len))); }
248			}
249			cout << p << endl;
250			}
251
252
253			// Timings of the above algorithms, on an i486 33 MHz, running Linux.
254	kreckel	1.5	// m = 390112. (For Jörg Arndt's formula (4.15).)
255	kreckel	1.1	// N 1a 1b 1c 1d 2a 2b 2c 2d 3
256			// 10 0.0027 0.0018 0.0019 0.0019 0.0032 0.0022 0.0019 0.0019 0.0042
257			// 25 0.0085 0.0061 0.0058 0.0061 0.0095 0.0069 0.0056 0.0061 0.028
258			// 50 0.024 0.018 0.017 0.017 0.026 0.020 0.016 0.017 0.149
259			// 100 0.075 0.061 0.057 0.054 0.079 0.065 0.052 0.052 0.71
260			// 250 0.41 0.33 0.32 0.26 0.42 0.36 0.28 0.24 3.66
261			// 500 1.57 1.31 1.22 0.88 1.57 1.36 1.10 0.83 13.7
262			// 1000 6.08 5.14 4.56 2.76 6.12 5.36 4.06 2.58 45.5
263			// 2500 36.5 32.2 25.8 10.2 38.4 33.6 22.2 9.1 191
264			// 5000
265			// 10000
266			// asymp. N^2 N^2 N^2 FAST N^2 N^2 N^2 FAST FAST
267			//
268	kreckel	1.5	// m = 319. (For Jörg Arndt's formula (4.7).)
269	kreckel	1.1	// N 1a 1b 1c 1d 2a 2b 2c 2d 3
270			// 1000 6.06 4.40 9.17 3.82 5.29 3.90 7.50 3.53 50.3
271			//
272	kreckel	1.5	// m = 18. (For Jörg Arndt's formula (4.4).)
273	kreckel	1.1	// N 1a 1b 1c 1d 2a 2b 2c 2d 3
274			// 1000 11.8 9.0 22.3 6.0 10.2 7.7 17.1 5.7 54.3
Christian Bauer">Christian Bauer	ViewVC Help
Powered by ViewVC 1.1.15