// x:存放給定n個節點的值，且x0<x1<...<xn-1
// y:存放n個給定節點上的函數值y0,y1,...,yn-1
// 給定節點的個數
// 指定插值點
// 函數反回指定插值點t處的函數近似值z=f(t)

  double lgr(x,y,n,t)
  int n;
  double t,x[],y[];
  { int i,j,k,m;
    double z,s;
    z=0.0;
    if (n<1) return(z);
    if (n==1) { z=y[0];return(z);}
    if (n==2)
      { z=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
        return(z);
      }
    i=0;
    while ((x[i]<t)&&(i<n)) i=i+1;
    k=i-4;
    if (k<0) k=0;
    m=i+3;
    if (m>n-1) m=n-1;
    for (i=k;i<=m;i++)
      { s=1.0;
        for (j=k;j<=m;j++)
          if (j!=i) s=s*(t-x[j])/(x[i]-x[j]);
        z=z+s*y[i];
      }
    return(z);
  }