/*
 *  Spline.c: Spline Interpolation スプライン補間
 */
#include <stdio.h>
#include <math.h>
#define ND  50
#define MD  20

int Hokan(double xp[], double yp[], double u[], 
double A2[], double A3[], double y[], int np, int mm)
/*
    スプライン補間
 */
{
    int i, k;
    double xmin, xmax, dx, x, t;

    xmin = xp[0];
    xmax = xp[np-1];
    dx = (xmax-xmin)/mm;
    k = 1;
    printf("      x            y\n");
    for (i=0; i<=mm; i++) {
        x = xmin+dx*i;
        if (x>=xp[k]) k = k+1;
        t = x-xp[k-1];
        y[i] = ((A3[k]*t+A2[k])*t+u[k-1])*t+yp[k-1];
        printf(" %10.6f %e\n",x,y[i]);
    }
    return 0;
}

int Graph(double y[], int mm)
/*
    グラフを描く
 */
{
    int i, j, iy;
    double ymin, ymax, sf;

    ymin = 999.99e30;
    ymax = -ymin;
    for (i=0; i<=mm; i++) {
      if (y[i]<ymin) ymin = y[i];
      if (y[i]>ymax) ymax = y[i];
    }
    sf = 60.0/(ymax-ymin);
    putchar('\n');
/*  プロット */
    for (i=0; i<=mm; i++) {
        iy = (int)((y[i]-ymin)*sf+0.5);
        for (j=1; j<=iy; j++) putchar(' ');
        printf("*\n");
    }
    return 0;
}

int Spline(double xp[], double yp[], int np, double s[],
double u[], double A2[], double A3[])
/*
    スプライン係数の計算
 */
{
    int k, n;
    double A[ND+1], B[ND+1], C[ND+1], D[ND+1];
    double hl, hr, C1, C2, h1;

    n = np-2;
/*  連立方程式の係数の計算 */
    for (k=1; k<=n; k++) {
        hl = xp[k]-xp[k-1];
        hr = xp[k+1]-xp[k];
        A[k] = hl/(hl+hr);
        B[k] = 1-A[k];
        C1 = (yp[k+1]-yp[k])/hr;
        C2 = (yp[k]-yp[k-1])/hl;
        C[k] = 6*((C1-C2)/(hl+hr));
        D[k] = 2;
    }
/*  連立方程式を解く(前進消去) */
    for (k=1; k<=n-1; k++) {
        B[k] = B[k]/D[k];
        C[k] = C[k]/D[k];
        D[k+1] = D[k+1]-A[k+1]*B[k];
        C[k+1] = C[k+1]-A[k+1]*C[k];
    }
/*  (後退代入) */
    s[n] = C[n]/D[n];
    for (k=n-1; k>=1; k--) s[k] = C[k]-B[k]*s[k+1];
/*  微分係数u[k],A2[k],A3[k]の計算 */
    h1 = xp[1]-xp[0];
    u[0] = (yp[1]-yp[0])/h1-h1*s[1]/6;
    for (k=1; k<=n+1; k++) {
        hl = xp[k]-xp[k-1];
        u[k] = u[k-1]+hl*(s[k-1]+s[k])/2;
        A3[k] = (s[k]-s[k-1])/(6*hl);
        A2[k] = s[k-1]/2;
    }
    return 0;
}
/*    ------------------    Main    ------------------    */
int main()
{
    int i, np, mm;
    double xp[ND+1],yp[ND+1],s[ND+1],u[ND+1],a2[ND+1],a3[ND+1];
    double y[MD+1];
    char c[135];

    mm = MD;
/*  入力 */
    printf("点の数を入れてください　"); gets(c); sscanf(c,"%d",&np);
    for (i=0; i<=np-1; i++) {
        printf(" x[%2d], y[%2d] = ",i,i);
        gets(c); sscanf(c,"%lf%lf", &xp[i], &yp[i]);
    }
/*  計算 */
    Spline(xp, yp, np, s, u, a2, a3);
/*  出力 */
    Hokan(xp, yp, u, a2, a3, y, np, mm);
    Graph(y, mm);
}
/* -- 実行例:戸川隼人「数値計算」岩波書店 p.160 プログラム7.4 -- */
/*
Z:\src\C>Spline
点の数を入れてください　5
 x[ 0], y[ 0] = 1 1
 x[ 1], y[ 1] = 2 1
 x[ 2], y[ 2] = 3 3
 x[ 3], y[ 3] = 4 1
 x[ 4], y[ 4] = 5 1
      x            y
   1.000000 1.000000e+00
   1.200000 8.354286e-01
   1.400000 7.120000e-01
   1.600000 6.708571e-01
   1.800000 7.531429e-01
   2.000000 1.000000e+00
   2.200000 1.427429e+00
   2.400000 1.950857e+00
   2.600000 2.460571e+00
   2.800000 2.846857e+00
   3.000000 3.000000e+00
   3.200000 2.846857e+00
   3.400000 2.460571e+00
   3.600000 1.950857e+00
   3.800000 1.427429e+00
   4.000000 1.000000e+00
   4.200000 7.531429e-01
   4.400000 6.708571e-01
   4.600000 7.120000e-01
   4.800000 8.354286e-01
   5.000000 1.000000e+00

        *
    *
 *
*
  *
        *
                   *
                                 *
                                              *
                                                        *
                                                            *
                                                        *
                                              *
                                 *
                   *
        *
  *
*
 *
    *
        *

*/