//総和の計算
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

int main(void)
{
	int i;
	double a;
	a=0.0;

	for(i=1;i<=10;i++){
		a+=2.21;
		printf("i=\% d a=\% f\yen n",i,a);
	}
	printf("結果 a=\% f\yen n",a);

	system("PAUSE");
	return 0;
}


//偶数奇数の判定
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

int main(void)
{
	int i,a;

	a=0;

	for(i=0;i<=10;i++){
		a+=1.0;
		if(a%2==0){
			printf("偶数	i=%d	a=%d\n",i,a);
		}
		if(a%2==1){
			printf("奇数	i=%d	a=%d\n",i,a);
		}
	}
	
	system("PAUSE");
	return 0;
}


//1次元配列
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

int main(void)
{
	int i;
	double a[5],c;

	for(i=1;i<=4;i++){
		a[i]=2.0*i;
	}

	for(i=1;i<=4;i++){
		printf("i=%d	%f\n",i,a[i]);
	}

	c=a[2]+a[3];	
	printf("c=%f\n",c);

	system("PAUSE");
	return 0;
}


//2次元配列
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

int main(void)
{
	int i,j;
	double a[3][4],b[3][4],c[3][4];

	for(i=0;i<=2;i++){
		for(j=0;j<=3;j++){
			a[i][j]=2.0*i+5.0*j;
			printf("i=%d	j=%d	%f\n",i,j,a[i][j]);
		}
	}

	printf("a=\n");
	for(i=0;i<=2;i++){
		for(j=0;j<=3;j++){
			printf("%f	",a[i][j]);
		}
		printf("\n");
	}

	printf("\n");
	printf("b=\n");
	for(i=0;i<=2;i++){
		for(j=0;j<=3;j++){
			b[i][j]=2.5*i+4.0*j+1.2;
			printf("%f	",b[i][j]);
		}
		printf("\n");
	}

	printf("\n");
	printf("c=\n");
	for(i=0;i<=2;i++){
		for(j=0;j<=3;j++){
			c[i][j]=a[i][j]+b[i][j];
			printf("%f	",c[i][j]);
		}
		printf("\n");
	}

	system("PAUSE");
	return 0;
}


//シュウィンガー・ダイソン方程式
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define xc 50	//xの取る数
#define yc 150	//積分の分割数
#define alpha 0.3	//結合定数α
#define IR 0.0001	//積分の下限
#define UV 1.0	//積分の上限
#define cmax 100	//繰り替えす回数の上限

double funcB(double y,double x,double B,int cond);	//被積分関数
double sim(double a,double b,int n,double *sum);	//シンプソン積分公式のプログラム
double least(double *x,double *BB,int f);	//最小2乗法のプログラム


int main(void)
{
	int i,count,j;
	double x[xc+1],y1[2*yc+1],y2[2*yc+1],B1[2*yc+1],B2[2*yc+1],BB[xc+1],X1[2*yc+1],X2[2*yc+1];
	double a,b,B0;
	FILE *out;

	remove("./B.txt");remove("./ab.txt");

	printf("計算開始\n");

	for(i=0;i<=xc;i++){
		x[i]=IR+(UV-IR)*i/xc;	//変数xの値を入れている
	}

	for(count=1;count<=cmax;count++){

		for(i=0;i<=xc;i++){
			BB[i]=0.0;	//値を0に初期化するために最初に0の値を入れている
			for(j=0;j<=2*yc;j++){
				B1[j]=0.0;B2[j]=0.0;	//値を0に初期化するために最初に0の値を入れている
				y1[j]=IR+(x[i]-IR)*j/(2.0*yc);	//第一項の積分のyの値		
				y2[j]=x[i]+(UV-x[i])*j/(2.0*yc);	//第二項の積分のyの値	
				if(count==1){
					B1[j]=B2[j]=1.0;	//最初のBの定数
				}
				else{
					B1[j]=a*y1[j]+b;	//最小二乗法による第一項の積分範囲でのB
					B2[j]=a*y2[j]+b;	//最小二乗法による第二項の積分範囲でのB
				}
				X1[j]=funcB(y1[j],x[i],B1[j],1);	//第一項の各yでの被積分関数の値
				X2[j]=funcB(y2[j],x[i],B2[j],2);	//第二項の各yでの被積分関数の値
			}
			BB[i]=alpha*(sim(IR,x[i],yc,X1)+sim(x[i],UV,yc,X2))/4.0;	//積分実行
		}		

		printf("%d	%9f\n",count,BB[0]);

		a=least(x,BB,1);	//最小二乗法	B=ax+b
		b=least(x,BB,2);	//最小二乗法	B=ax+b
	}

	out=fopen("B.txt","a");
	for(i=0;i<=xc;i++){
		fprintf(out,"%f	%f\n",x[i],BB[i]);
	}
	fclose(out);

	out=fopen("ab.txt","a");
		fprintf(out,"a=%f	b=%f\n",a,b);
	for(i=0;i<=xc;i++){
		fprintf(out,"%f	%f\n",x[i],a*x[i]+b);
	}
	fclose(out);

	return 0;
}

double funcB(double y,double x,double B,int cond)
{
	double X;

	if(cond==1){
		X=3.0*y*B/((y+B*B)*x);	//第一項の被積分関数
	}
	if(cond==2){
		X=3.0*B/(y+B*B);	//第二項の被積分関数
	}

	return(X);
}

double sim(double a,double b,int n,double *sum)
{
// aは積分の下限、bは上限、nは分割数、sumは被積分関数の各値
	double S, h;
	int i;

	h=(b-a)/(2.0*n);

	S=(sum[0]+sum[2*n]);

	for(i=1;i<n;i++){
		S+=4.0*sum[2*i-1]+2.0*sum[2*i];
	}
		
	S+=4.0*sum[2*n-1];
	S*=h/3.0;
	return S;
}

double least(double *x,double *BB,int f)
{

	int n,i;
	double X1,X2,Y,XY,c;

	n=xc+1;
	X1=X2=Y=XY=0.0;
	for(i=0;i<=xc;i++){
		X1+=x[i];
		X2+=x[i]*x[i];
		Y+=BB[i];
		XY+=x[i]*BB[i];
	}

//a,bの場合わけ
	if(f==1){
		c=(n*XY-X1*Y)/(n*X2-X1*X1);	//a
	}
	if(f==2){
		c=(X2*Y-XY*X1)/(n*X2-X1*X1);	//b
	}

	return c;
}