#include <iostream.h>
#include <windows.h>
#include <math.h>

char bufRus[256];

char* Rus(const char* text) {
	CharToOem(text, bufRus);
	return bufRus;
}

float fun1(float* X) {
	float Y;
	Y = 3 * pow(X[0], 3) + 4 * pow(X[1], 2) - 2 * X[2];
	return (Y);
}

float fun2(float* X) {
	float Y;
	Y = 5 * X[2] * X[0] - 2 * X[1];
	return (Y);
}

float fun3(float* X) {
	float Y;
	Y = 4 * pow(X[0], 2) + 2 * pow(X[1], 3) - X[2];
	return (Y);
}

void Gauss(float** a, float* b, float* x, int n)
// Метод Гаусса с выбором главного элемента по столбцу
{
	float tmpValue;
	int i, j, k, z;
	float dblLeadElement;
// Прямой ход
	for (i = 0; i < n; i++) {
		dblLeadElement = a[i][i];
// Находим строку, в которой элемент, стоящий под ведущим - наибольший
		float tmpMax = dblLeadElement;
		int tmpMaxNumber = i;

		for (z = i; z < n; z++) {
			if (a[z][i] > tmpMax) {
				tmpMax = a[z][i];
				tmpMaxNumber = z;
			}
		}
// Меняем местами i-ю строку и строку tmpMaxNumber
		for (z = i; z < n; z++) {
			tmpValue = a[i][z];
			a[i][z] = a[tmpMaxNumber][z];
			a[tmpMaxNumber][z] = tmpValue;
		}
		tmpValue = b[i];
		b[i] = b[tmpMaxNumber];
		b[tmpMaxNumber] = tmpValue;
		dblLeadElement = tmpMax;

		for (j = i; j < n; j++) {
			a[i][j] /= dblLeadElement;
		}
		b[i] /= dblLeadElement;

		for (k = i + 1; k < n; k++) {
			float dblToDivide = a[k][i] / a[i][i];
			for (z = i; z < n; z++) {
				a[k][z] -= a[i][z] * dblToDivide;
			}
			b[k] -= b[i] * dblToDivide;
		}
	}
// Обратный ход
	x[n - 1] = b[n - 1];
	for (k = n - 2; k >= 0; k--) {
		float sum = b[k];
		for (j = k + 1; j < n; j++) {
			sum -= a[k][j] * x[j];
		}
		x[k] = sum;
	}
}

void Newton(int N, float** A, float* B, float* H, float* X, float E) {
//Метод Ньютона
	float xa;
	bool R;
	do {
		R = 0;
		B[0] = -fun1(X);
		B[1] = -fun2(X);
		B[2] = -fun3(X);
		for (int j = 0; j < N; j++) {
			xa = X[j];
			X[j] = xa + E;
			for (int i = 0; i < N; i++) {
				if (i == 0)
					A[i][j] = (fun1(X) + B[i]) / E;
				if (i == 1)
					A[i][j] = (fun2(X) + B[i]) / E;
				if (i == 2)
					A[i][j] = (fun3(X) + B[i]) / E;
			}
			X[j] = xa;
		}
		Gauss(A, B, H, N);
		for (int j = 0; j < N; j++) {
			if (fabs(H[j]) >= E)
				R = 1;
			X[j] = X[j] + H[j];
		}
	} while (R == 1);
}

int main() {
	system("cls");
	int N;
	float E;
	cout << Rus("Введите порядок системы уравнений: ");
	cin >> N;
	cout << Rus("Введите показатель степени точности: ");
	cin >> E;
	E = pow(10, E);
	float** A;
	A = (float**) malloc(N * sizeof(float*));
	float* B;
	B = (float*) malloc(N * sizeof(float));
	float* H;
	H = (float*) malloc(N * sizeof(float));
	float* X;
	X = (float*) malloc(N * sizeof(float));
	float* Y;
	Y = (float*) malloc(N * sizeof(float));
	for (int i = 0; i < N; i++) {
		A[i] = (float*) malloc(N * sizeof(float));
		cout << Rus("Введите начальное приближение X в уравнении ") << i + 1
				<< ": ";
		cin >> X[i];
	}
	Newton(N, A, B, H, X, E);
	Y[0] = fun1(X);
	Y[1] = fun2(X);
	Y[2] = fun3(X);
	cout << Rus("Корни системы уравнений:") << endl;
	for (int i = 0; i < N; i++)
		cout << "X[" << i + 1 << "]=" << X[i] << endl;
	cout << Rus("Проверка:") << endl;
	for (int i = 0; i < N; i++)
		cout << Rus("Y[") << i + 1 << "]=" << Y[i] << endl;
	free(A);
	free(B);
	free(H);
	free(X);
	free(Y);
	cin.sync();
	cin.get();
	return 0;
}

float f1(float* x)
{return pow(x[0],2)+pow(x[1],2)+pow(x[2],2)-1;}

float f2(float* x)
{return 2*pow(x[0],2)+pow(x[1],2)-4*x[2];}

float f3(float* x)
{return 3*pow(x[0],2)-4*x[1]+pow(x[2],2);}