/*
 *  Nwtint.c: Newton Interpolation ニュートン補間
 */
#include <stdio.h>
#include <math.h>
#define ND 100
#define MD  10

int Kaisa(double yy[], double deltas[][ND+1], int m, int n)
/*
    階差表を作る
 */
{
    int j, k;

    for (k=1; k<=n; k++) deltas[0][k] = yy[k];
    for (j=0; j<=m-1; j++)
        for (k=1; k<=n-j-1; k++)
            deltas[j+1][k] = deltas[j][k+1]-deltas[j][k];
    printf("階差テーブル\n");
    for (k=1; k<=n; k++) {
        for (j=1; j<=m; j++)
          if (k<=(n-j)) printf("%11.7f",deltas[j][k]);
        putchar('\n');
    }
    return 0;
}

double Nwtint(double x, double deltas[][ND+1], double xx[],
int k, int m)
/*
    ニュートン補間
 */
{
    int i;
    double dx, s, f, g, y;

    dx = xx[k+1]-xx[k];
    s = (x-xx[k])/dx;
    f = 1.0; g = s;
    y = deltas[0][k];
    for (i=1; i<=m; i++) {
        f *= (double)i;
        y += g*deltas[i][k]/f;
        g *= s-(double)i;
    }
    return y;
}

int Suuhyou(double xx[], double yy[])
/*
    数表の擬似データ作成
    y = sin(x) (x=10,...,180)
 */
{
    int i, n;
    double t;

    n = 18;
    printf("   i     xx[i]     yy[i]\n");
    for (i=1; i<=n; i++) {
        xx[i] = 10.0*(double)i;
        t = 3.14159265358979*xx[i]/180.0;
        yy[i] = sin(t);
        printf("%4d%10.3f%10.6f\n",i,xx[i],yy[i]);
    }
    return n;
}

/*    ------------------    Main    ------------------    */
int main()
{
    int i, k, m, n;
    double x, y, xx[ND+1], yy[ND+1], deltas[MD+1][ND+1];
    char s[135];
/*  数表を読み込む */
    n = Suuhyou(xx, yy);
/*  階差表を作る */
    m = 7;
    Kaisa(yy, deltas, m, n);
/*  xを入れてp(x)を計算してみる */
    for (i=1; i<=3; i++) {
        printf(" x,k,m= "); gets(s); sscanf(s,"%lf%d%d",&x,&k,&m);
        y = Nwtint(x, deltas, xx, k, m);
        printf(" y= %e\n", y);
    }
    return 0;
}
/* -- 実行例:戸川隼人「数値計算」岩波書店 p.155 プログラム7.3 -- */
/*
Z:\src\C>Nwtint
   i     xx[i]     yy[i]
   1    10.000  0.173648
   2    20.000  0.342020
   3    30.000  0.500000
   4    40.000  0.642788
   5    50.000  0.766044
   6    60.000  0.866025
   7    70.000  0.939693
   8    80.000  0.984808
   9    90.000  1.000000
  10   100.000  0.984808
  11   110.000  0.939693
  12   120.000  0.866025
  13   130.000  0.766044
  14   140.000  0.642788
  15   150.000  0.500000
  16   160.000  0.342020
  17   170.000  0.173648
  18   180.000  0.000000
階差テーブル
  0.1683720 -0.0103921 -0.0048001  0.0004616  0.0001318 -0.0000180 -0.0000035
  0.1579799 -0.0151922 -0.0043385  0.0005934  0.0001138 -0.0000215 -0.0000028
  0.1427876 -0.0195308 -0.0037451  0.0007072  0.0000923 -0.0000243 -0.0000021
  0.1232568 -0.0232759 -0.0030379  0.0007995  0.0000680 -0.0000264 -0.0000013
  0.0999810 -0.0263137 -0.0022383  0.0008675  0.0000417 -0.0000276 -0.0000004
  0.0736672 -0.0285521 -0.0013708  0.0009092  0.0000140 -0.0000281  0.0000004
  0.0451151 -0.0299229 -0.0004616  0.0009232 -0.0000140 -0.0000276  0.0000013
  0.0151922 -0.0303845  0.0004616  0.0009092 -0.0000417 -0.0000264  0.0000021
 -0.0151922 -0.0299229  0.0013708  0.0008675 -0.0000680 -0.0000243  0.0000028
 -0.0451151 -0.0285521  0.0022383  0.0007995 -0.0000923 -0.0000215  0.0000035
 -0.0736672 -0.0263137  0.0030379  0.0007072 -0.0001138 -0.0000180  0.0000040
 -0.0999810 -0.0232759  0.0037451  0.0005934 -0.0001318 -0.0000140
 -0.1232568 -0.0195308  0.0043385  0.0004616 -0.0001458
 -0.1427876 -0.0151922  0.0048001  0.0003158
 -0.1579799 -0.0103921  0.0051159
 -0.1683720 -0.0052762
 -0.1736482

 x,k,m= 45 3 5
 y= 7.071069e-01
 x,k,m= 45 3 2
 y= 7.068574e-01
 x,k,m= 45 2 3
 y= 7.071285e-01
*/