You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

274 lines
7.9 KiB

// 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;