#include<stdio.h>
#include<math.h>
#include<stdlib.h>

/*********** 表2.1のスプライン補間の誤差を出力するプログラム **************

	・実行すると，スプライン補間の節点の数N+1に対応する
		N(textのNに対応)の入力を求めます．

	・コマンドラインの第一引数にNを与えることでも
		即時に誤差を出力することができます．
		例）	>spline 12

	・L_2ノルムの2乗の計算には積分の台形公式(式(3.4))を用います．

**************************************************************************/

#define nn 100			//小区間中のきざみ数(L2nの精度に関わる)

double f(double x){
	double y;
	y=1./(1.+25.*pow(x,2));
	return y;
}

int main(int argc, char **argv){
	int i, j, jmax, N;
	double h, dx, temp, x, p, L2n;
	double *a, *b, *c;
	double *xi, *fi, *Y, *C;

		/********* 入力処理(始) ***********/

	if(!(argc==1 || argc==2)){
		printf("引数のエラーです\n");
		return 1;
	}
	if(argc == 1){
		printf("スプライン補間の節点の数N+1に対応する						Nを入力してください．\n N = ");
		scanf("%d",&N);
	}else{
		argv++;
		N = atoi(*argv);
	}

		/********* 入力処理(終) **********/

	/*** 右側コメント文と同等のメモリ領域を確保 ***/
	/*** a:α  b:β  c:γ ，求める連立一次方程式のXiはC[2][] ***/
	a  = malloc(sizeof(double)*(N+2));		//a[N+2]
	b  = malloc(sizeof(double)*(N+2));		//b[N+2]
	c  = malloc(sizeof(double)*(N+1));		//c[N+1]
	xi = malloc(sizeof(double)*(N+2));		//xi[N+2]
	fi = malloc(sizeof(double)*(N+2));		//fi[N+2]
	Y  = malloc(sizeof(double)*(N+2));		//Y[N+2]
	C  = malloc(sizeof(double)*5*(N+2));		//C[5][N+2]

	if(xi==NULL || fi==NULL || a==NULL || b==NULL || c==NULL || 				Y==NULL || C==NULL){
		printf("can't malloc!!\n");
		return 2;
	}

	dx = 2./N;	//小区間の幅(節点間の幅) dx
	h = dx/nn;	//小区間のきざみ幅 h(台形公式のh)

	for(i=0; i<N+2; i++){
		xi[i] = -1.+dx*(double)i;
		fi[i] = f(xi[i]);
	}

	for(i=1; i<N+1; i++){
//		C[1][i] = f(xi[i-1]);
		C[1*(N+2)+i] = f(xi[i-1]);
	}

	a[1] = 2.*dx;
	a[N+1] = 2.*dx;
	c[1] = dx;
	Y[1] = 3.*(fi[1]-fi[0]);
	Y[N+1] = 3.*(fi[N]-fi[N-1]);

	for(i=2; i<N+1; i++){
		a[i] = 4.*dx;
		c[i] = dx;
		Y[i] = 3.*(fi[i]-fi[i-2]);
	}

	for(i=2; i<N+2; i++)
		b[i] = dx;

	/***** 行列の係数を置き換えて上三角行列に直す *****/
	for(i=2; i<N+2; i++){
		temp = b[i]/a[i-1];
		a[i] = a[i]-c[i-1]*temp;
		Y[i] = Y[i]-Y[i-1]*temp;
	}

	/***** 上三角行列を解いてC[2][](=Xi)を求める *****/
//	C[2][N+1] = Y[N+1]/a[N+1];
	C[2*(N+2)+N+1] = Y[N+1]/a[N+1];
	for(i=N; 0<i; i--){
//		C[2][i] = (Y[i]-c[i]*C[2][i+1])/a[i];
		C[2*(N+2)+i] = (Y[i]-c[i]*C[2*(N+2)+i+1])/a[i];
	}

	/***** C[3][],C[4][]を求める *****/
	for(i=1; i<N+1; i++){
//		C[3][i] = (-2.*C[2][i]-C[2][i+1]+3.*(fi[i]-fi[i-1])/dx)/dx;
		C[3*(N+2)+i] = (-2.*C[2*(N+2)+i]-C[2*(N+2)+i+1]								+3.*(fi[i]-fi[i-1])/dx)/dx;
//		C[4][i] = (C[2][i]+C[2][i+1]-2.*(fi[i]-fi[i-1])/dx)/pow(dx,2);
		C[4*(N+2)+i] = (C[2*(N+2)+i]+C[2*(N+2)+i+1]								-2.*(fi[i]-fi[i-1])/dx)/pow(dx,2);
	}

	/***** 区間ごとの自然スプライン関数を求める *****/
	L2n = 0.;
	for(i=1; i<N+1; i++){
		for(j=0, x=xi[i-1]; j<nn+1; j++, x+=h){
//			p = C[1][i]+C[2][i]*(x-xi[i-1])+C[3][i]			//				*pow((x-xi[i-1]),2)+C[4][i]*pow((x-xi[i-1]),3);
			p = C[1*(N+2)+i]+C[2*(N+2)+i]*(x-xi[i-1])						+C[3*(N+2)+i]*pow((x-xi[i-1]),2)							+C[4*(N+2)+i]*pow((x-xi[i-1]),3);

		if(i==0 || i==nn*N)
			L2n += .5*h*pow((p-f(x)),2);
		else
			L2n += h*pow((p-f(x)),2);

		/**********************************************************
		下のコメント"//"を消去すると，P_i(x)を出力します．
		Excelなどでグラフ化するなら，>spline 12 > data.csv
		などとコマンドライン上でファイルにリダイレクトしてください．
		***********************************************************/
//			printf("%lf, %11.3e\n", x, p);
		}
	}

	printf("L2n=%11.3e\n",L2n);

	/***** メモリ領域の開放 *****/
	free(xi);free(fi);free(a);free(b);free(c);free(Y);free(C);
	return 0;
}