#include <iostream>
#include <cmath>
using namespace std;

#include "Newton_polynomial.h"
using namespace Data_structures;

/*
 *  This function prints out Maple code which plots
 *  the function sin(x), the polynomial interpolating
 *  sin(x) at up to five points, and the Newton
 *  polynomial evaluated at various points on the
 *  interpolated range: [low, high].
 *
 *  The output can be cut-and-pasted directly into Maple.
 */

int main() {
	int const N = 5;
	Newton_polynomial<N> p;
	cout.precision( 20 );

	double const dx = 0.00006103515625;
	double const low =  0.0;
	double const high = 1.0;

	cout << "Digits := 100:" << endl;
	cout << "pts := [" << endl;

	// 8.92e-17, 9.68e-17, 1.03e-16

	for ( int i = 0; i < N; ++i ) {
		double x = static_cast<double>( i ) / static_cast<double>( N - 1 );
		double y = std::exp( x );
		cout << "[" << x << "," << y << "]," << endl;
		p.insert( x, y );
	}

	cout << "NULL]:" << endl;

	cout << "f := CurveFitting:-PolynomialInterpolation( pts, x ):" << endl;

	/////////////////////////////
	// The Expanded Polynomial //
	/////////////////////////////

	cout << "Lopt := [";

	for ( double xi = low; xi <= high; xi += dx ) {
		cout << "[" << xi << "," << (0.069415674220267*(xi*xi)*(xi*xi) + 0.140276004141752*(xi*xi)*xi + 0.509787138019482*(xi*xi) + 0.998803012077544*xi + 1) << "],";
	}

	cout << "NULL]:" << endl;

	cout << "N := nops( Lopt ):" << endl;

	cout << "Eopt := [seq( [i[1], abs( i[2] - eval( f, x = i[1] ) )], i = Lopt )]:" << endl;
	cout << "LL := plots[pointplot]( Eopt, style = point, symbol = box, color = black, symbolsize = 3 ):" << endl;
	cout << "MSEopt := add( Eopt[i][2]^2, i = 1..N ):" << endl;

	/////////////////////////////////////////////////
	// The Optimal Interpolating Newton Polynomial //
	/////////////////////////////////////////////////

	cout << "Lopt := [";

	for ( double xi = low; xi <= high; xi += dx ) {
		cout << "[" << xi << "," << p.evaluate( xi, OPTIMAL ) << "],";
	}

	cout << "NULL]:" << endl;

	cout << "N := nops( Lopt ):" << endl;

	cout << "Eopt := [seq( [i[1], abs( i[2] - eval( f, x = i[1] ) )], i = Lopt )]:" << endl;
	cout << "LL := plots[pointplot]( Eopt, style = point, symbol = box, color = black, symbolsize = 3 ):" << endl;
	cout << "MSEopt := add( Eopt[i][2]^2, i = 1..N ):" << endl;

	/////////////////////////////////////////////////
	// The Forward Interpolating Newton Polynomial //
	/////////////////////////////////////////////////

	cout << "Lfwd := [";

	for ( double xi = low; xi <= high; xi += dx ) {
		cout << "[" << xi << "," << p.evaluate( xi, FORWARD ) << "],";
	}

	cout << "NULL]:" << endl;

	cout << "Efwd := [seq( [i[1], abs( i[2] - eval( f, x = i[1] ) )], i = Lfwd )]:" << endl;
	cout << "LL := plots[pointplot]( Efwd, style = point, symbol = cross, color = aquamarine, symbolsize = 3 ), LL:" << endl;
	cout << "MSEfwd := add( Efwd[i][2]^2, i = 1..N ):" << endl;

	//////////////////////////////////////////////////
	// The Backward Interpolating Newton Polynomial //
	//////////////////////////////////////////////////

	cout << "Lbak := [";

	for ( double xi = low; xi <= high; xi += dx ) {
		cout << "[" << xi << "," << p.evaluate( xi, BACKWARD ) << "],";
	}

	cout << "NULL]:" << endl;

	cout << "Ebak := [seq( [i[1], abs( i[2] - eval( f, x = i[1] ) )], i = Lbak )]:" << endl;
	cout << "LL :=  plots[pointplot]( Ebak, style = point, symbol = diagonalcross, color = red, symbolsize = 3 ), LL:" << endl;
	cout << "MSEbak := add( Ebak[i][2]^2, i = 1..N ):" << endl;
	cout << "plots[display]( plot( f, x = " << low << ".." << high << ", colour = blue ), plots[pointplot]( pts, symbol = circle ) );" << endl << endl;
	cout << "plots[display]( LL );" << endl;
	cout << "plots[display]( LL );" << endl;
	cout << "plots[display]( LL, view = [DEFAULT, 0..1e-15] );" << endl;
	cout << "plots[display]( LL, view = [DEFAULT, 0..1e-15] );" << endl;
	cout << "sqrt( MSEopt/N ), sqrt( MSEfwd/N ), sqrt( MSEbak/N );" << endl;

	return 0;
}