#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>
typedef complex double Complex;

/* Implement Mueler's method from Burden, Faires
   and Burden, Numerical Analysis, 10th Edition,
   Chapter 2.6, Algorithm 2.8, page 97. */

Complex f(Complex z){
    return z*(z*(z*(z-3)+1)+1)+1;
}

Complex mueller(Complex p0,Complex p1,Complex p2,
    double TOL, int N0){
    Complex h1=p1-p0,h2=p2-p1;          // Step 1
    Complex delta1=(f(p1)-f(p0))/h1;
    Complex delta2=(f(p2)-f(p1))/h2;
    Complex d=(delta2-delta1)/(h2+h1);
    int i=3;
    while(i<=N0){                       // Step 2
        Complex b=delta2+h2*d;          // Step 3
        Complex D=csqrt(b*b-4*f(p2)*d);
        Complex E;                      // Step 4
        if(cabs(b-D)<cabs(b+D)) E=b+D;
        else E=b-D;
        Complex h=-2*f(p2)/E;           // Step 5
        Complex p=p2+h,z=f(p);
        printf("p=(%g,%g) f(p)=(%g,%g)\n",
        creal(p),cimag(p),creal(z),cimag(z));
        if(cabs(h)<TOL) return(p);      // Step 6
        p0=p1; p1=p2; p2=p;             // Step 7
        h1=p1-p0; h2=p2-p1;
        delta1=(f(p1)-f(p0))/h1;
        delta2=(f(p2)-f(p1))/h2;
        d=(delta2-delta1)/(h2+h1);
        i=i+1;
    }
    fprintf(stderr,
    "Method failed after %d iterations!\n",N0);
    exit(1);
}

int main(){
    Complex p0=0.5,p1=-0.5,p2=0;
    Complex r=mueller(p0,p1,p2,1e-5,15);
    printf("The root is (%g,%g)\n",creal(r),cimag(r));
    return 0;
}