#pragma once

#include "./numerics/interpolation1d/interpolation1d_base.h"

#include "./numerics/abs.h"
#include "./utils/vector.h"


namespace numerics{
	
	template <typename T>
	struct interp_polynomial : Base_interp<T> {
		using Base   = Base_interp<T>;
	  	// bring base data members into scope (or use this->xx / this->yy below)
	  	using Base::xx;
	  	using Base::yy;
	  	using Base::mm;

	  	T dy;


		interp_polynomial(const utils::Vector<T> &xv, const utils::Vector<T> &yv, uint64_t m)
		: Base_interp<T>(xv, &yv[0], m), dy(T{0}){}

		T rawinterp(int64_t jl, T x) override{

			int64_t ns=0;
			T y, den, dif, dift, ho, hp, w;
			const T *xa = &xx[jl], *ya = &yy[jl];
			utils::Vector<T> c(mm,0), d(mm,0);
			dif = numerics::abs(x-xa[0]);

			for (int64_t i = 0; i < mm; ++i){					// Here we find the index ns of the closest table entry,
				dift = numerics::abs(x-xa[i]);
				if (dift < dif){
					ns = i;
					dif=dift;
				}
				c[i]=ya[i];									// and initialize the tableau of c’s and d’s.
				d[i]=ya[i];
			}
			y = ya[ns];										// This is the initial approximation to y.
			ns -= 1;

			for (int64_t m = 1; m < mm; ++m){				// For each column of the tableau,
				for (int64_t i = 0; i < mm-m; ++i){		// we loop over the current c’s and d’s and update them.
					ho = xa[i]-x;
					hp = xa[i+m]-x;
					w = c[i+1]-d[i];
					den = ho-hp;
					if (den == T{0.0}){
						throw std::invalid_argument("interp_polynomial error"); // This error can occur only if two input xa’s are (to within roundoff identical.
					}
					den = w/den;							// Here the c’s and d’s are updated.
					d[i] = hp*den;
					c[i] = ho*den;
				}
				bool take_left = 2 * (ns + 1) < (mm - m);

				if (take_left) {
				    dy = c[ns + 1];
				    y += dy;
				} else {
				    dy = d[ns];   
				    y += dy;
				    ns -= 1;
				}
			}
			return y;
		}

	};

} // namespace numerics