/****************************************************
*     Program to demonstrate Mueller's method       *
* ------------------------------------------------- *
* Reference: BASIC Scientific Subroutines, Vol. II  *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. *
*                                                   *
*                C++ version by J-P Moreau, Paris.  *
*                       (www.jpmoreau.fr)           *
* ------------------------------------------------- *
* Example:  Find a real root of f(x)=(x+1)^5        *
*                                                   *
* Sample run:                                       *
*                                                   *
* Input the initial guess:                          *
*    X0 = 0                                         *
* What is the bond of this guess: 3                 *
* Error criterion: 1e-6                             *
* Maximum number of iterations: 100                 *
*                                                   *
* The calculated zero is X = -1.000476              *
* The associated value is Y =  0.000000             *
* The number of steps was: 47                       *
*                                                   *
****************************************************/
#include <stdio.h>
#include <math.h>

// global variables
int    k,n;
double d,e,x,x0;


//***************************************************
//  Function subroutine                          
double Y(double x) {
  return (1+5*x+10*x*x+10*x*x*x+5*x*x*x*x+x*x*x*x*x);
}
//***************************************************


/**********************************************
*         Mueller's method subroutine         *
* ------------------------------------------- *
* This routine iteratively seeks the root of  *
* a function Y(x) by fitting a parabola to    *
* three points and calculating the nearest    *
* root as described in Becket and Hurt, nume- *
* rical calculations and algorithms.          *
* writeS:                                     *
* The routine requires an initial guess, x0,  *
* a bound on the error on this guess, d and a *
* convergence criteria, e. Also required is a *
* limit on the number of iterations, n.       *
* OUTPUTS:                                    *
* The routine returns the value of the root   *
* found, x and the number of iterations per-  *
* formed, k.                                  *
**********************************************/
void Mueller() {
  //Labels: e100,e200
  double a1,b,c1,d1,e1,e2,e3,x1,x2,x3,xl,xl1;
  // Set up the three evaluation points
  k=1; x3=x0; x1=x3-d; x2=x3+d;
  // Calculate Mueller parameters and guard against divide by zero
  if (x2-x1==0) x2=x2*1.0000001;
e100: if (x2-x1==0) x2=x2+0.0000001;
  xl1=(x3-x2)/(x2-x1);
  d1=(x3-x1)/(x2-x1);
  if (k>1) goto e200;
  // Get values of function
  e1=Y(x1); e2=Y(x2);
e200: e3=Y(x3);
  a1=xl1*xl1*e1-d1*d1*e2+(xl1+d1)*e3;
  c1=xl1*(xl1*e1-d1*e2+e3);
  b=a1*a1-4.0*d1*c1*e3;
  // Test for complex root, meaning the parabola is inverted
  if (b<0) b=0.0;
  // Choose closest root
  if (a1<0) a1=a1-sqrt(b);
  if (a1>0) a1=a1+sqrt(b);
  // Guard against divide by zero
  if (a1==0) a1=1e-7;
  xl=-2.0*d1*e3/a1;
  // Calculate next estimate
  x=x3+xl*(x3-x2);
  // Test for convergence
  if (fabs(x-x3)<e) return;
  // Test for number of iterations
  if (k>=n) return;
  // Otherwise, make another pass
  k++;
  // Some saving
  x1=x2; x2=x3; x3=x; e1=e2; e2=e3;
  goto e100;
} // Mueller()


void main() {
  printf("\n Input the initial guess:\n");
  printf("\n    X0 = "); scanf("%lf",&x0);
  printf("\n What is the bond of this guess: "); scanf("%lf",&d);
  printf("\n Error criterion: "); scanf("%lf",&e);
  printf("\n Maximum number of iterations: "); scanf("%d",&n);

  Mueller();  // Call Mueller's routine

  printf("\n\n The calculated zero is X = %f\n", x);
  printf("\n The associated Y value is Y = %f\n", Y(x));
  printf("\n The number of steps was: %d\n\n", k);
}

// End of file Mueller.cpp