// Computation of arctan(1/m) (m integer) to high precision.
-#include "cl_integer.h"
-#include "cl_rational.h"
-#include "cl_real.h"
-#include "cl_lfloat.h"
+#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 "cl_timing.h"
+#include <cstdlib>
+#include <cstring>
+#include "cln/timing.h"
#undef floor
-#include <math.h>
+#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))
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_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 uintL n = 1; fterm >= eps; n++) {
+ for (var uintC n = 1; fterm >= eps; n++) {
fterm = fterm/m2;
fterm = cl_LF_shortenwith(fterm,eps);
if ((n % 2) == 0)
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++) {
+ 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*(sintL)actuallen);
+ 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 sintL N = (sintL)(0.69314718*intDsize/2*actuallen/log(cl_double_approx(m))) + 1;
+ 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 (sintL n = N-1; n>=0; n--) {
+ for (sintC n = N-1; n>=0; n--) {
// Multiply sum with 1/m^2:
den = den * m2;
// Add (-1)^n/(2n+1):
{
var uintC actuallen = len + 1;
var cl_I m2 = m*m;
- var uintL N = (uintL)(0.69314718*intDsize/2*actuallen/log(cl_double_approx(m))) + 1;
+ 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 uintL n;
+ 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);
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_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 uintL n = 1; fterm >= eps; n++) {
+ 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);
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++) {
+ 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*(sintL)actuallen);
+ 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 uintL N = (uintL)(0.69314718*intDsize*actuallen/log(cl_double_approx(m2))) + 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 (uintL n = N; n>0; n--) {
+ 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);
{
var uintC actuallen = len + 1;
var cl_I m2 = m*m+1;
- var uintL N = (uintL)(0.69314718*intDsize*actuallen/log(cl_double_approx(m2))) + 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 uintL n;
+ var uintC n;
new (&pv[0]) cl_I (m);
new (&qv[0]) cl_I (m2);
for (n = 1; n < N; n++) {
if (argc < 2)
exit(1);
cl_I m = (cl_I)argv[1];
- uintL len = atoi(argv[2]);
+ uintC len = atol(argv[2]);
cl_LF p;
ln(cl_I_to_LF(1000,len+10)); // fill cache
// Method 1.
// Timings of the above algorithms, on an i486 33 MHz, running Linux.
-// m = 390112. (For Jörg Arndt's formula (4.15).)
+// 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
// 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).)
+// 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).)
+// 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