#include <iostream>
#include <cmath>
#include <GL\glut.h>
#include <GL\GL.h>
#include <vector>

using namespace std;

double f(double x, double y)
{
	double f;
	f = y*y - 2 * x*y - 3 * x*x;
	return f;
}

double g(double x, double y)
{
	double g;
	g = y*y - x*y - 2 * x*x - 4;
	return g;
}

double f_diff_x(double x, double y)
{
	double f_diff_x;
	f_diff_x = -2 * y - 6 * x;
	return f_diff_x;
}

double f_diff_y(double x, double y)
{
	double f_diff_y;
	f_diff_y = 2 * y - 2 * x;
	return f_diff_y;
}

double g_diff_x(double x, double y)
{
	double g_diff_x;
	g_diff_x = -y - 4 * x;
	return g_diff_x;
}

double g_diff_y(double x, double y)
{
	double g_diff_y;
	g_diff_y = 2 * y - x;
	return g_diff_y;
}

void GetMinor(vector <vector<double>> A, vector <vector<double>> &B, int x, int y, int size)
{
	int _x = 0;
	int _y;

	for (int i = 0; i < size; ++i)
	{
		if (i != x)
		{
			_y = 0;
			for (int j = 0; j < size; ++j)
			{
				if (j != y)
				{
					B[_x][_y] = A[i][j];
					_y++;
				}
			}
			_x++;
		}
	}
}

//Finding determinant of matrix (size x size)
double FindDeterminant(vector <vector<double>> A, int size)
{
	if (size == 1)
	{
		return A[0][0];
	}
	else
	{
		double det = 0;
		vector <vector<double>> Minor(size - 1, vector <double>(size - 1));
		for (int i = 0; i < size; ++i)
		{
			GetMinor(A, Minor, 0, i, size);
			det += pow(-1, i) * A[0][i] * FindDeterminant(Minor, size - 1);
		}
		return det;
	}
}

double InsertFreeCoeff(vector <vector<double>> A, vector<double> B, int size, int j)
{
	for (int i = 0; i < size; ++i)
	{
		A[i][j] = B[i];
	}
	double det = FindDeterminant(A, size);
	return det;
}

const int N = 10;

double  x[N] = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
double y[N] = { -2.5, - 6.7, - 7.9, - 8.5, - 6.7, - 6.3, - 3.3, 0.3, 4.0, 9.2 };
int t = 20;
void mnk()
{
	double *a = NULL, *b = NULL, *roots = NULL, **sum = NULL;
	int K, i, j, k;

	printf("Xi \t   \t Yi\n");
	for (int i = 0; i < N; i++)
	{
		printf("%.1f \t | \t %.1f\n", x[i], y[i]);
	}
	cout << endl <<"rank k= ";
	cin >> K;
	cout << endl;
	b = new double[K + 1];
	a = new double[K + 1];
	sum = new double *[K + 1];
	for (int i = 0; i<K + 1; i++)
		sum[i] = new double[K + 1];


	for (i = 0; i<K + 1; i++)
	{
		for (j = 0; j<K + 1; j++)
		{
			sum[i][j] = 0;
			for (k = 0; k<N; k++)
			{
				sum[i][j] += pow(x[k], i + j);
			}
		}
	}
	for (i = 0; i<K + 1; i++)
	{
		b[i] = 0;
		for (k = 0; k<N; k++)
		{
			b[i] += pow(x[k], i) * y[k];
		}
	}
	for (int i = 0; i<K + 1; i++)
	{
		for (int j = 0; j < K + 1; j++)
		{
			if (j < K)
			{
				cout << "A" << j << "*(" << sum[i][j] << ") \t+ ";
			}
			else
			{
				cout << "A" << j << "*(" << sum[i][j] << ") \t ";
			}
		}
			
		cout << " = " << b[i];
		cout << endl;
	}
	cout << endl;
	//пихаем все в вектор для дин. использования
	vector <vector<double> > A(K + 1, vector <double>(K + 1));
	vector <double> B(K + 1);
	roots = new double[K + 1];
	for (int i = 0; i < K + 1; i++)
	{
		for (int j = 0; j < K + 1; j++)
		{
			A[i][j] = sum[i][j];
		}
		B[i] = b[i];
	}

	//решаем методом Крамера
	double det = FindDeterminant(A, K + 1);
	if (det == 0)
	{
		cout << ">>>No roots!<<<" << endl;
	}
	else
	{
		cout << "det_A: " << det << endl;
		cout << "\n";

		for (int i = 0; i < K + 1; ++i)
		{
			printf("det_%.d = %.2f\n", i+1, InsertFreeCoeff(A, B, K + 1, i));
		}
		cout << endl;
		cout << "Coefficients: " << endl;
		for (int i = 0; i < K + 1; ++i)
		{
			roots[i] = InsertFreeCoeff(A, B, K + 1, i) / det;
			cout << "A" << i  << " = " << roots[i] << "\n";
		}
		cout << endl << "Final equation: " << endl;
		cout << "F(X) = ";
		for (int i = 0; i < K + 1; ++i)
		{
			if (i < K)
			{
				printf("(%.3f)*X^%d + ", roots[i], i);
			}
			else
			{
				printf("(%.3f)*X^%d = 0\n", roots[i], i);
			}
			
		}

	}
	double val[50];
	double q = -5;
	int M = 50;
	for (int i = 0; i < M; i++)
	{
		val[i] = 0;
		for (int p = 0; p < K+1; p++)
		{
			val[i] += roots[p] * pow((double)q, p);
		}
		q=q+0.5;
	}
	
	q = -5;
	glLineWidth(2);
	glBegin(GL_LINE_STRIP);
	glColor3f(1.0, 0.0, 0.0);
	for (int i = 0; i <= M; i++)
	{
		glVertex2d(q * t, val[i] * t);
		q = q + 0.5;
	}
	glEnd();
}


void display()
{
	//координатная плоскость
	glColor3f(0.0, 1.0, 0.0);
	glBegin(GL_LINES);
	glVertex2f(-550, 0);
	glVertex2f(550, 0);
	glEnd();

	glBegin(GL_LINES);
	glVertex2f(0, 150);
	glVertex2f(0, -150);
	glEnd();
	glPointSize(2.0);
	for (int i = -550; i < 550; i = i + 10)
	{
		if (i == 0)
		{
			glPointSize(4.0);
			glBegin(GL_POINTS);
			glColor3f(1.0, 1.0, 1.0);
			glVertex2d(i, 0);
			glEnd();
		}
		else
		{
			glPointSize(2.0);
			glBegin(GL_POINTS);
			glColor3f(1.0, 1.0, 1.0);
			glVertex2d(i, 0);
			glEnd();
		}
	}
	for (int i = -150; i < 150; i = i + 10)
	{
		if (i == 0)
		{
			glPointSize(4.0);
			glBegin(GL_POINTS);
				glColor3f(1.0, 1.0, 1.0);
				glVertex2d(0, i);
			glEnd();
		}
		else
		{
			glPointSize(2.0);
			glBegin(GL_POINTS);
				glColor3f(1.0, 1.0, 1.0);
				glVertex2d(0, i);
			glEnd();
		}
	}
	
	//отмечаем точки и соединяем
	glPointSize(5.0);
	for (int i = 0; i < N; i++)
	{
		glBegin(GL_POINTS);
		glColor3f(0.0, 0.0, 1.0);
		glVertex2d(x[i] * t, y[i] * t);
		glEnd();
	}

	glBegin(GL_LINE_STRIP);
	glColor3f(1.0, 0.0, 1.0);
	for (int i = 0; i <N; i++)
	{
		glVertex2d(x[i] * t, y[i] * t);
	}
	glEnd();

	mnk();
	
	glFlush();
}

int main(int argc, char **argv) {
	glutInit(&argc, argv); 
	glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB); 
	glutInitWindowSize(600*2, 165*2);
	glutInitWindowPosition(0, 0);
	glutCreateWindow("Sin");
	glClearColor(1.0, 1.0, 1.0, 1.0);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(-600.0, 600.0, -200.0, 200.0, -1.0, 1.0);
	glutDisplayFunc(display);
	glutMainLoop();
	return 0;
}