/*****************************************************/
/*                                                   */
/*    Generator of Correctly Rounded Coefficients    */
/*        of the polynomial which stem from          */
/*          Chebychev Integration Formula            */
/*  based on Classical Error Estimation (CEE) Method */
/*                                                   */
/*      Copyright (c) 2008- Tomonori Kouya           */
/*                                                   */
/* Version 0.0: 2008-06-11 First published Version   */
/*                                                   */
/* This program is free software: you can            */
/* redistribute it and/or modify it under the terms  */
/* of the GNU General Public License as published by */
/* the Free Software Foundation, either version 3 of */
/* the License, or any later version.                */
/*                                                   */
/* This program is distributed in the hope that it   */
/* will be useful, but WITHOUT ANY WARRANTY; without */
/* even the implied warranty of MERCHANTABILITY or   */
/* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU    */
/* General Public License for more details.          */
/*                                                   */
/* You should have received a copy of the GNU        */
/* General Public License along with this program.   */
/* If not, see <http://www.gnu.org/licenses/>.       */
/*                                                   */
/* Please contact me via http://na-inet.jp/          */
/* if you have any comments or requests on this      */
/* program.                                          */
/*                                                   */
/* How to compile:                                   */
/*   $ cc [options] chevcoef.c -lgmp -lm             */
/*                                                   */
/* Compile option:                                   */
/*   -DUSE_MPFR ... compute with mpfr library        */
/*                                (link with -lmpfr) */
/*   -DCHECK_CR ... with regorous check              */
/*            whether or not correctly rounded coefs */
/*                                                   */
/*   -DPRINT_WARNING ... print relative errors       */
/*                              and history of retry */
/*                                                   */
/* How to use:                                       */
/*   $ ./chevcoef degree correct_decimal_digits      */
/*                                                   */
/*****************************************************/
#include <stdio.h>
#include <math.h>
#include "gmp.h"

#ifdef USE_MPFR
	#include "mpfr.h"
	#include "mpf2mpfr.h"
#endif

// Maximum degree of polynomial
#define MAX_DEGREE 32768

// for Arrays to keep coefficients of polynomial
#define MAX_DEGREE_P1 (MAX_DEGREE + 1)

// Max macro
#define MAX(i, j) (((i) > (j)) ? (i) : (j))

// using IEEE754 double ... so, it is unuseful.
/* coef[deg] * x^deg + ... + coef[1] * x + coef[0] */
void get_chevcoef(double coef[], long int deg)
{
	long int i, j, index;
	double sum;

	coef[deg] = 1.0;
	for(i = 1; i <= deg / 2; i++)
	{
		coef[deg - (i * 2 - 1)] = 0.0;
		sum = 0;
		for(j = 1; j <= i; j++)
			sum += (1.0 / (2 * j + 1)) * coef[deg - 2 * (i - j)];
		sum *= -deg / (2 * i);
		coef[deg - i * 2] = sum;
	}
}

// 
/* coef[deg] * x^deg + ... + coef[1] * x + coef[0] */
void get_mpf_chevcoef(mpf_t coef[], long int deg)
{
	long int i, j;
	mpf_t sum, tmp;

	mpf_init2(sum, mpf_get_prec(coef[0]));
	mpf_init2(tmp, mpf_get_prec(coef[0]));

	mpf_set_ui(coef[deg], 1UL);
	for(i = 1; i <= deg / 2; i++)
	{
		mpf_set_ui(coef[deg - (i * 2 - 1)], 0UL);
		mpf_set_ui(sum, 0UL);
		for(j = 1; j <= i; j++)
		{
			mpf_set_ui(tmp, 1UL);
			mpf_div_ui(tmp, tmp, 2 * j + 1);
			mpf_mul(tmp, tmp, coef[deg - 2 * (i - j)]);
			mpf_add(sum, sum, tmp);
		}
		mpf_mul_ui(sum, sum, deg);
		mpf_div_ui(sum, sum, 2 * i);
		mpf_neg(sum, sum);
		mpf_set(coef[deg - i * 2], sum);
	}

	mpf_clear(sum);
	mpf_clear(tmp);
}

// log10(a)
void _mpf_log10(mpf_t ret, mpf_t a)
{
	mpf_t tmp;
	long int a_2exp;
	double a_double;

	mpf_init2(tmp, mpf_get_prec(a));

	// tmp = relerr * 2^relerr_2exp
#ifndef USE_MPFR
	a_double = mpf_get_d_2exp(&a_2exp, a);
	mpf_set_d(ret, log10(fabs(a_double)));
	mpf_set_d(tmp, (double)a_2exp * log10(2.0));
	mpf_add(ret, ret, tmp);
#else
	mpfr_log10(ret, a, GMP_RNDN);
#endif
	mpf_clear(tmp);
}


// CEE-based chevcoef 
void print_cee_chevcoef(long int deg, unsigned long dprec)
{
	unsigned long s, l, bs, bl, c, p[2], prec;
	long int stop_flag, i, relerr_2exp;
	double log10_relerr_max, relerr;

	mpf_t mpfcoef_s[MAX_DEGREE_P1], mpfcoef_l[MAX_DEGREE_P1];
	mpf_t tmp, rounded_s, rounded_l, tmp1;

	// deg <= 33768
	if((deg > MAX_DEGREE) || (dprec <= 0))
		return;

	// User-required precision
	prec = (unsigned long)ceil(dprec / log10(2.0));

#ifdef CHECK_CR
	mpf_init2(rounded_s, prec);
	mpf_init2(rounded_l, prec);
#endif

	mpf_init2(tmp1, prec);

	// precision
	log10_relerr_max = -(double)dprec;
	c = MAX(dprec / 10, 10);
	s = dprec + c;
	l = s + c;

	// additional precision
	p[0] = 0; p[1] = 0;

	// main loop
	do {

		stop_flag = 0;

		// to binary
		bs = (unsigned long)ceil(s / log10(2.0));
		bl = (unsigned long)ceil(l / log10(2.0));

		mpf_init2(tmp, bs);
		for(i = 0; i <= deg; i++)
		{
			mpf_init2(mpfcoef_s[i], bs);
			mpf_init2(mpfcoef_l[i], bl);
		}

		fprintf(stderr, "Additional digits: %d (%d bits)\n", MAX(p[0], p[1]), (unsigned long)ceil(MAX(p[0], p[1])/log10(2.0)));
		fprintf(stderr, "Short digits: %d (%d bits)\n", s, bs);
		get_mpf_chevcoef(mpfcoef_s, deg);
		fprintf(stderr, "Long  digits: %d (%d bits)\n", l, bl);
		get_mpf_chevcoef(mpfcoef_l, deg);

		// reset additional digits
		p[0] = 0; p[1] = 0;

		// check relative errors
		for(i = 0; i <= deg; i++)
		{
			mpf_reldiff(tmp, mpfcoef_s[i], mpfcoef_l[i]);
			mpf_abs(tmp, tmp);
			if(mpf_cmp_ui(tmp, 0UL) == 0)
				continue;
			else
				_mpf_log10(tmp, tmp);

			// check relative error
			if((mpf_cmp_ui(tmp, 0UL) != 0) && (mpf_cmp_d(tmp, log10_relerr_max) > 0))
			{
#ifdef PRINT_WARNING
				fprintf(stderr, "Too big relerr! :coef[%d] %7.1e < ", i, log10_relerr_max);
				mpf_out_str(stderr, 10, 3, tmp);
				fprintf(stderr, "\n");
#endif
				stop_flag++;

				// p[0] = max(log10(relerr(coef[i])) - log10(reltol))
				mpf_set_d(tmp1, log10_relerr_max);
				mpf_sub(tmp, tmp, tmp1);
				if(p[0] < mpf_get_ui(tmp))
					p[0] = mpf_get_ui(tmp);
			}

#ifdef CHECK_CR
			// check if correctly rounded 
			mpf_set(rounded_s, mpfcoef_s[i]);
			mpf_set(rounded_l, mpfcoef_l[i]);
			if(mpf_cmp(rounded_s, rounded_l) != 0)
			{
#ifdef PRINT_WARNING
				fprintf(stderr, "Not Correctly rounded!: coef[%d] reldiff = ", i);
				mpf_reldiff(tmp, rounded_s, rounded_l);
				mpf_out_str(stderr, 10, 2, tmp);
				printf("\n");
#endif
				stop_flag++;

				// p[1] = max(log10(relerr(coef[i])) - log10(reltol))
				mpf_set_d(tmp1, log10_relerr_max);
				_mpf_log10(tmp, tmp);
				mpf_sub(tmp, tmp, tmp1);
				if(p[1] > mpf_get_ui(tmp))
					p[1] = mpf_get_ui(tmp);
			}
#endif
		}

		// retry ?
		if(stop_flag <= 0)
			break;

		// reset
		c += c + MAX(p[0], p[1]);
		s += c;
		l += c;

		mpf_clear(tmp);
		for(i = 0; i <= deg; i++)
		{
			mpf_clear(mpfcoef_s[i]);
			mpf_clear(mpfcoef_l[i]);
		}

	} while(stop_flag > 0);

	// print & clear
	for(i = deg; i >= 0; i--)
	{
		printf("%d ", i);
		mpf_out_str(stdout, 10, dprec, mpfcoef_l[i]);
		printf("\n");

		mpf_clear(mpfcoef_l[i]);
		mpf_clear(mpfcoef_s[i]);
	}

	// clear
	mpf_clear(tmp);
	mpf_clear(tmp1);

#ifdef CHECK_CR
	mpf_clear(rounded_s);
	mpf_clear(rounded_l);
#endif

}

// usage
void usage(const char *progname)
{
	printf("USAGE: %s DEGREE DECIMAL_DIGITS\n", progname);
	printf("   ex: %s 100 50 -> 100th degree, 50 decimal digits\n", progname);
}

// Main function
int main(int argc, char *argv[])
{
	//
	if(argc <= 2)
	{
		usage(argv[0]);
		exit(0);
	}

	print_cee_chevcoef((long int)atoi(argv[1]), (unsigned long)atoi(argv[2]));

	return 0;
}