/*T Math/CS 466/666 Midterm Solutions
\medskip\noindent
Problem 1(ii).  Write or modify a computer program to implement Mueller's
method and use it to approximate the solution to $1+z+z^2+z^3+z^4=0$
starting with an initial guess of $p_0=1$, $p_1=2$ and $p_2=3$.  Print the
first 6 iterations of the moethod, or in otherwords print $p_n$ for
$n=3,\ldots,8$. */

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

typedef complex double Complex;
Complex f(Complex z){
    //T Factor $1+z+z^2+z^3+z^4$ for efficiency and accuracy.
    return (((z+1)*z+1)*z+1)*z+1;
}
void printnp(int i,Complex p){
    printf("%3d %24.14e %24.14e\n",i,creal(p),cimag(p));
}

/*T The following code, based on Mueler's method given in Burden, Faires
and Burden, {\it Numerical Analysis}, 10th Edition, Chapter 2.6, Algorithm
2.8, page 97, has been modified to perform exactly 6 iterations. */
    
void mueller(Complex p0,Complex p1,Complex p2,
    double TOL, int N0){
    Complex h1=p1-p0,h2=p2-p1;                  //T Step 1
    Complex delta1=(f(p1)-f(p0))/h1;
    Complex delta2=(f(p2)-f(p1))/h2;
    Complex d=(delta2-delta1)/(h2+h1);
    for(int i=3;i<=8;i++){                      //T Step 2
        Complex b=delta2+h2*d;                  //T Step 3
        Complex D=csqrt(b*b-4*f(p2)*d);
        Complex E;                              //T Step 4
        if(cabs(b-D)<cabs(b+D)) E=b+D;
        else E=b-D;
        Complex h=-2*f(p2)/E;                   //T Step 5
        Complex p=p2+h,z=f(p);
        printnp(i,p);
        if(cabs(h)<TOL) return(p);              //T Step 6
        p0=p1; p1=p2; p2=p;                     //T 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);
    }
}

int main(){
    printf(
        "Math/CS 466/666 Midterm\nProblem 1(ii).\n\n"
        "%3s %24s %24s\n","n","real(pn)","imag(pn)");
    Complex p0=1.0,p1=2.0,p2=3.0;
    printnp(0,p0);
    printnp(1,p1);
    printnp(2,p2);
    mueller(p0,p1,p2,1e-5,15);
    return 0;
}