36 lines
710 B
C
36 lines
710 B
C
#include "solve.h"
|
|
|
|
#include <float.h>
|
|
#include <math.h>
|
|
#include <stdio.h>
|
|
|
|
|
|
// Interpolation polynomial according to Aitken's scheme
|
|
double t3_solve (const double x_0, const int n, const double * restrict X, double * restrict Y)
|
|
{
|
|
for (int k = 0; k < n-1; ++k)
|
|
{
|
|
double x_j;
|
|
double y_j = Y[n-1];
|
|
|
|
// printf ("------- K = %d -------\n", k);
|
|
|
|
for (int i = n-2; i >= k; --i)
|
|
{
|
|
const double x_i = X[i-k];
|
|
const double y_i = Y[i];
|
|
x_j = X[i+1];
|
|
|
|
if (fabs(x_j - x_i) < DBL_EPSILON)
|
|
return DBL_MAX;
|
|
|
|
Y[i+1] = ((x_0-x_i) * y_j - (x_0-x_j) * y_i) / (x_j - x_i);
|
|
// printf ("I = %d, f(x%d, ... , x%d) = %lf\n", i, i-k+1, i+2, Y[i+1]);
|
|
|
|
y_j = y_i;
|
|
}
|
|
}
|
|
|
|
return Y[n-1];
|
|
}
|
|
|