/* This is a C program to find the root of the 
non-linear equation x*sin(PI*x)-exp(-x)=0 using the Bisection 
algorithm, given an initial interval x1<= x <= x2. */    

#include<stdio.h>
#include<math.h>
#define PI (acos(-1.))

int main (void)
{
	double x1,x2,x3,f1,f2,f3;   
        int ib,n;
	double epsilon,error;
        printf("Enter the lower and upper bounds of the interval:\n");
        scanf("%lf%lf",&x1,&x2);
        printf("Specify convergence criterion as a fraction of initial "
               "interval\n");
        scanf("%lf",&epsilon);
        f1 = x1*sin(PI*x1)-exp(-x1);
        f2 = x2*sin(PI*x2)-exp(-x2); 

/* Check starting interval */ 
/* note: if reading input from file, could use getchar()!=EOF to test for
         end of user input.  See pg 51 of ABC */
        while(f1*f2>0) {
        	printf("f1*f2>0, try another interval\n");
                printf("Enter the lower and upper bounds of the interval:\n");
        	scanf("%lf%lf",&x1,&x2); 
                f1 = x1*sin(PI*x1)-exp(-x1); 
       	 	f2 = x2*sin(PI*x2)-exp(-x2);
	}   

	x3 = (x2+x1)/2.0;
        f3 = x3*sin(PI*x3)-exp(-x3);
	error = fabs(x2-x1);
	n = log(error/x3/epsilon)/log(2.0);
	if(fabs(x3)<epsilon)     /* divide by zero situation */
	        n = log(error/epsilon)/log(2.0);  
 /* Start bisections */ 
        for(ib=0;ib<=n;++ib) {
                if (f3*f1>0.0) {
		       x1=x3;
		       f1=f3;
		}
                else {
		       x2=x3;
		       f2=f3;
		} 
	 	x3 = (x1+x2)/2.0;
                f3 = x3*sin(PI*x3)-exp(-x3);
		error = fabs(x2-x1);
	}  
        printf("Iteration converged. Root is %e\n",x3);
        printf("The function value at root is %e\n",f3);
	printf("The convergence error is %e\n",error);

	return(0);
}