/*    測試數據來自課本55頁例5        《數值分析》清華大學出版社第四版*/
//輸入
3
27.7    4.1
28    4.3
29    4.1
30    3.0
1
3.0    -4.0
//輸出
輸出三次樣條插值函數：
1: [27.7 , 28]
        13.07*(x - 28)^3 + 0.22*(x - 27.7)^3
        + 14.84*(28 - x) + 14.31*(x - 27.7)

2: [28 , 29]
        0.066*(29 - x)^3 + 0.1383*(x - 28)^3
        + 4.234*(29 - x) + 3.962*(x - 28)

3: [29 , 30]
        0.1383*(30 - x)^3 - 1.519*(x - 29)^3
        + 3.962*(30 - x) + 4.519*(x - 29)


//三次樣條插值函數
#include<iostream>
#include<iomanip>
using namespace std;

const int MAX = 50;
float x[MAX], y[MAX], h[MAX];
float c[MAX], a[MAX], fxym[MAX];

float f(int x1, int x2, int x3){
    float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);
    float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);
    return (a - b)/(x[x3] - x[x1]);
}    //求差分
void cal_m(int n){        //用追趕法求解出彎矩向量M……
    float B[MAX];
    B[0] = c[0] / 2;
    for(int i = 1; i < n; i++)
            B[i] = c[i] / (2 - a[i]*B[i-1]);

    fxym[0] = fxym[0] / 2;
    for(i = 1; i <= n; i++)
        fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);
    for(i = n-1; i >= 0; i--)
        fxym[i] = fxym[i] - B[i]*fxym[i+1];
}
void printout(int n);

int main(){
    int n,i; char ch;
    do{
        cout<<"Please put in the number of the dots:";
        cin>>n;
        for(i = 0; i <= n; i++){
            cout<<"Please put in X"<<i<<':';
            cin>>x[i];    //cout<<endl;
            cout<<"Please put in Y"<<i<<':';
            cin>>y[i]; //cout<<endl;
        }
        for(i = 0; i < n; i++)            //求 步長
            h[i] = x[i+1] - x[i];
        cout<<"Please 輸入邊界條件\n 1: 已知兩端的一階導數\n 2:兩端的二階導數已知\n 默認:自然邊界條件\n";
        int t;
        float f0, f1;
        cin>>t;
        switch(t){
        case 1:cout<<"Please put in Y0\' Y"<<n<<"\'\n";
            cin>>f0>>f1;
            c[0] = 1; a[n] = 1;
            fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];
            fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];
            break;
        case 2:cout<<"Please put in Y0\" Y"<<n<<"\"\n";
            cin>>f0>>f1;
            c[0] = a[n] = 0;
            fxym[0] = 2*f0; fxym[n] = 2*f1;
            break;
        default:cout<<"不可用\n";//待定
        };//switch
        for(i = 1; i < n; i++)
            fxym[i] = 6 * f(i-1, i, i+1);
        for(i = 1; i < n; i++){
            a[i] = h[i-1] / (h[i] + h[i-1]);
            c[i] = 1 - a[i];
        }
        a[n] = h[n-1] / (h[n-1] + h[n]);
        cal_m(n);
        cout<<"\n輸出三次樣條插值函數：\n";
        printout(n);
        
        cout<<"Do you to have anther try ? y/n :";
        cin>>ch;
    }while(ch == 'y' || ch == 'Y');
    return 0;
}
void printout(int n){
    cout<<setprecision(6);
    for(int i = 0; i < n; i++){
        cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";
        /*
        cout<<fxym[i]/(6*h[i])<<" * ("<<x[i+1]<<" - x)^3 + "<<<<" * (x - "<<x[i]<<")^3 + "
            <<(y[i] - fxym[i]*h[i]*h[i]/6)/h[i]<<" * ("<<x[i+1]<<" - x) + "
            <<(y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i]<<"(x - "<<x[i]<<")\n";
        cout<<endl;*/
        float t = fxym[i]/(6*h[i]);
        if(t > 0)cout<<t<<"*("<<x[i+1]<<" - x)^3";
        else cout<<-t<<"*(x - "<<x[i+1]<<")^3";
        t = fxym[i+1]/(6*h[i]);
        if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")^3";
        else cout<<" - "<<-t<<"*(x - "<<x[i]<<")^3";
        cout<<"\n\t";
        t = (y[i] - fxym[i]*h[i]*h[i]/6)/h[i];
        if(t > 0)cout<<"+ "<<t<<"*("<<x[i+1]<<" - x)";
        else cout<<"- "<<-t<<"*("<<x[i+1]<<" - x)";
        t = (y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i];
        if(t > 0)cout<<" + "<<t<<"*(x - "<<x[i]<<")";
        else cout<<" - "<<-t<<"*(x - "<<x[i]<<")";
        cout<<endl<<endl;
    }
    cout<<endl;
}