/******************************************************
*  Program to demonstrate the Newton root subroutine  *
*             in the complex domain                   *
* --------------------------------------------------- *
*  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)            *
* --------------------------------------------------- *
* SAMPLE RUN:                                         *
*                                                     *
* (Example: find a complex root of  z^2 + 1 = 0)      *
*                                                     *
* What is the initial guess:                          *
*                                                     *
*    X0 = .2                                          *
*    Y0 = 1.2                                         *
*                                                     *
* Convergence criterion: 1e-8                         *
* Maximum number of iterations: 10                    *
*                                                     *
* The root estimate is:                               *
*                                                     *
*    X0 =  -0.000000                                  *
*    Y0 =  1.000000                                   *
*                                                     *
* The number of iterations performed was: 5           *
*                                                     *
******************************************************/
#include <stdio.h>
#include <math.h>

#define TINY  1e-12  //small number

int    k,n;
double a,e,x0,yy0,u,v,u1,v1,u2,v2,x,y;


//*********************************************************
//  Functions subroutine
void Eval(double x, double y, double *u, double *v,
		  double *u1, double *v1, double *u2, double *v2) {
  *u=x*x-y*y+1.0; *v=2.0*x*y;
  *u1=2.0*x; *u2=-2.0*y;
  *v1=2.0*y; *v2=2.0*x;
}
//*********************************************************

/************************************************
*  Complex root seeking using Newton"s method   *
* --------------------------------------------- *
* This routine uses the complex domain form of  *
* Newton"s method for iteratively searching     *
* for roots. The complex function and its first *
* partial derivatives must be available in the  *
* form; F(Z) = U(X,Y) + I V(X,Y). The required  *
* derivatives are DU/DX and DU/DY.              *
* --------------------------------------------- *
* INPUTS; initial guess, x0 and y0, convergence *
* criteria, e, maximum number of iterations, n. *
* OUTPUTS; approximation to the root, x and y,  *
* number of performed iterations, k.            *
************************************************/
void ZNewton()  {
  // Label: e100
  k=0;
e100: k++;
  // Get u,v and the derivatives u1,v1,u2,v2
  Eval(x0,yy0,&u,&v,&u1,&v1,&u2,&v2);
  a=u1*u1+u2*u2;
  // Guard against a=0
  if (a < TINY) {
    printf("\n ZERO DIVIDE ERROR - u1*u1+u2*u2 must be <> 0.\n\n");
    return;
  }
  x=x0+(v*u2-u*u1)/a;
  y=yy0-(v*u1+u*u2)/a;
  // Check for convergence in euclidean space
  if ((x0-x)*(x0-x)+(yy0-y)*(yy0-y)<=e*e) return;
  if (k>=n) return;
  x0=x ; yy0=y;
  goto e100;
} // ZNewton()


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

  ZNewton();     // Call ZNewton subroutine

  if (a!=0) {
    printf("\n\n The root estimate is:\n\n");
    printf("    X0 = %f\n",x);
    printf("    Y0 = %f\n",y);
    printf("\n The number of iterations performed was: %d\n\n",k);
  }
}


// End of file Znewton.cpp