#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 = 6;
	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;

	// srand( 783512 );  // 17  6.57e-17, 7.72e-17, 2.94e-16
	// srand( 95135 );  // 18  7.88e-17, 1.10e-16, 2.59e-16
	// srand( 9531254 );  // 19  8.38e-17, 1.50e-16, 5.61e-16

	double y = 1.0 + static_cast<double>( rand() )/RAND_MAX;
	cout << "[" << 0.0 << "," << y << "]," << endl;
	p.insert( 0.0, y );

	for ( int i = 0; i < N - 2; ++i ) {
		double x = static_cast<double>( rand() )/RAND_MAX/(N - 2) + static_cast<double>(i)/(N - 2);
		y = 1.0 + static_cast<double>( rand() )/RAND_MAX;
		cout << "[" << x << "," << y << "]," << endl;
		p.insert( x, y );
	}

	y = 1.0 + static_cast<double>( rand() )/RAND_MAX;
	cout << "[" << 1.0 << "," << y << "]," << endl;
	p.insert( 1.0, y );

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

	cout << "f := CurveFitting:-PolynomialInterpolation( pts, x ):" << 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;
}