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

/****************** 表2.2の誤差を出力するプログラム *************************

	・表2.2のN一行のみの出力を行います．

	・実行すると，1，2，3のいずれかのデータの入力を求めます．

	・補間点を0点とするなら1，(2.78)とするなら2，(2.80)とする
		なら3を入力してください．

	・コマンドラインの第一引数に1,2,3のいずれかを与えることでも
		即時に誤差を出力することができます．
		例）	>cheby 1

	・ただし，下記のNを変更する場合はコンパイルしなおしてください．

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

	・N元連立一次方程式の解法には，
		ガウスの単純消去法(前進消去，後退代入)と
			ピボット選択法を用います(5.2節)．

	・N元連立一次方程式の枢軸要素(ピボット)の絶対値を最大となるように
		選択して誤差が広がるのを防ぎます．この処理がなければ，
			Nを大きくするにつれて，誤差が目に見えて広がります．

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

#define N 49			//連立方程式の元数(textのNに対応)
#define Pi acos(-1.)		//πの値
#define n 1000
#define h 2./n			//台形公式(式(3.4))のh
#define eps pow(10, -15)	//ピボットの許容最小値，0除算の判定に使用

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

double fTN(double x, double ck[]){		//式(2.77)
	int k;
	double fN=0.;
	for(k=0; k<N; k++){
		fN += ck[k]*cos(k*acos(x));
	}
	return fN;
}

/******************************************************************
 **** ピボット交換処理と上三角行列への変換(前進消去)を行う関数 ****
 *****************************************************************/

int PIV(double a[N+1][N], double yi[]){
	int i, ii, k, ip[N+1];
	double pivot, w, temp;

	/***** 交換情報を保持する配列ip[]の初期化 *****/
	for(i=1; i<N+1; i++)
		ip[i] = 0;

	for(i=1; i<N; i++){

		/***** ピボットに最初の比較対象である対角要素を代入 *****/
		pivot = fabs(a[i][i-1]);
		for(ii=i; ii<N+1; ii++){
			if(pivot < fabs(a[ii][i-1])){
				pivot = fabs(a[ii][i-1]);

/***** i行目のピボットの交換位置(同じ列で絶対値が最大となる行)を記録 *****/
				ip[i] = ii;
			}
		}

		/***** ピボット交換処理 *****/
		if(ip[i]>0){		
			for(k=i-1; k<N; k++){
				w = a[i][k];
				a[i][k] = a[ip[i]][k];
				a[ip[i]][k] = w;
			}
			w = yi[i];
			yi[i] = yi[ip[i]];
			yi[ip[i]] = w;
		}			

		/***** 行列を上三角行列に直す処理 *****/
		for(ii=i+1; ii<N+1; ii++){			//行変数

			/***** ピボットが0でないことを確認 *****/
			if(fabs(a[i][i-1]) < eps){
				printf("divided by zero!\n");
				return 0;
			}
			temp = a[ii][i-1]/a[i][i-1];
			for(k=i-1; k<N; k++)			//列変数
				a[ii][k] += -a[i][k]*temp;
			yi[ii] += -yi[i]*temp;
		}
	}
	return 1;	//正常に処理を終了したときのみ，1(真)を返す
}

int main(int argc, char **argv){
	int i, k, swtch;
	double theta[N+1], xi[N+1], Tk[N+1], ck[N];
	double wi, lambdak, x, p, L2n;
	double a[N+1][N], yi[N+1];

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

	if(!(argc ==1 || argc ==2)){
		printf("引数のエラーです\n");
		return 1;
	}
	if(argc == 1){
		printf("補間点を0点とするなら1，(2.78)とするなら2，(2.80)とする				なら3を入力してください\n switch = ");
		scanf("%d",&swtch);
	}else{
		argv++;
		swtch = atoi(*argv);
	}
	if(!(swtch == 1 || swtch == 2 || swtch ==3)){
		printf("入力エラーです\n");
		return 2;
	}

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

	for(i=1; i<N+1; i++){
		switch (swtch){
			case 1 :	//補間点0点
				theta[i] = Pi*(2.*(double)i-1.)/(2.*N);
				break;
			case 2 :	//補間点(2.78)
				theta[i] = Pi*(double)i/(N+1.);
				break;
			case 3 :	//補間点(2.80)
				theta[i] = Pi*((double)i-1)/(N-1.);
				break;
			default :
				return 3;
		}
		xi[i] = cos(theta[i]);
		yi[i] = f(xi[i]);
	}

	if(swtch == 1){				//補間点を0点とするとき
		wi = Pi/N;
		for(k=0; k<N; k++){
			if(k==0)
				lambdak = Pi;
			else
				lambdak = Pi/2;
			ck[k] = 0.;
			for(i=1; i<N+1; i++){
				Tk[i] = cos(k*theta[i]);
				ck[k] += wi*Tk[i]*yi[i];	//式(2.74)部分
			}
			ck[k] /= lambdak;			//式(2.74)
		}
	}else{

	/***** N元一次連立方程式の係数行列a[i][k] *****/
	/***** i=1,2,3,...,N,	k=0,1,2,...,N-1   *****/
	for(i=1; i<N+1; i++)
		for(k=0; k<N; k++)
			a[i][k] = cos((double)k*theta[i]);

	/***** ピボット交換と前進消去 *****/
	if( PIV(a, yi) != 1 )	return 4;	
	/***** 関数の返り値が1(真)以外のとき，エラーコードを返して中断 *****/

	/***** 後退代入を繰り返し，ck[]を求める *****/
	if(fabs(a[N][N-1]) < eps){
		printf("divided by zero!\n");
		return 5;
	}
	ck[N-1] = yi[N]/a[N][N-1];
	for(i=N-1; 0<i; i--){
		p = 0.;
		for(k=i; k<N; k++)
			p += a[i][k]*ck[k];
		ck[i-1] = (yi[i]-p)/a[i][i-1];
	}
	}

/***** L_2ノルムの2乗(式(2.59))を積分の台形公式(式(3.4))を用いて求める *****/
	L2n = 0.;

	/************************************************************
	下のコメント"//"を消去すると，一定間隔でf^T_N(x)を出力します．
	Excelなどでグラフ化するなら，>cheby 2 > data.csv　などと
	コマンドライン上でファイルにリダイレクトしてください．
	*************************************************************/
	for(i=0; i<n+1; i++){
		x = -1.+h*(double)i;
//		if(fmod(i,(int)(n/10))==0)					//			printf("%6.3lf,%11.3e\n", x, fTN(x,ck));

		if(i == 0 || i == n)
			L2n += .5*h*pow((fTN(x,ck)-f(x)),2);
		else
			L2n += h*pow((fTN(x,ck)-f(x)),2);

	}

	printf("switch=%d, N=%4d, L2n=%11.3e\n", swtch, N, L2n);
	return 0;
}