// LagranreInterpolate.cpp: implementation of the LagranreInterpolate class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "LagranreInterpolate.h"
#include <assert.h>

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

LagranreInterpolate::LagranreInterpolate(int n, double x[], double y[])
{
	assert(n>0);
	m_pdNodeX = new double[n+1];
	m_pdNodeY = new double[n+1];
	m_pdBaseDeno = new double[n+1];
	m_nNum = n;
	for(int i =0; i<=n; i++)
	{
		m_pdNodeX[i] = x[i];
		m_pdNodeY[i] = y[i];
	}

	CalculateBaseDeno();
}

LagranreInterpolate::~LagranreInterpolate()
{
	delete[]  m_pdNodeX;
	delete[]  m_pdNodeY;
}

double LagranreInterpolate::Evaluate(double x)
{
	double omega = 1;
	int i;
	for(i=0; i<=m_nNum; i++)
		omega *= x-m_pdNodeX[i];
	
	double value = 0;
	for(i=0; i<=m_nNum; i++)
		value += m_pdNodeY[i]/((x-m_pdNodeX[i])*m_pdBaseDeno[i]);
	value *= omega;

	return value;
}

void LagranreInterpolate::CalculateBaseDeno()
{
	//初始化分母為1
	int i, j;
	for(i=0; i<=m_nNum; i++)
		m_pdBaseDeno[i] = 1;
	
	//xi-xj會出現在兩個分母中，把它乘入
	int sign = 1;
	for(i=0; i<=m_nNum; i++)
	{
		for(j=i+1; j<=m_nNum; j++)
		{
			m_pdBaseDeno[i] *= m_pdNodeX[i]-m_pdNodeX[j];
			m_pdBaseDeno[j] *= m_pdNodeX[i]-m_pdNodeX[j];
		}
		m_pdBaseDeno[i] *= sign;
		sign *= -1;
	}
}