#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() {
	Newton_polynomial<5> p;
	cout.precision( 20 );
	int const N = 5;

	cout << p << endl;

	for ( double x = 1; x <= 2*N; ++x ) {
		p.insert( x, std::sin(x) );
		cout << "p := " << p << ":" << endl;

		double low = std::max( 1.0, x - N + 1 );
		double high = x;

		cout << "L := [";

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

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

		cout << "plots[display]( plots[pointplot]( L, symbol = circle ), plot( p, x = " << low << ".." << high << " ), plot( sin(x), x = " << low << ".." << high << ", colour = blue ) );" << endl << endl;
	}

	return 0;
}