/*
 *  Lagint.c: Lagrange Interpolation ラグランジュ補間
 */
#include <stdio.h>
#include <math.h>
#define ND 100

double Nix(int i, int k, int m, double x, double xx[])
/*
    Ni(x)の計算
 */
{
    int j;
    double p;

    p = 1;
    for (j=k; j<=k+m; j++)
        if (i!=j) p *= (x-xx[j])/(xx[i]-xx[j]);
    return p;
}
double Lagint(double x, double xx[], double yy[], int k, int m)
/*
    ラグランジュ補間
 */
{
    int i;
    double s;

    s = 0;
    for (i=k; i<=k+m; i++)
        s += yy[i]*Nix(i, k, m, x, xx);
    return s;
}
int Datain(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*i;
        t = 3.14159265*xx[i]/180;
        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,yt, xx[ND+1], yy[ND+1];
    char s[135];
/* 数表を読み込む */
    n = Datain(xx, yy);
/* xを入れてp(x)を計算してみる */
    for (i=1; i<=4; i++) {
        printf(" x,k,m= "); gets(s); sscanf(s,"%lf%d%d",&x,&k,&m);
        y = Lagint(x, xx, yy, k, m);
        yt= sin(3.14159265*x/180);
        printf(" y= %e, y(true)=%e\n", y,yt);
    }
    return 0;
}
/* -- 実行例:戸川隼人「数値計算」岩波書店 p.152 プログラム7.2 -- */
/*
Z:\src\C>Lagint
   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
 x,k,m= 45 3 5
 y= 7.071069e-01, y(true)=7.071068e-01
 x,k,m= 45 1 2
 y= 7.174846e-01, y(true)=7.071068e-01
 x,k,m= 45 3 2
 y= 7.174846e-01, y(true)=7.071068e-01
 x,k,m= 45 2 3
 y= 7.071285e-01, y(true)=7.071068e-01
 y= 7.071285e-01
*/