#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gmp.h>

void
taylorSin(mpf_t y, const mpf_t x, const int degree)
{
	mpf_t top, bottom_f, term;
	mpz_t bottom_z;

	mpf_set_str(y, "0", 10);

	mpf_inits(top, bottom_f, term, NULL);
	mpz_inits(bottom_z, NULL);

	for (int n = 0; n <= degree; n++) {
		mpf_pow_ui(top, x, (unsigned long int)(2 * n + 1));

		if (n % 2 == 1) {
			mpf_neg(top, top);
		}

		mpz_fac_ui(bottom_z, (unsigned long int)(2 * n + 1));
		mpf_set_z(bottom_f, bottom_z);

		mpf_div(term, top, bottom_f);
		mpf_add(y, y, term);
	}

	mpf_clears(top, bottom_f, term, NULL);
	mpz_clears(bottom_z, NULL);
}

void
firstDerivativeTaylorSin(mpf_t y, const mpf_t x, const int degree)
{
	mpf_set_str(y, "0", 10);

	mpf_t     h, h2, f1, f2, x1, x2;
	mpf_inits(h, h2, f1, f2, x1, x2, NULL);

	mpf_set_str(h,  "0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001", 10);
	mpf_set_str(h2, "0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002", 10);

	mpf_add(x1, x, h);
	mpf_sub(x2, x, h);

	taylorSin(f1, x1, degree);
	taylorSin(f2, x2, degree);

	mpf_sub(y, f1, f2);
	mpf_div(y, y, h2);
	
	mpf_clears(h, h2, f1, f2, x1, x2, NULL);
}

void
rootTaylorSin(mpf_t x, int degree, int iterations)
{
	mpf_t y, dydx, temp;
	mpf_inits(y, dydx, temp, NULL);

	for (int n = 0; n < iterations; n++) {
		taylorSin(y, x, degree);
		firstDerivativeTaylorSin(dydx, x, degree);

		mpf_div(temp, y, dydx);
		mpf_sub(x, x, temp);

		gmp_printf("%.100Ff\n", x);
	}

	mpf_clears(y, dydx, temp, NULL);
}

int
main(int argc, char *argv[])
{
	mpf_set_default_prec(65536);

	mpf_t result;
	mpf_t x;

	mpf_inits(result, x, NULL);

	mpf_set_str(x, "3", 10);

	/*taylorSin(result, x, 10);
	gmp_printf("%.100Ff\n", result);
	firstDerivativeTaylorSin(result, x, 10);
	gmp_printf("%.100Ff\n", result);*/
	rootTaylorSin(x, 3000, 16);
	gmp_printf("%.16000Ff\n", x);

	mpf_clears(result, x, NULL);

	return 0;
}