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


double* create_grid_even(int n, double a, double b) {
	double* h = malloc(sizeof(double) * n);
	if (h == NULL) {
		exit(16);
	}
	for (int i = 0; i < n;i++) {
		h[i] = a + ((b - a) / (n - 1)) * i;
	}
	return h;
}

double* euler_method(double y0, double a, double b, double(*f)(double,double),int n) {
	//(x[0], y0) - первая, имеющаяся у нас точка табличной функции
	double h = (b-a)/(n-1);
	double* y = malloc(sizeof(double) * n);
	if (y == NULL) {
		exit(16);
	}
	y[0] = y0;
	for (int i = 0;i < n-1;i++) {
		y[i + 1] = y[i] + h * f(a + i*h, y[i]);
	}
	return y;
}

double max_dif(double* r2, double* r1,int n) {
	double m = 0;
	double dif;
	for (int i = 0;i < (int)n/2-1;i++) {
		dif = fabs(r2[2 * i + 2] - r1[i + 1]);
		if (m < dif) {
			m = dif;
		}
	}
	return m;
}

double** runge(double eps, double y0, double a, double b, double(*f)(double, double), int* n_, int* len_h) {
	double* x = malloc(sizeof(double));
	double* y = malloc(sizeof(double));
	long double* h_all = malloc(sizeof(long double));// регистрирует все значения шагов
	long double* h_made = malloc(sizeof(long double));// регистрирует значения шага, соответствующие сделаным шагам
	long double h = (b - a) / 2;

	y[0] = y0;
	x[0] = a;
	h_all[0] = h;
	double y_i_1 = 0;
	double y_i_2 = 0;
	double y_i_2_2 = 0;

	int k = 1, i = 1;
	double f_prev,dif;

	printf("eps = %e\n", eps);
	
	eps = pow(eps, (double)3 / 2);
	while (x[i-1] + h <= b) {
		h_all = realloc(h_all, sizeof(long double) * (k + 1));
		h_all[k] = h;

		// подход, который я оставил ниже в комменте для достижения заданной точности не совсем подходит: для относительно большого значения шага значения табличной функции почти не меняются

		/*f_prev = f(a + (i - 1) * h, y[i - 1]);

		y_i_1 = y[i - 1] + h * f_prev;
		y_i_2 = y[i - 1] + (h / 2) * f_prev;
		y_i_2_2 = y_i_2 + (h / 2) * f(a + (i - 1 / 2) * h, y_i_2);*/

		// вместо него пользуюсь подходом ниже: сравниваю значение решения метода без разбиения и с разбиением на 9 точек (h -> h/8) 
		y_i_1 = euler_method(y[i - 1], x[i-1], x[i-1]+h, f, 2)[1];
		y_i_2_2 = euler_method(y[i - 1], x[i - 1], x[i - 1] + h, f, 9)[8];

		dif = fabs(y_i_2_2 - y_i_1);

		if (dif > eps) {
			// также, для ускорения работы метода h делю на 8 в случае невыполнения условия
			h = h/8;
		}


		else if (dif < eps/100) {
			
			x = realloc(x, sizeof(double) * (i + 1));
			x[i] = x[i - 1] + h;
			y = realloc(y, sizeof(double) * (i + 1));
			y[i] = y_i_2_2;
			
			h_made = realloc(h_made, sizeof(long double) * (i + 1));
			h_made[i - 1] = h;
			if (x[i] + h * 2 <= b) {
				h *= 2;
			}
			i += 1;

		
		}//шаг выполнен, удваивается


		else if(dif < eps) {
			x = realloc(x, sizeof(double) * (i + 1));
			x[i] = x[i-1] + h;
			

			y = realloc(y, sizeof(double) * (i + 1));
			y[i] = y_i_2_2;

			h_made = realloc(h_made, sizeof(long double) * (i + 1));
			h_made[i - 1] = h;

			i += 1;
		}//шаг выполнен, не изменяется
		k++;

	} 
	h_all[k-1] = h;
	h_made[i - 1] = h;
	*len_h = k;
	*n_ = i;
	double** ans = malloc(sizeof(double*) * 4);
	ans[0] = x;
	ans[1] = y;
	ans[2] = h_all;
	ans[3] = h_made;

	return ans;
}

double uglify_initial(double initial, int max_procent,int min_procent) {
	double procent = min_procent + (double)(max_procent-min_procent)*((double)rand()/(RAND_MAX));

	return initial - procent * initial / 100;
}