#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 2

void Vector(double*, double[N][N], double);
void Output(const char*, double*);
void YMT(double[N][N], double*, double*, double);
double Inaccuracy(double*, double, int);
int ComputingInaccuracy(double*, double, int);

int main() {
	double AD[N][N];
	double X0[N], X[N];
	double X_analit[N];
	double h[] = { 0.01, 0.001, 0.0001 };

	AD[0][0] = -1;
	AD[0][1] = 100;
	AD[1][0] = -100;
	AD[1][1] = -1;

	for (int i = 0; i < 3; i++) {
		printf("h = %f\n%dx%d LAE systems\nn = %d\n", h[i], N, N, N);
		Vector(X0, AD, 0.0);
		Output("X0", X0);
		YMT(AD, X0, X, h[i]);
		Vector(X_analit, AD, h[i]);
		printf("analitic value of X:\n");
		Output("X", X_analit);
		for (int j = 0; j < N; j++)	{
			printf("X[%d]:\neps=%32.16f\nepsComp=%11d\n", j, Inaccuracy(X, X_analit[j], j),	ComputingInaccuracy(X, X_analit[j], j));
		}
		printf("\n\n");
	}
	return 0;
}

double Inaccuracy(double* X, double Z, int i) {
	return fabs((X[i] - Z) / Z);
}

int ComputingInaccuracy(double* X, double Y, int k) {
	int i = 0, state = 0, y_i, z_i;
	double y, z;
	y = Y;
	z = X[k];
	y_i = (int)y;
	z_i = (int)z;
	while (y_i == z_i) {
		if (y_i && z_i && !state)
			state = 1;
		if ((y_i == z_i) && (y_i || (!y_i && state == 1)))
			i++;
		y = y * 10 - y_i * 10;
		z = z * 10 - z_i * 10;
		y_i = (int)y;
		z_i = (int)z;
	}
	return i;
}

void Vector(double* Z, double AD[N][N], double t) {
	Z[0] = (3*cos(100*t) + sin(100*t))*exp(-t);
	Z[1] = (cos(100*t) - 3*sin(100*t))*exp(-t);
}

void Output(const char* c, double* Y) {
	printf("vector %s:\n", c);
	for (int i = 0; i < N; i++) {
		printf("%sij = %32.16f\n", c, Y[i]);
	}
}

void YMT(double AD[N][N], double X0[N], double X[N], double h) {
	double X1[N];
	for (int i = 0; i < N; i++) {
		X1[i] = X0[i] + h * (AD[i][0] * X0[0] + AD[i][1] * X0[1]);
		X[i] = X0[i] + (h / 2) * (AD[i][0] * X0[0] + AD[i][1] * X0[1] + AD[i][0] * X1[0] + AD[i][1] * X1[0]);
	}
	Output("X", X);
}