#ifndef CA_UWATERLOO_ALUMNI_DWHARDER_NEWTON_POLYNOMIALS
#define CA_UWATERLOO_ALUMNI_DWHARDER_NEWTON_POLYNOMIALS

#include <cassert>
#include <ostream>
#include <cmath>

// Author:  Douglas Wilhelm Harder
// Copyright (c) 2010 by Douglas Wilhelm Harder.  All rights reserved.
// Newton polynomial divided difference table with evaluation using
// dynamically constructed polynomials using Horner's rule.

namespace Data_structures {
	/****************************************************
	 * ************************************************ *
	 * *              Newton Polynomials              * *
	 * ************************************************ *
	 ****************************************************/

	/*
	 *  Newton Polynomials
	 *
	 *  This is a class for allowing the user to create
	 *  a Newton polynomial of degree up to N.  The
	 *  user adds pairs (x[k], y[k]) and this algorithm
	 *  will store the divided difference table with the
	 *  coefficients of the Newton polynomial passing
	 *  through the last N points.
	 *
	 *  If N points have already been added and an
	 *  (N + 1)st point is inserted, the oldest point
	 *  is removed.
	 *
	 *  Memory:    O(N^2)
	 */

	enum direction {FORWARD, BACKWARD, OPTIMAL};

	template <int N, typename Type = double>
	class Newton_polynomial {
		public:
			Newton_polynomial();
			Type evaluate( Type const &, direction = OPTIMAL ) const;
			void insert( Type const &, Type const & );
			Type x( int k ) const;
			Type y( int j, int k ) const;
			void clear();

		private:
			Type x_values[N];
			Type table[N][N];
			int count;
			int begin;
			int end;

			static inline int next( int );
			static inline int previous( int );
			static inline int remainder( int );

		template <int M, typename T>
		friend std::ostream &operator<<( std::ostream &, const Newton_polynomial<M, T> & );

	};

	/****************************************************
	 * ************************************************ *
	 * *        Newton Polynomial Definitions         * *
	 * ************************************************ *
	 ****************************************************/

	template <int N, typename Type>
	Newton_polynomial<N, Type>::Newton_polynomial():
	count( 0 ),
	begin( 0 ),
	end( N - 1 ) {
		// Empty constructor
	}

	/*
	 * Evaluate the Newton Polynomial at x
	 *
	 *   BACKWARD    Evaluate so as to most efficiently evaluate the
	 *               Newton polynomial at the point near x[1]
	 *   FORWARD     Evaluate so as to most efficiently evaluate the
	 *               Newton polynomial at the point near x[N]
	 *   OPTIMAL     Evaluate so as to most efficiently evaluate the
	 *               Newton polynomial at any point for the given
	 *               divided-difference table.
	 *
	 *  This function uses Horner's rule.
	 *  Run time:  O(N)
	 *  Memory:    O(1)
	 */

	template <int N, typename Type>
	Type Newton_polynomial<N, Type>::evaluate( Type const &x, direction dir ) const {
		// Deal with the 0 function and the constant function

		if ( count == 0 ) {
			return Type();
		} else if ( count == 1 ) {
			return table[end][0];
		}

		Type y;
		int b, e;
		bool forward;

		switch ( dir ) {
			case BACKWARD:
				// Evaluate the upper edge of the divided-difference
				// table as follows:
				//
				//      y1
				//         y21
				//      *      y321
				//          *     y4321
				//      *       *
				//          *
				//      *
				//
				//   y1 + y21*(x - x1) + y321*(x - x1)*(x - x2)
				//      + y4321*(x - x3)*(x - x2)*(x - x1)

				y = table[end][count - 1];
				b = next( begin );

				for ( int k = count - 2; k >= 0; --k ) {
					y = y*(x - x_values[b]) + table[end][k];
					b = next( b );
				}

				return y;
			case FORWARD:
				// Evaluate the lower edge of the divided-difference
				// table as follows:
				//
				//      *
				//          *
				//      *       *
				//          *     y4321
				//      *      y432
				//         y43
				//      y4
				//
				//   y4 + y43*(x - x4) + y432*(x - x4)*(x - x3)
				//      + y4321*(x - x4)*(x - x3)*(x - x2)

				y = table[end][count - 1];
				e = previous( end );

				for ( int k = count - 2; k >= 0; --k ) {
					y = y*(x - x_values[e]) + table[e][k];
					e = previous( e );
				}

				return y;
			case OPTIMAL:
				// Starting with the tip of the divided-difference table, y4321,
				// move away from the x[k]-value further from the point which is
				// being evaluated.
				//
				//   x[begin]
				//            *
				//      *           *
				//            *         y4321
				//      *           *
				//            *
				//    x[end]

				y = table[end][count - 1];
				b = begin;
				e = end;

				forward = std::fabs( x - x_values[b] ) > std::fabs( x - x_values[e] );

				for ( int k = count - 2; k >= 0; --k ) {
					if ( forward ) {
						b = next( b );
					} else {
						e = previous( e );
					}

					forward = std::fabs( x - x_values[b] ) > std::fabs( x - x_values[e] );

					if ( forward ) {
						y = y*(x - x_values[b]) + table[e][k];
					} else {
						y = y*(x - x_values[e]) + table[e][k];
					}
				}

				assert( b == e );

				return y;
			default:
				assert( false );
		}
	}

	/*
	 * Insert Another Point to Interpolate
	 *
	 *  Update the divided difference table.
	 *  Run time:  O(N)
	 *  Memory:    O(1)
	 */

	template <int N, typename Type>
	void Newton_polynomial<N, Type>::insert( Type const &x, Type const &y ) {
		end = next( end );

		x_values[end] = x;
		table[end][0] = y;

		if ( count == N ) {
			begin = next( begin );
		} else {
			++count;
		}

		for ( int m = 1; m < count; ++m ) {
			int ix =  begin;
			int jx = remainder( begin + m );

			int iy = remainder( begin + m - 1 );
			int jy = remainder( begin + m );

			for ( int n = m; n < count; ++n ) { 
				table[jy][m] = (table[jy][m - 1] - table[iy][m - 1])/(x_values[jx] - x_values[ix]); 

				ix = next( ix );
				jx = next( jx );

				iy = next( iy );
				jy = next( jy );
			}
		}
	}

	/*
	 *  Clear the Divided-Difference Table
	 *   void Newton_polynomial<N, Type> :: clear() 
	 *
	 *  Empty the divided difference table.
	 *  O(1)
	 */

	template <int N, typename Type>
	void Newton_polynomial<N, Type>::clear() {
		count = 0;
		begin = 0;
		end = N - 1;
	}

	template <int N, typename Type>
	Type Newton_polynomial<N, Type>::x( int k ) const {
		assert ( 0 <= k && k < N );

		return x_values[begin + k >= N ? begin + k - N : begin + k];
	}

	template <int N, typename Type>
	Type Newton_polynomial<N, Type>::y( int j, int k ) const {
		if ( k < j ) {
			return y( k, j );
		}

		assert ( 0 <= j && k < N );

		int column = k - j;

		return table[begin + k >= N ? begin + k - N : begin + k][column];
	}

	/*
	 *  The Next Index in the Divided Difference Table
	 *   int Newton_polynomial<N, Type> :: next( n ) 
	 *
	 *  Increment n and return 0 if the result would be N.
	 *  O(1)
	 */

	template <int N, typename Type>
	inline int Newton_polynomial<N, Type>::next( int n ) {
		return (n == N - 1) ? 0 : n + 1;
	}

	/*
	 *  The Previous Index in the Divided Difference Table
	 *   int Newton_polynomial<N, Type> :: previous( n ) 
	 *
	 *  Decrement n and return N - 1 if the result would be -1.
	 *  O(1)
	 */

	template <int N, typename Type>
	inline int Newton_polynomial<N, Type>::previous( int n ) {
		return (n == 0) ? N - 1 : n - 1;
	}

	/*
	 *  The Remainder of n/N for n in [0, 2*N)
	 *   int Newton_polynomial<N, Type> :: previous( n ) 
	 *
	 *  Return the remainder of n/N for values of n in 0 <= n < 2*N.
	 *  O(1)
	 */

	template <int N, typename Type>
	inline int Newton_polynomial<N, Type>::remainder( int n ) {
		return (n >= N) ? n - N : n;
	}

	/****************************************************
	 * ************************************************ *
	 * *                   Friends                    * *
	 * ************************************************ *
	 ****************************************************/


	/*********************************************************************
	 * out << poly;
	 *
	 *  Print the Newton polynomial in the classical format:
	 *    y1 + y21*(x - x[1]) + y321*(x - x[1])*(x - x[2]) + y4321*(x - x[3])*(x - x[2])*(x - x[1])
	 *********************************************************************/

	template <int N, typename Type>
	std::ostream &operator<<( std::ostream &out, const Newton_polynomial<N, Type> &poly ) {
		if ( poly.count == 0 ) {
			out << 0;
			return out;
		}

		out << poly.table[poly.begin][0];

		int m = ( poly.begin == N - 1 ) ? 0 : poly.begin + 1;

		for ( int i = 1; i < poly.count; ++i ) {
			if ( poly.table[m][i] < 0 ) {
				out << " - " << (-poly.table[m][i]);
			} else if ( poly.table[m][i] == 0 ) {
				continue;
			} else {
				out << " + " << poly.table[m][i];
			}

			int n = poly.begin;

			for ( int j = 0; j < i; ++j ) {
				if ( poly.x_values[n] < 0 ) {
					out << "*(x + " << (-poly.x_values[n]) << ")";
				} else if ( poly.x_values[n] == 0 ) {
					out << "*x";
				} else {
					out << "*(x - " << poly.x_values[n] << ")";
				}

				n = (n == N - 1) ? 0 : n + 1;
			}

			m = (m == N - 1) ? 0 : m + 1;
		}

		return out;
	}
}

#endif