#include <stdio.h>
#include <math.h>

// FALTA: 'SPLINE (PRECISA SPLINE?)

/*
Funções básicas de interpolação, do tipo linear,
bilinear, polinomial e spline.
Para as interpolações lineares, polinomiais e spline,
o funcionamento é:

	y0 = interpol_xxx(n, x, y, x0)

onde x e y são vetores de double de tamanho n, que 
descrevem em apenas alguns pontos uma função y = f(x).
O x0 é um double que é o ponto que queremos obter a
interpolação em questão.

A interpolação bilinear possui o seguinte funcionamento:

	z0 = interpol_bilinear(n, m, x, y, z, x0, y0)

onde x é um vetor de tamanho n e y é um vetor de tamanho m,
z é uma matriz n por m que descreve uma função z = f(x,y)
em alguns pontos. Importante notar que o algoritmo realiza
a interpolação primeiro no eixo y e depois no x.

Exemplo:

	Função y(x) = x^2 + x - 1.
	Definimos x = {1,2,3,4,5} e n = 5.
	Logo obtemos y = {1,5,11,19,29}, aplicando y(x) em um
	loop.
	-- A partir de agora, só temos conhecimento de x e y,
	-- a expressão de y(x) não é mais conhecida.
	Queremos y0 = y(x0) em x0 = 3.5. Claramente não temos
	o valor exato nos vetores, logo necessitamos realizar
	uma interpolação:
	x0 = 3.5
	y0 = interpol_xxx(n, x, y, x0)
	No final temos que y0 = 15, muito próximo do resultado
	real, y(3.5) = 14.75

Observações:

Um polinômio de grau N-1 pode ser descrito perfeitamente por N pontos, e vice-versa.

Para muitos pontos (n > 5), a interpolação polinomial se
torna instável, então não é recomendada.
*/

double interpol_linear(int n, double x[n], double y[n], double x0)
{
	int k, i;
	for(i = 1; i < n; i++) // Encontrar índices vizinhos
	{
		if(x0 <= x[i])
		{
			k = i - 1;
			break;
		}
	}
	if(i == n) k = i - 2;

	return (y[k+1] - y[k])/(x[k+1] - x[k])*(x0 - x[k]) + y[k];
}

double interpol_polinomial(int n, double x[n], double y[n], double x0)
{
	double soma = 0;
	int i, j;

	for(i = 0; i < n; i++)
	{
		double produto = 1;
		for(j = 0; j < n; j++)
		{
			if(j != i)
			{
				produto *= (x0 - x[j])/(x[i] - x[j]);
			}
		}
		soma += produto * y[i];
	}

	return soma;
}

double interpol_bilinear(int n, int m, double x[n], double y[m], double z[n][m], double x0, double y0)
{
	double zn1, zn2;

	int i,kx;
	for(i = 1; i < n; i++) // Encontrar índices vizinhos
	{
		if(x0 <= x[i])
		{
			kx = i - 1;
			break;
		}
	}

	zn1 = interpol_linear(m, y, z[kx], y0);
	zn2 = interpol_linear(m, y, z[kx+1], y0);

	double xpos[2] = {x[kx],x[kx+1]};
	double zpos[2] = {zn1,zn2};
	return interpol_linear(2, xpos, zpos,x0);
}

/* ----------- COMEÇO DO EXEMPLO ----------- */

double f(double x)
{
	return 0.1*x*x*x + x*x + x - 1;
}

int main(void)
{
	double x[4] = {1,2,3,4}, y[4], x0 = 0.9;

	int i;
	for(i = 0; i < 4; i++) y[i] = f(x[i]);

	printf("x\t=\t%f\nf(x)\t=\t%f\nl(x)\t=\t%f\np(x)\t=\t%f\n", x0, f(x0), interpol_linear(4,x,y,x0), interpol_polinomial(4,x,y,x0));

	double zn[2][2] ={{0,1},{2,3}};
	double xn[2] = {0,1};
	double yn[2] = {0,1};
	printf("%f",interpol_bilinear(2,2,xn,yn,zn,0.5,0.5));

	// Interpolação polinomial retorna ponto exato pois a
	// função é de 3° grau, e são 4 pontos.
}


/* ------------ FIM DO EXEMPLO ------------ */