{******************************************************
*  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].  *
*                                                     *
*               Pascal 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           *
*                                                     *
******************************************************}
PROGRAM DEMO_ZNEWTON;
Uses WinCrt;

Var
        k,n : INTEGER;
        a,e,x0,y0,u,v,u1,v1,u2,v2,x,y : DOUBLE;


{*****************************************************
  Functions subroutine                               }
PROCEDURE Eval(x,y:double;VAR u,v,u1,v1,u2,v2:double);
Begin
  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
End;
{****************************************************}

{************************************************
*  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.            *
************************************************}
PROCEDURE ZNewton;
Label 100,fin;
Const TINY = 1e-12;
Begin
  k:=0;
100: k:=k+1;
  {Get u,v and the derivatives u1,v1,u2,v2}
  Eval(x0,y0,u,v,u1,v1,u2,v2);
  a:=u1*u1+u2*u2;
  {Guard against a=0}
  if a < TINY then
  begin
    writeln(' ZERO DIVIDE ERROR - u1*u1+u2*u2 must be <> 0.');
    goto fin
  end;
  x:=x0+(v*u2-u*u1)/a;
  y:=y0-(v*u1+u*u2)/a;
  {Check for convergence in euclidean space}
  if (x0-x)*(x0-x)+(y0-y)*(y0-y)<=e*e then goto fin;
  if k>=n then goto fin;
  x0:=x ; y0:=y;
  goto 100;
fin: End;


{main program}
BEGIN
  clrscr;
  writeln;
  writeln(' What is the initial guess:');
  writeln;
  write('    X0 = '); read(x0);
  write('    Y0 = '); read(y0);
  writeln;
  write(' Convergence criterion: '); read(e);
  write(' Maximum number of iterations: '); read(n);
  writeln;

  ZNewton;     {Call ZNewton subroutine}

  if a<>0 then
  begin
    writeln;
    writeln(' The root estimate is:');
    writeln;
    writeln('    X0 = ',x:9:6);
    writeln('    Y0 = ',y:9:6);
    writeln;
    writeln(' The number of iterations performed was: ',k);
    writeln
  end;
  ReadKey; DoneWinCrt
END.

{End of file znewton.pas}