/*--------------------------------------------------------*
* Main Program                                            *
*    This program evaluates root of a polynomial using    *
*    Muller's method                                      *
*---------------------------------------------------------*
* Functions invoked                                       *
*    NIL                                                  *
*---------------------------------------------------------*
* Subroutines used                                        *
*    F(x)                                                 *
*---------------------------------------------------------*
* Variables used                                          *
*    x1,x2,x3 - initial values                            *
*    f1,f2,f3 - function values at x1,x2,x3               *
*    a0,a1,a2 - coefficients of quadratic polynomial      *
*    hi - xi-x3                                           *
*    di - function difference fi-f3                       *
*---------------------------------------------------------*
* Constants used                                          *
*    EPD  -  Error bound                                  *
*---------------------------------------------------------*/

#include<math.h>
#define EPS 0.000001

main()
{
    float F(float x);
    float x1,x2,x3,x4,f1,f2,f3,f4,h1,h2,d1,d2,a0,a1,a2,h;
    printf("\nInput three initial points \n");
    scanf("%f %f %f",&x1,&x2,&x3);

    f1 = F(x1);
    f2 = F(x2);
    f3 = F(x3);

    begin:
    h1 = x1-x3;
    h2 = x2-x3;

    d1 = f1-f3;
    d2 = f2-f3;
/* Compute parameters a0,a1,a2  */

    a0 = f3;
    a1 = (d2*h1*h1-d1*h2*h2)/(h1*h2*(h1-h2));
    a2 = (d1*h2-d2*h1)/(h*h2*(h1-h2));

/* Compute h  */

    if(a1>0.0)
        h=(-2.0*a0)/(a1+sqrt(a1*a1-4*a2*a0));
    else
        h=(-2.0*a0)/(a1-sqrt(a1*a1-4*a2*a0));

/* Compute x4 and f4 */

    x4 = x3+h;
    f4 = F(x4);

/* Test for accuracy  */

    if(f4<=EPS)  /*  root obtained */
    {
        printf("\n\nROOT BY MULLER'S METHOD\n\n");
        printf("Root = %f\n",x4);
        printf("\n");
    }
    else
    {
        x1 = x2;
        x2 = x3;
        x3 = x4;
        f1 = f2;
        f2 = f3;
        f3 = f4;
        goto begin;
    }
}
/* End of main() program  */
/* -----------------------------------------------------------*/
/* Defining the subroutine F(x)   */
float F(float x)
{
    float f;
    f = x*x*x+2*x*x+10*x-20;
    return (f);
}
/* End of subroutine ()  */
/* -----------------------------------------------------------*/