// decent.cpp : Defines the entry point for the console application.
//

//#include "stdafx.h"
#include <iostream>
#include <cmath>
using namespace std;

int n = 2;
double e = 0.0001;
double K = 8;
int max_iter = 1000;

double f(double* x, const int n)
{
	double y = pow(x[0], 4) + 2*K*pow(x[0], 3) - 2*pow(K, 2)*pow(x[0], 2) - 2*pow(K, 3)*x[0] + 
	+ pow(x[1], 4) - pow(x[1], 2)*(2*pow(K, 2) + 2*K) + 
	2*pow(K, 4) + pow(K, 3) + pow(K, 2);
	return y;
}
double* grad(double* x, int n)
{
	double* g = new double[n];
	g[0] = 4*pow(x[0], 3) + 6*K*pow(x[0], 2) - 4*pow(K, 2)*x[0] - 2*pow(K, 3);
	g[1] = 4*pow(x[1], 3) - 2*(2*pow(K, 2) + 2*K)*x[1];
	return g;
}
double goldens(double _a, double _b, double* tmp1, double* tmp2, int ia, int ib)
{
	double a = _a; double b = _b;
	double phi = (sqrt(5.0) + 1)/2;
	double x1 = b - (b-a)/phi;
	double x2 = a + (b-a)/phi;
	double min_l;
	int iter = 0;
	while(!(b - a < e && b - a > - e) && iter < max_iter){
		iter++;
		tmp1[ia] = x1; tmp2[ib] = x2;
		if(f(tmp1, n) < f(tmp2, n)){
			b = x2;
			x2 = x1;
			x1 = b - (b-a)/phi;
		}
		else{
			a = x1;
			x1 = x2;
			x2 = a + (b-a)/phi;
		}
		min_l = (f(tmp1, n)<f(tmp2, n))?x1:x2;
		cout.width(4); cout.setf(ios::left);
		cout << iter;
		for(int i = 0; i < n; i++){
			cout.width(10);
			cout << ((x1==min_l)?tmp1[i]:tmp2[i]) << " ";
		}
	cout << endl;
	}
	return min_l;
}
double goldens(double _a, double _b, double* gr, double * x0, double* tmp1, double* tmp2)
{
	double a = _a; double b = _b;
	double phi = (sqrt(5.0) + 1)/2;
	double x1 = b - (b-a)/phi;
	double x2 = a + (b-a)/phi;
	double min_l;
	int iter = 0;
	while( !(b - a < e && b - a > -e) && iter < max_iter ){
		iter++;
		for(int i = 0; i < n; i++){
			tmp1[i] = x0[i] + x1*gr[i];
			tmp2[i] = x0[i] + x2*gr[i];
		}
		if(f(tmp1, n) < f(tmp2, n)){
			b = x2;
			x2 = x1;
			x1 = b - (b-a)/phi;
		}
		else{
			a = x1;
			x1 = x2;
			x2 = a + (b-a)/phi;
		}

		min_l = (f(tmp1, n)<f(tmp2, n))?x1:x2;
		/*for(int i = 0; i < n; i++)
			cout << ((x1==min_l)?tmp1[i]:tmp2[i]) << " ";
		cout << endl;*/
	}
	return min_l;
}
double mps(double* x0, double a, double b){
	double* xk = new double[n];
	double* xk1 = new double[n];
	for(int i=0; i < n; i++)
		xk[i] = xk1[i] = x0[i];
	for(int i = 0; i < n; i++){
		for(int j = 0; j < n; j++)
			xk[j] = xk1[j];
		xk1[i] = goldens(a, b, xk, xk1, i, i);
		}
	//delete tmp1, tmp2;
	return f(xk1, n);
}
double gm(double* x0, double a, double b){
	double* xk = new double[n];
	double* xk1 = new double[n];
	double* tmp1 = new double[n];
	double* tmp2 = new double[n];
	for(int i = 0; i < n; i++){
		xk1[i] = xk[i] = x0[i];
	}
	int iter = 0;
	do
	{
		iter++;
		for(int i = 0; i < n; i++)
			xk[i] = xk1[i];
		double* gr = grad(xk1, n);
		for(int i = 0; i < n; i++)
			tmp1[i] = tmp2[i] = xk[i];
		double lamda = goldens(a, b, gr, xk, tmp1, tmp2);
		cout.width(4); cout.setf(ios::left);
		cout << iter;
		for(int i = 0; i < n; i++){
			xk1[i] = xk[i] + lamda*gr[i];
			cout.width(10);
			cout << xk1[i] << " ";
		}
		cout << endl;
		delete gr;
	}
	while(!((f(xk1, n) - f(xk, n) < e) && (f(xk1, n) - f(xk, n) > -e)) && iter < max_iter);
	return f(xk1, n);
}
//int _tmain(int argc, _TCHAR* argv[])
int main(int argc, char* argv[])
{
	if(argc != n+1){ cout << "Usage: descent x1 x2" << endl; return 0; }
	if ( (!atof(argv[1]) && argv[1][0] != '0' ) && (!atof(argv[2]) && argv[2][0] != '0' ) ) 
	{
		cout << "Wrong input" << endl;
		return 0;  
	}
	cout << "x^4 + 2*k*x^3 - 2*k^2*x^2 - 2^3*x + y^4 - (2*k^2 + 2*k)*y^2 + 2*k^4 + k^3 + k^2 -> min" << endl;
	cout << "k = 8" << endl;

	double* x0 = new double[n];
	for(int i = 0; i < n; i++){
		x0[i] = atof(argv[i+1]);
	}
	cout << "MPS" << endl;
	cout <<"Final: " << mps(x0, -1000, 1000) << endl;
	cout << "Gradient" << endl;
	cout << "Final: " << gm(x0, -1000, 1000) << endl;

	return 0;
}