// 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 <stdlib.h>
#include <string.h>
#include "cln/timing.h"
#undef floor
#include <math.h>
#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*(sintL)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 uintL 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 uintL 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*(sintL)actuallen);
}
const cl_LF atan_recip_1c (cl_I m, uintC len)
{
var uintC actuallen = len + 1;
var cl_I m2 = m*m;
var sintL N = (sintL)(0.69314718*intDsize/2*actuallen/log(double_approx(m))) + 1;
var cl_I num = 0, den = 1; // "lazy rational number"
for (sintL 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 uintL N = (uintL)(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 uintL 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*(sintL)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 uintL 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 uintL 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*(sintL)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 uintL N = (uintL)(0.69314718*intDsize*actuallen/log(double_approx(m2))) + 1;
var cl_I num = 0, den = 1; // "lazy rational number"
for (uintL 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 uintL N = (uintL)(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 uintL 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];
uintL len = atoi(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
