/****************************************************
*    Program to demonstrate the complex domain      *
*            Mueller's subroutine                   *
* ------------------------------------------------- *
* 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 complex root of zē + 1 = 0      *
*                                                   *
* SAMPLE RUN:                                       *
*                                                   *
* Write the initial guesses and their bounds:       *
*                                                   *
*    X0 = 0                                         *
*    Bond on X0 = 3                                 *
*    Y0 = 0                                         *
*    Bond on Y0 = 3                                 *
*                                                   *
* Convergence criterion: 0                          *
*                                                   *
* Maximum number of iterations: 100                 *
*                                                   *
* The calculated zero is (X,Y) =                    *
* -0.000014  -1.000034                              *
*                                                   *
* The associated Z value is (U,V) =                 *
* -0.000067   0.000028                              *
*                                                   *
* The number of steps was: 100                      *
****************************************************/
#include <stdio.h>
#include <math.h>

int    k,n;
double b1,b2,e,u,v,x,y,x0,yy0;


//**********************************************
//  Function subroutine                
void Z(double x,double y,double *u,double *v) {
  *u = x*x-y*y+1;
  *v = 2.0*x*y;
}
//**********************************************


/************************************************
* Mueller's method for complex roots subroutine *
* --------------------------------------------- *
* This routine uses the parabolic fitting tech- *
* nique associated with Mueller's method, but   *
* does it in the complex domain.                *
* --------------------------------------------- *
* INPUTS:                                       *
*   x0, y0 - The initial guess                  *
*   b1, b2 - A bound on the error in this guess *
*   e - The convergence criteria                *
*   n - The maximum number of iterations        *
* OUTPUTS:                                      *
*   x,y - The value of the complex root found   *
*   k - The number of iterations performed,     *
************************************************/
void ZMueller()  {
  //Labels: e100,e200,e300
  double d,d1,d2,x1,x2,x3,y1,y2,y3;
  double u1,u2,u3,v1,v2,v3,xl1,xl2;
  double a,a1,a2,b,c1,c2,e1,e2;
  // Start calculations
  k=1;
  x3=x0; y3=yy0; x1=x3-b1; y1=y3-b2; x2=x3+b1; y2=y3+b2;
e100: d=(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);
  // Avoid divide by zero
  if (d==0) d=1e-7;
  xl1=(x3-x2)*(x2-x1)+(y3-y2)*(y2-y1);
  xl1=xl1/d;
  xl2=(x2-x1)*(y3-y2)+(x3-x2)*(y2-y1);
  xl2=xl2/d;
  d1=(x3-x1)*(x2-x1)+(y3-y1)*(y2-y1);
  d1=d1/d;
  d2=(x2-x1)*(y3-y1)+(x3-x1)*(y2-y1);
  d2=d2/d;
  // Get function values
  Z(x1,y1,&u1,&v1);
  Z(x2,y2,&u2,&v2);
  Z(x3,y3,&u3,&v3);
  // Calculate Mueller parameters
  e1=u1*(xl1*xl1-xl2*xl2)-2.0*v1*xl1*xl2-u2*(d1*d1-d2*d2);
  e1=e1+2.0*v2*d1*d2+u3*(xl1+d1)-v3*(xl2+d2);
  e2=2.0*xl1*xl2*u1+v1*(xl1*xl1-xl2*xl2)-2.0*d1*d2*u2-v2*(d1*d1-d2*d2);
  e2=e2+u3*(xl2+d2)+v3*(xl1+d1);
  c1=xl1*xl1*u1-xl1*xl2*v1-d1*xl1*u2+xl1*d2*v2+u3*xl1;
  c1=c1-u1*xl2*xl2-v1*xl1*xl2+u2*xl2*d2+v2*d1*xl2-v3*xl2;
  c2=xl1*xl2*u1+xl1*xl1*v1-d2*xl1*u2-xl1*d1*v2+v3*xl1;
  c2=c2+u1*xl1*xl2-v1*xl2*xl2-u2*xl2*d1+v2*d2*xl2+u3*xl2;
  b1=e1*e1-e2*e2-4.0*(u3*d1*c1-u3*d2*c2-v3*d2*c1-v3*d1*c2);
  b2=2.0*e1*e2-4.0*(u3*d2*c1+u3*d1*c2+v3*d1*c1-v3*d2*c2);
  // Guard against divide by zero
  if (b1==0) b1=1e-7;
  a=atan(b2/b1); a=a/2.0;
  b=sqrt(sqrt(b1*b1+b2*b2));
  b1=b*cos(a); b2=b*sin(a);
  a1=(e1+b1)*(e1+b1)+(e2+b2)*(e2+b2);
  a2=(e1-b1)*(e1-b1)+(e2-b2)*(e2-b2);
  if (a1>a2) goto e200;
  a1=e1-b1; a2=e2-b2; goto e300;
e200: a1=e1+b1; a2=e2+b2;
e300: a=a1*a1+a2*a2;
  xl1=a1*d1*u3-a1*d2*v3+a2*u3*d2+a2*v3*d1;
  // Guard against divide by zero
  if (a==0) a=1e-7;
  xl1=-2.0*xl1/a;
  xl2=-d1*u3*a2+d2*v3*a2+a1*u3*d2+a1*v3*d1;
  xl2=-2.0*xl2/a;
  // Calculate new estimate
  x=x3+xl1*(x3-x2)-xl2*(y3-y2);
  y=y3+xl2*(x3-x2)+xl1*(y3-y2);
  // Test for convergence
  if (fabs(x-x0)+fabs(y-yy0)<e) return;
  // Test for number of iterations
  if (k>=n) return;
  // Continue
  k=k+1; x0=x; yy0=y; x1=x2; y1=y2; x2=x3; y2=y3; x3=x; y3=y;
  goto e100;
} // ZMueller()


void main()  {
  printf("\n Write the initial guesses and their bounds:\n\n");
  printf("    X0 = "); scanf("%lf",&x0);
  printf("    Bond on X0 = "); scanf("%lf",&b1);
  printf("\n");
  printf("    yy0 = "); scanf("%lf",&yy0);
  printf("    Bond on yy0 = "); scanf("%lf",&b2);
  printf("\n Convergence criterion: "); scanf("%lf",&e);
  printf("\n Maximum number of iterations: "); scanf("%d",&n);

  ZMueller();   // Call complex domain Mueller subroutine

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

// End of file zmueller.cpp