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

#include "cln/integer.h"
#include "cln/rational.h"
#include "cln/real.h"
#include "cln/complex.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


// Method 1: atanh(1/m) = sum(n=0..infty, 1/(2n+1) * 1/m^(2n+1))
// Method 2: atanh(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 4: atanh(x) = 1/2 ln((1+x)/(1-x)),
// using the general built-in algorithm      [FAST]


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

const cl_LF atanh_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);
		fsum = fsum + LF_to_LF(fterm/(2*n+1),actuallen);
	}
	return shorten(fsum,len);
}

const cl_LF atanh_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);
		fsum = fsum + floor1(fterm,2*n+1);
	}
	return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintL)actuallen);
}

const cl_LF atanh_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/(2n+1):
		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 atanh_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 ((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: atanh(1/m) = sum(n=0..infty, (-4)^n*n!^2/(2n+1)! * m/(m^2-1)^(n+1))

const cl_LF atanh_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);
		if ((n % 2) == 0)
			fsum = fsum + LF_to_LF(fterm,actuallen);
		else
			fsum = fsum - LF_to_LF(fterm,actuallen);
	}
	return shorten(fsum,len);
}

const cl_LF atanh_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);
		if ((n % 2) == 0)
			fsum = fsum + fterm;
		else
			fsum = fsum - fterm;
	}
	return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintL)actuallen);
}

const cl_LF atanh_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 atanh_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 (-(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 = atanh_recip_1a(m,len); }
	}
	cout << p << endl;
	{ CL_TIMING;
	  for (int rep = repetitions; rep > 0; rep--)
	    { p = atanh_recip_1b(m,len); }
	}
	cout << p << endl;
	{ CL_TIMING;
	  for (int rep = repetitions; rep > 0; rep--)
	    { p = atanh_recip_1c(m,len); }
	}
	cout << p << endl;
	{ CL_TIMING;
	  for (int rep = repetitions; rep > 0; rep--)
	    { p = atanh_recip_1d(m,len); }
	}
	cout << p << endl;
	// Method 2.
	{ CL_TIMING;
	  for (int rep = repetitions; rep > 0; rep--)
	    { p = atanh_recip_2a(m,len); }
	}
	cout << p << endl;
	{ CL_TIMING;
	  for (int rep = repetitions; rep > 0; rep--)
	    { p = atanh_recip_2b(m,len); }
	}
	cout << p << endl;
	{ CL_TIMING;
	  for (int rep = repetitions; rep > 0; rep--)
	    { p = atanh_recip_2c(m,len); }
	}
	cout << p << endl;
	{ CL_TIMING;
	  for (int rep = repetitions; rep > 0; rep--)
	    { p = atanh_recip_2d(m,len); }
	}
	cout << p << endl;
	// Method 3.
	{ CL_TIMING;
	  for (int rep = repetitions; rep > 0; rep--)
	    { p = The(cl_LF)(atanh(cl_RA_to_LF(1/(cl_RA)m,len))); }
	}
	cout << p << endl;
	// Method 4.
	{ CL_TIMING;
	  for (int rep = repetitions; rep > 0; rep--)
	    { p = The(cl_LF)(scale_float(ln(cl_RA_to_LF((cl_RA)(m+1)/(cl_RA)(m-1),len)),-1)); }
	}
	cout << p << endl;
}


// Timings of the above algorithms, on an i486 33 MHz, running Linux.
// m = 3 -> 1/2 ln(2)
//   N    1a     1b     1c     1d     2a     2b     2c     2d      3
//   10   0.021  0.014  0.019  0.012  0.029  0.015  0.023  0.015   0.015
//   25   0.060  0.041  0.073  0.041  0.082  0.046  0.086  0.051   0.066
//   50   0.164  0.110  0.258  0.120  0.203  0.124  0.295  0.142
//  100   0.49   0.35   1.05   0.37   0.60   0.35   1.19   0.42
//  250   2.5    1.9    7.2    1.7    2.9    1.9    8.0    1.8
//  500  10.1    7.2   33.4    5.5   10.7    7.3   36.5    5.9
// 1000  38     30    145     16.1   39     29    158     16.8
// 2500 231    188    976     53    237    186   1081     58
// asymp. N^2    N^2    N^2    FAST   N^2    N^2    N^2    FAST
//
// m = 9 -> 1/2 ln(5/4)
//   N    1a     1b     1c     1d     2a     2b     2c     2d     3      4
//   10   0.0106 0.0072 0.0084 0.0061 0.0139 0.0073 0.0098 0.0073 0.0140 0.0211
//   25   0.031  0.021  0.029  0.019  0.039  0.022  0.031  0.022  0.063  0.081
//   50   0.083  0.057  0.091  0.056  0.098  0.058  0.098  0.060  0.232  0.212
//  100   0.25   0.17   0.32   0.16   0.28   0.17   0.34   0.17   0.60   0.59
//  250   1.28   0.94   2.11   0.77   1.40   0.91   2.18   0.76   2.76   2.76
//  500   5.1    3.6    9.4    2.5    5.2    3.4    9.3    2.4   10.4    9.7
// 1000  19.1   14.7   42      7.8   18.5   13.6   42      7.4   31     30
// 2500 116     93    279     29.6  113     86    278     30.0  129    125
// asymp. N^2    N^2    N^2    FAST   N^2    N^2    N^2    FAST   FAST   FAST