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

#define tolerance 1.0e-8
#define n 100
#define PI 3.14159265358

void gauss_seidel(double x[],double r[]);
double norm(double x[]);

int main(void)
{
	int i,m;
	double residual;
	double x[n],r[n];
	
	//initial guess
	for (i=0;i<n;i++)
		x[i]=(double) rand()/RAND_MAX;

	//the first iteration to get the first residaul
	gauss_seidel(x,r);
        residual=norm(r);

	m=0;
	while(residual >tolerance)
	{
		gauss_seidel(x,r);
		residual=norm(r);
		m=m+1;
	}	
	for (i=0;i<n;i++)
		printf("x[%d]=%.6f\n",i,x[i]);	
	
	printf("number of iteration = %d \n", m);
	return 0;
}
    

void gauss_seidel(double x[],double r[])
{
	int i;
	double bi;
	double aii=2.01;
	
	//compute the next iterate and update
	for (i=0;i<n;i++)
	{
		bi=sin(2.*PI*(i+1)/50.)/100.;

		if (i==0)
			x[i]=(bi+x[i+1])/aii;
		else if (i==n-1)
			x[i]=(bi+x[i-1])/aii;
		else
			x[i]=(bi+x[i-1]+x[i+1])/aii;
	}
	
	//compute the residual
	for(i=0;i<n;i++)
	{
		bi=sin(2.*PI*(i+1)/50.)/100.;
		if (i==0)
			r[i]=bi+x[i+1]-aii*x[i];
		else if (i==n-1)
			r[i]=bi+x[i-1]-aii*x[i];
		else
			r[i]=bi+x[i-1]+x[i+1]-aii*x[i];
	}
}          
double norm(double x[])
{
	int i;
	double dis=0.0;
	for (i=0;i<n;i++)
		dis=dis+x[i]*x[i];
	return sqrt(dis);
}