#define _CRT_SECURE_NO_WARNINGS
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "header.h"

double integrate_left_rects(double(*f)(double), double prev_S, long int N, double a, double b, int* f_counter) {
	double S;
	long double H = (long double)(b - a) / N;

	S = 0;
	for (int i = 2;i <= N;i += 2) {
		S += f(a + (i - 1) * H, f_counter);
	}

	return H * S + (0.5) * prev_S;
}
void runge_left_rects(double eps, double(*f)(double), double a, double b,int* f_counter) {
	double S, S2 = 0;
	long int N = 1;

	S2 = (b - a) * f(0,f_counter);

	do {
		S = S2;
		N *= 2;
		S2 = integrate_left_rects(f, S, N, a, b, f_counter);

	} while (fabs(S2 - S) > eps);
	return S2;
}



double polyval(double* coeffs, double x, int n) {
	double ans = 0;
	for (int i = 0;i < n + 1;i++) {
		ans += coeffs[i] * pow(x, n - i);
	}
	return ans;
}
double* polyder(double* coeffs, int n) {
	double* der = malloc(sizeof(double) * (n));
	for (int i = 0;i < n;i++) {
		der[i] = coeffs[i] * (n - i);
	}
	return der;
}
double newton_solve(double* coeffs, int n, double a, double b) {
	double x_next = b, x_prev = b;
	double* der = polyder(coeffs, n);
	do {
		x_prev = x_next;
		x_next = x_prev - polyval(coeffs, x_prev, n) / polyval(der, x_prev, n - 1);
	} while (fabs(x_next - x_prev) >= 10e-10);
	free(der);
	return x_next;
}
double* solve_poly(double* coeffs, int n, double a, double b) {
	double* xx = calloc(n, sizeof(double));
	int c = 0;
	double x1 = a, x2 = a;
	while (x1 < b) {
		while (x2 < x1) {
			if (polyval(coeffs, x1, n) * polyval(coeffs, x2, n) < 0) {
				xx[c] = newton_solve(coeffs, n, x2, x1);
				c++;
				x1 = x2;
				x2 += (b - a) / 10000;
				break;
			}
			x2 += (b - a) / 10000;
		}
		x1 += (b - a) / 10000;
	}
	return xx;
}
double** table_xx() {
	double** ret = malloc(sizeof(double*) * 4);

	for (int n = 2;n <= 5;n++) {
		double* coeffs = malloc(sizeof(double) * (n + 1));
		double* sk = malloc(sizeof(double) * (n + 1));
		if (coeffs == NULL || sk == NULL) {
			exit(16);
		}
		coeffs[0] = 1;

		double A = 2.0 / n;
		for (int i = 1;i < (n + 1);i++) {
			sk[i] = (1.0 / (i + 1) - pow(-1, i + 1) / (i + 1)) / A;
			coeffs[i] = sk[i];
			for (int j = 1;j < i;j++) {
				coeffs[i] += coeffs[j] * sk[i - j];
			}
			coeffs[i] *= -1.0 / i;
		}

		ret[n - 2] = solve_poly(coeffs, n, -1, 1);

		free(sk);
		free(coeffs);
	}
	return ret;
}


double chebyshev_method(double(*f)(double), double a1, double b1, int n, double** xx_tab, int* f_counter) {
	double s = 0, a, b, ans;

	a = a1;
	b = b1;
	double* xx = xx_tab[n - 2];
	double A = 2.0 / n;

	ans = 0;
	for (int k = 0;k < n;k++) {
		ans += f((a + b) / 2 + ((b - a) / 2) * xx[k], f_counter);
	}
	s = ((b - a) / 2) * ans * A;



	return s;
}


double adaptive_div(double(*f)(double), double** xx_tab, double a, double b, int n, double eps, double S_n, double S_n_plus_1, int* N, int* f_counter) {

	if (fabs(S_n_plus_1 - S_n) < eps / (*N)) {
		return S_n;
	}
	else {
		double S_n_plus_1_r = chebyshev_method(f, a, (a + b) / 2, n + 1, xx_tab, f_counter);
		double S_n_r = chebyshev_method(f, a, (a + b) / 2, n, xx_tab, f_counter);
		double S_n_plus_1_l = chebyshev_method(f, (a + b) / 2, b, n + 1, xx_tab, f_counter);
		double S_n_l = chebyshev_method(f, (a + b) / 2, b, n, xx_tab, f_counter);
		*N += 1;

		return (adaptive_div(f, xx_tab, a, (a + b) / 2, n, eps, S_n_l, S_n_plus_1_l,N, f_counter) + adaptive_div(f, xx_tab, (a + b) / 2, b, n, eps, S_n_r, S_n_plus_1_r, N, f_counter));
	}

}



double cheb_approx(double(*f)(double), double a, double b, int n, double eps,int* f_counter) {
	
	double** xx_table = table_xx(a, b);
	int N = 1;

	return adaptive_div(f, xx_table, a, b, n, eps, 100, 1,&N, f_counter);
}

double* counters(double(*f)(double), double a, double b, double eps) {
	int counter_recs = 0,counter_cheb = 0;

	cheb_approx(f, a, b, 3, eps, &counter_cheb);
	runge_left_rects(eps, f, a, b, &counter_recs);
	printf("eps = %e, counter_recs = %d, counter_cheb = %d\n",eps, counter_recs, counter_cheb);
	int* counter = malloc(sizeof(int) * 2);
	counter[0] = counter_recs;
	counter[1] = counter_cheb;
	
	return counter;
}