外插法

#include <math.h>
   
typedef double S8B;

//輸入: X, Y, n(元素數目), x, 輸出: y, dy(精確度)

void extrapolate (S8B X[], S8B Y[], int n, S8B x, S8B& y, S8B& dy)
{
    S8B  t, d; 
    S8B  dx = fabs (x - X[0]);
    S8B* C  = new S8B[n];
    S8B* D  = new S8B[n];
    int  i, m, ns = 0; 
    dy = y = 0;
   
    for (i=0; i < n; i++)
    {
        d = fabs (x - X[i]);
        if (d == 0) {y = Y[i]; dy=0; return;}
        else if (d < dx) {ns = i; dx = d;}
        C[i] = Y[i];
        D[i] = Y[i] + 1e-25;   //防止下面發生 0/0
    }
    y = Y [ns--];

    for (m=1; m<n; m++) 
    {
        for (i=0; i<n-m; i++)
        {
            t = (X[i]-x) * D[i] / (X[i+m]-x);
            d = t - C[i+1];
            if (d == 0) cout <<"x 處存在極點";
            d = (C[i+1] - D[i]) / d;
            C[i] = d*t;
            D[i] = C[i+1]*d;
        }
        y += (dy = ns+ns+2< n-m? C[ns+1]: D[ns--]);        
    }
    delete[] C;
    delete[] D;
}

int main()
{
    S8B X[] = {1920, 1930, 1940, 1950, 
               1960, 1970, 1980, 1990};
    S8B Y[] = {106.46, 123.08, 123.12, 152.27, 
               180.67, 205.05, 227.23, 249.46};
    S8B y, dy;
    extrapolate (X, Y, 8, 2000, y, dy);
    printf ("y = %lf ± %lf", y, dy);
}