/*
 * Bisection.c: Bisection Method 二分法
 */
#include <stdio.h>
#include <math.h>

double f(double x)
/*  f(x) = x^2-2  */
{return x*x-2;}

int main()
{
    double x, x1, x2, err;
    double epsa = 1e-70, epsr = 1e-6;

    x1 = 1;
    x2 = 2;
    do {
        x = (x1+x2)/2;
        if (f(x)>0) x2 = x;
        else        x1 = x;
        err = fabs(x2-x1)/(epsa+epsr*fabs(x));
        printf(" x =  %e err =  %e\n",x,err);
    } while (err >= 1);
    putchar('\n');
    printf(" sqrt(2) =  %e\n", x);
    return 0;
}
/*   ---   実行例: Pascalで学ぶプログラミング p.57   ---   */
/*
Z:\src\C>Bisection
 x =  1.500000e+00 err =  3.333333e+05
 x =  1.250000e+00 err =  2.000000e+05
 x =  1.375000e+00 err =  9.090909e+04
 x =  1.437500e+00 err =  4.347826e+04
 x =  1.406250e+00 err =  2.222222e+04
 x =  1.421875e+00 err =  1.098901e+04
 x =  1.414062e+00 err =  5.524862e+03
 x =  1.417969e+00 err =  2.754821e+03
 x =  1.416016e+00 err =  1.379310e+03
 x =  1.415039e+00 err =  6.901311e+02
 x =  1.414551e+00 err =  3.451847e+02
 x =  1.414307e+00 err =  1.726221e+02
 x =  1.414185e+00 err =  8.631852e+01
 x =  1.414246e+00 err =  4.315740e+01
 x =  1.414215e+00 err =  2.157916e+01
 x =  1.414200e+00 err =  1.078970e+01
 x =  1.414207e+00 err =  5.394820e+00
 x =  1.414211e+00 err =  2.697403e+00
 x =  1.414213e+00 err =  1.348700e+00
 x =  1.414214e+00 err =  6.743493e-01

 sqrt(2) =  1.414214e+00
*/