{ -------------------------- UNIT fbauhube.pas ----------------------------
* Ref.: "Numerical algorithms with C, By Gisela Engeln-Muellges and       *
*        Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].                  *
*                                                                         *
*                                       TPW Release By J-P Moreau, Paris. *
*                                               (www.jpmoreau.fr)         *
--------------------------------------------------------------------------}
UNIT FBauhube;

INTERFACE

Uses WinCrt;

Const
  ITERMAX = 1000;                    { Maximal number of function
                                       evaluations per root           }
  MACH_EPS = 1E-15;                  { Small real number              }
  EPSROOT = 31.622E6;                { SQRT(MACH_EPS)                 } 
  EPS = 64.0 * MACH_EPS;             { Accuracy for functional value  }
  BETA = 8.0 * EPS;
  QR = 0.1;                          { Real and imaginary parts of the}
  QI = 0.9;                          { spiralization constant         }

  NMAX = 100;

Type VEC = Array[0..NMAX] of Double;


  Function bauhub (         {Bauhuber's method for complex polynomials}
            real0,                 { Are the coefficients real ? .....}
            scale,                 { Scaling ? .......................}
            n: integer;            { degree of polynomial ............}
            ar,                    { Real parts of coefficients ......}
            ai: VEC;               { Imaginary parts, coefficients ...}
        Var rootr,                 { Real parts of roots .............}
            rooti,                 { Imaginary parts of roots ........}
            absf: VEC              { Absolute value of function values}
           ):Integer;

IMPLEMENTATION

  Procedure scpoly (n: Integer;      { length of vector .............}
                Var ar,              { Real part of the vector ......}
                    ai: VEC;         { Imaginary part of the vector .}
                Var scal: Double);   { Scaling factor ...............}
                    Forward;

  Function bauroot (n,                   { largest degree ...............}
                    iu: Integer;         { lowest degree ................}
                    ar,                  { Real parts of the coefficients}
                    ai: VEC;             { Imaginary parts, coefficients }
                Var x0r,                 { Real part of the root ........}
                    x0i:Double):Integer; { Imaginary part of the root ...}
                    Forward;

  Procedure chorner (n,              { highest degree in polynomial .}
                     iu: Integer;    { lowest degree in polynomial ..}
                     ar,             { Real parts of coefficients ...}
                     ai: VEC;        { Imaginary parts, coefficients }
                     xr,             { Real part of x ...............}
                     xi: Double;     { Imaginary part of x ..........}
                 Var pr,             { Real part of function value ..}
                     pi,             { Imaginary part of function v. }
                     p1r,            { Real part of first derivative }
                     p1i,            { Imaginary part, first deriv. .}
                     p2r,            { Real part, second derivative .}
                     p2i,            { Imaginary part, second deriv. }
                     rf1: Double);   { Error estimate for 1st deriv. }
                     Forward;

  Procedure polydiv (n,              { maximal degree ...............}
                     iu: Integer;    { minimal degree ...............}
                 Var ar,             { Real parts of coefficients ...}
                     ai: VEC;        { Imaginary parts, coefficients }
                     x0r,            { Real part of x ...............}
                     x0i: Double);   { Imaginary part of x ..........}
                     Forward;

 Function comabs        { Complex absolute value ....................}
              (
               ar,                { Real part .......................}
               ai: Double         { Imaginary part ..................}
              ): Double;
 {====================================================================*
 *                                                                    *
 *  Complex absolute value of   a                                     *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      ar,ai    REAL   ar, ai;                                       *
 *               Real, imaginary parts of  a                          *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      Absolute value of a                                           *
 *                                                                    *
 *   Macros used :    SQRT, ABS, SWAP                                 *
 *   =============                                                    *
 *                                                                    *
 *====================================================================}
 Var t: Double;
 Begin
 
  if (ar = 0.0) and (ai = 0.0) then
  begin
    comabs:=0.0;
    exit
  end;

  ar := ABS (ar);
  ai := ABS (ai);

  if ai > ar then                                  {Switch  ai and ar}
  begin
    t:=ai; ar:=ai; ai:=t
  end;

  if ai=0.0 then
    comabs:=ar
  else
    comabs:=ar * SQRT(1.0 + Sqr(ai/ar))
 End;

 Function comdiv         { Complex division ..........................}
           (
            ar,                    { Real part of numerator ..........}
            ai,                    { Imaginary part of numerator .....}
            br,                    { Real part of denominator ........}
            bi: Double;            { Imaginary part of denominator ...}
        Var cr,                    { Real part of quotient ...........}
            ci: Double             { Imaginary part of quotient ......}
           ): Integer;
 {====================================================================*
 *                                                                    *
 *  Complex division  c = a / b                                       *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      ar,ai    REAL   ar, ai;                                       *
 *               Real, imaginary parts of numerator                   *
 *      br,bi    REAL   br, bi;                                       *
 *               Real, imaginary parts of denominator                 *
 *                                                                    *
 *   Output parameters:                                               *
 *   ==================                                               *
 *      cr,ci    REAL   *cr, *ci;                                     *
 *               Real , imaginary parts of the quotient               *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      = 0      ok                                                   *
 *      = 1      division by 0                                        *
 *                                                                    *
 *   Macro used :     ABS                                             *
 *   ============                                                     *
 *                                                                    *
 *====================================================================}
Label fin;
Var
  tmp: Double;
Begin
  if (br = 0.0) and (bi = 0.0) then
  begin
    comdiv:=1;
    goto fin
  end;

  if ABS (br) > ABS (bi) then
  begin
    tmp  := bi / br;
    br   := tmp * bi + br;
    cr   := (ar + tmp * ai) / br;
    ci   := (ai - tmp * ar) / br
  end
  else
  begin
    tmp  := br / bi;
    bi   := tmp * br + bi;
    cr   := (tmp * ar + ai) / bi;
    ci   := (tmp * ai - ar) / bi
  end;

  comdiv:=0;
fin: End;

 Function Max(x,y: Double): Double;
 Begin
   If x>=y then
     Max:=x
   Else
     Max:=y
 End;

 {calculate y^x}
 FUNCTION Power(y,x: Double): Double;
 BEGIN
   IF x<0 THEN EXIT;
   Power:=Exp(x*Ln(y))
 END;


  Function bauhub (         {Bauhuber's method for complex polynomials}
            real0,                 { Are the coefficients real ? .....}
            scale,                 { Scaling ? .......................}
            n: integer;            { degree of polynomial ............}
            ar,                    { Real parts of coefficients ......}
            ai: VEC;               { Imaginary parts, coefficients ...}
        Var rootr,                 { Real parts of roots .............}
            rooti,                 { Imaginary parts of roots ........}
            absf: VEC              { Absolute value of function values}
           ):Integer;
 {====================================================================*
 *                                                                    *
 *  bauhub uses Bauhuber's Method to find all real or complex roots   *
 *  of a polynomial of degree n :                                     *
 *                                               n-1           n      *
 *      P(x) = a[0] + a[1] * x + ... + a[n-1] * x    + a[n] * x ,     *
 *                                                                    *
 *  where a[i], i=0, ..., n, can be complex.                          *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Application:                                                     *
 *   ===========                                                      *
 *      Find roots of arbitrary polynomials with complex coefficients.*
 *      If the polynomial roots are ill-condi=ditioned, i.e., if small*
 *      changes in the coefficients lead to large changes in the roots*
 *      the polynomial should not be scaled. Otherwise scaling helps  *
 *      with stability and performance.                               *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      real0    integer;                                             *
 *        = 0    Polynomial coefficients are complex                  *
 *       <> 0    Polynomial coefficients are real                     *
 *      scale    integer;                                             *
 *        = 0    no scaling                                           *
 *       <> 0    Scaling, see procedure scpoly()                      *
 *      n        integer;                                             *
 *               degree of the polynomial (must be >= 1)              *
 *      ar, ai   type VEC;                                            *
 *               Real and imaginary parts of the polynomial           *
 *               coefficients (ar[0], ..., ar[n])                     *
 *                                                                    *
 *   Output parameters:                                               *
 *   =================                                                *
 *      rootr    type VEC;   (Vector of length  n+1 )                 *
 *               rootr[0],..,rootr[n-1] are the real parts of the n   *
 *               roots                                                *
 *      rooti    type VEC;   (Vector of length n+1 )                  *
 *               rooti[0],..,rooti[n-1] are the imaginary parts       *
 *               of the roots                                         *
 *      absf     type VEC;                                            *
 *               absf[0],..,absf[n-1] are the magnitudes of the       *
 *               polynomial values at the computed roots              *
 *                                                                    *
 *   Return value:                                                    *
 *   ============                                                     *
 *      = 0      all is ok                                            *
 *      = 1      n < 1 or invalid input parameter                     *
 *      = 2      ar[n] = 0.0 and ai[n] = 0.0                          *
 *      = 3      Iteration maximum ITERMAX exceeded                   *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Functions or Procedures used:                                    *
 *   ============================                                     *
 *     bauroot():   determines one root of the polynomial             *
 *     scpoly():    Scales the polynomial                             *
 *     chorner():   Evaluates the polynomial                          *
 *     polydiv():   factors off a root                                *
 *     comabs():    magnitude of a complex number                     *
 *                                                                    *
 *====================================================================}
Label fin; {return label}
Var
  i, res: Integer;
  x0r, x0i, tempr, tempi, t1, t2, t3, t4, t5: Double;
  scalefak: Double;
Begin

  scalefak := 1.0; 
  if n < 1 then     {n is to small!}
  begin
    bauhub:=1;
    goto fin
  end;

  if (ar[n] = 0.0) and (ai[n] = 0.0) then  { Leading coefficient must}
  begin                                    { differ from zero        }
    bauhub:=2;
    goto fin
  end;                              

  for i := 0 to n do              { store given coefficients in root }
  begin
    rootr[i] := ar[i];
    rooti[i] := ai[i];
    if i < n then absf[i] := 0.0
  end;

  scalefak := 1.0;
  if scale <> 0 then                  { Scale polynomial, if desired }
    scpoly (n, rootr, rooti, scalefak);

  x0r := 0.0;                                       { Starting value }
  x0i := 0.0;

  for i := 0 to n-1 do
  begin                                       { compute the ith root }
    res := bauroot (n, i, rootr, rooti, x0r, x0i);

    rootr[i] := scalefak * x0r;                         { store root }
    rooti[i] := scalefak * x0i;

    if res<>0 then
    begin
      bauhub := 3;                      { Iteration maximum reached? }
      goto fin
    end;

    { Polynomial value of input polynomial at (rootr[i], rooti[i])   }

    chorner (n, 0, ar, ai, rootr[i], rooti[i],tempr, tempi, t1, t2, t3, t4, t5);

    absf[i] := comabs(tempr, tempi);            { store error        }

    polydiv (n, i, rootr, rooti, x0r, x0i);     { reduce degree      }

    if real0<>0 then                    { New starting value         }
      x0i := -x0i                       { depending on real x..      }
    else
    begin
      x0r := 0.0;
      x0i := 0.0
    end
  end;

  bauhub := 0;                                         { normal exit }

fin: End; {bauhub}


  Procedure scpoly (n: Integer;      { length of vector .............}
                Var ar,              { Real part of the vector ......}
                    ai: VEC;         { Imaginary part of the vector .}
                Var scal: Double);   { Scaling factor ...............}
 {====================================================================*
 *                                                                    *
 *  scalpoly scales the polynomial P :                                *
 *                                               n-1           n      *
 *      P(x) = a[0] + a[1] * x + ... + a[n-1] * x    + a[n] * x ,     *
 *                                                                    *
 *  where all a[i], i=0, ..., n, can be complex.                      *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Eingabeparameter:                                                *
 *   ================                                                 *
 *      n        integer;                                             *
 *               degree of the polynomila (>=1)                       *
 *      ar, ai   type VEC;                                            *
 *               Real and imaginary parts of the coefficients         *
 *               a[0],..,a[n]                                         *
 *                                                                    *
 *   Output parameters:                                               *
 *   =================                                                *
 *      ar, ai   type VEC;                                            *
 *               Real and imaginary parts of the coefficients         *
 *               a[0],..,a[n] of the scaled polynomial.               *
 *      scal     REAL                                                 *
 *               Scaling factor                                       *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Functions used:                                                  *
 *   ==============                                                   *
 *      comabs:  Modulus of a complex number                          *
 *      Power :  y^x                                                  *
 *      Max   :  Maximum of two real numbers                          *
 *                                                                    *
 *====================================================================}
Var
  p, pot: Double;
  i: Integer;
Begin

  scal := 0.0;
                              { scal =                               }
  p := comabs (ar[n], ai[n]); {              a[i]  1/(n-i)           }
  for i := 0 to n-1 do        {   max{ cabs( ---- )       , i=0..n-1 } 
                              {              a[n]                    }
    if (ar[i] <> 0.0) or (ai[i] <> 0.0) then
    begin
      ai[i] := ai[i] / p;
      ar[i] := ar[i] / p;

      pot := Power(comabs (ar[i],ai[i]), 1.0/(n-i));
      scal := Max(scal, pot)
    end;

  ar[n] := ar[n]/p;                    { Absolute value of a[n] = 1  }
  ai[n] := ai[n]/p;

  if scal = 0.0 then scal := 1.0;

  p:=1.0;
  for i := n-1 Downto 0 do
  begin
    p := p*scal;               {                    n-i              }
    ar[i]:=ar[i]/p;            { a[i] = a[i] / (scal    ), i=0..n-1  }
    ai[i]:=ai[i]/p             {                                     }
  end
End;


  Procedure chorner (n,              { highest degree in polynomial .}
                     iu: Integer;    { lowest degree in polynomial ..}
                     ar,             { Real parts of coefficients ...}
                     ai: VEC;        { Imaginary parts, coefficients }
                     xr,             { Real part of x ...............}
                     xi: Double;     { Imaginary part of x ..........}
                 Var pr,             { Real part of function value ..}
                     pi,             { Imaginary part of function v. }
                     p1r,            { Real part of first derivative }
                     p1i,            { Imaginary part, first deriv. .}
                     p2r,            { Real part, second derivative .}
                     p2i,            { Imaginary part, second deriv. }
                     rf1: Double);   { Error estimate for 1st deriv. }
 {====================================================================*
 *                                                                    *
 *  Horner scheme for polynomial with complex coefficients.           *
 *  We compute :                                                      *
 *    1. Polynomial value of the polynomial P (complex) of degree     *
 *       n - iu,                                                      *
 *    2. value of first derivative at x,                              *
 *    3. value of 2nd derivative at x,                                *
 *    4. an error estimate for the first derivative.                  *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;                                               *
 *               Maximal degree of the polynomial ( >= 1 )            *
 *      ar, ai   REAL   ar[], ai[];                                   *
 *               Real and imaginary parts of the coefficients of the  *
 *               polynomial with the coefficients a[iu], ..., a[n]    *
 *      x0r,x0i  REAL   x0r, x0i;                                     *
 *               Real and imaginary parts of the point of evaluation  *
 *                                                                    *
 *   Ausgabeparameter:                                                *
 *   ================                                                 *
 *      pr, pi   REAL   *pr, *pi;                                     *
 *               Real and imaginary part of the polynomial            *
 *      p1r, p1i REAL   *p1r, *p1i;                                   *
 *               Real and imaginary parts of the 1st derivative there *
 *      p2r, p2i REAL   *p2r, *p2i;                                   *
 *               Real and imaginary parts of the 2nd derivative       *
 *      rf1      REAL   *rf1;                                         *
 *               Error estimate for the first derivative              *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Functions used:                                                  *
 *   ==============                                                   *
 *     comabs():    modulus of a complex number                       *
 *     pow ():      Power function                                    *
 *====================================================================*
 *                                                                    *
 *   Constant used:  EPS                                              *
 *   =============                                                    *
 *                                                                    *
 *====================================================================}
Label continue, 10;
Var
    i,i1,j: Integer;
    temp, temp1: Double;
Begin

  p2r := ar[n];
  p2i := ai[n];

  pr := p2r;
  p1r := p2r;
  pi := p2i;
  p1i := p2i;

  rf1 := comabs(pr, pi);
  i1 := n - iu;

  j:=n-iu;
  for i := n - 1 Downto iu do
  begin
    temp := pr;                        { Polynomial value (pr,pi)    }
    pr := pr * xr - pi * xi + ar[i];
    pi := pi * xr + temp * xi + ai[i];
    if i = iu then goto continue;                      { exit i loop }

    temp := p1r;                       { 1st derivative (p1r,p1i)    }
    p1r := p1r * xr - p1i * xi;
    p1i := p1i * xr + temp * xi;

    temp := comabs(p1r, p1i);          { Error estimate for the 1st  }
    p1r := p1r + pr;                   { derivative of P             }
    p1i := p1i + pi;

    temp1 := comabs (pr, pi);
    temp := Max(temp, temp1);

    if temp > rf1 then
    begin
      rf1 := temp;
      i1 := j - 1
    end;

    if i - iu <= 1 then goto 10;                    {go to next i, j }

    temp := p2r;                         { 2nd derivative (p2r,p2i)  }
    p2r := p2r * xr - p2i * xi + p1r;
    p2i := p2i * xr + temp * xi + p1i;
10: Dec(j)
  end;

continue: temp := comabs(xr, xi);

  if temp <> 0.0 then
    rf1 := rf1 * Power(temp, 1.0*i1) * (i1 + 1)
  else
    rf1 := comabs(p1r, p1i);

  rf1 := rf1 * EPS;

  p2r := p2r + p2r;
  p2i := p2i + p2i

End; {chorner}

Procedure quadsolv       { Complex quadratic equation ................}
             (
               ar,                  {second degree coefficient .......}
               ai,
               br,                  {linear coefficient ..............}
               bi,   
               cr,                  {polynomial constant .............}
               ci: Double;
           Var tr,                  {solution ........................}
               ti: Double
             );
 {====================================================================*
 *                                                                    *
 *  Compute the least magnitude solution of the quadratic equation    *
 *  a * t**2 + b * t + c = 0. Here a, b, c and t are complex.         *
 *                                         2                          *
 *  Formeula used: t = 2c / (-b +/- sqrt (b  - 4ac)).                 *
 *  This formula is valid for a=0 .                                   *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *  Input parameters:                                                 *
 *  ================                                                  *
 *      ar, ai   coefficient of t**2             REAL   ar, ai;       *
 *      br, bi   coefficient of t                REAL   br, bi;       *
 *      cr, ci   constant term                   REAL   cr, ci;       *
 *                                                                    *
 *  Output parameter:                                                 *
 *  ================                                                  *
 *      tr, ti   complex solution of minimal magnitude                *
 *                                               REAL   *tr, *ti;     *
 *                                                                    *
 *  Macro used :     SQRT                                             *
 *  ============                                                      *
 *                                                                    *
 *====================================================================}
Var
   pr, pi, qr, qi, h: Double;
Begin
  pr := br * br - bi * bi;
  pi := 2.0 * br * bi;                       {  p = b * b             }

  qr := ar * cr - ai * ci;
  qi := ar * ci + ai * cr;                   {  q = a * c             }

  pr := pr - 4.0 * qr;       
  pi := pi - 4.0 * qi;                       { p = b * b - 4 * a * c  }

  h  := SQRT (pr * pr + pi * pi);            { q = sqrt (p)           }

  qr := h + pr;
  if qr > 0.0 then
    qr := SQRT (qr * 0.5)
  else
    qr := 0.0;

  qi := h - pr;
  if qi > 0.0 then
    qi := SQRT (qi * 0.5)
  else
    qi := 0.0;

  if pi < 0.0 then qi := -qi;

  h := qr * br + qi * bi;      { p = -b +/- q, choose sign for large }
                               { magnitude  p                        }
  if h > 0.0 then
  begin
    qr := -qr;
    qi := -qi
  end;

  pr := qr - br;
  pi := qi - bi;
  h := pr * pr + pi * pi;                      { t = (2 * c) / p     }

  if h = 0.0 then
  begin
    tr := 0.0;
    ti := 0.0
  end
  else
  begin
    tr := 2.0 * (cr * pr + ci * pi) / h;
    ti := 2.0 * (ci * pr - cr * pi) / h
  end
End; {quadsolv}

  Function bauroot (n,                   { largest degree ...............}
                    iu: Integer;         { lowest degree ................}
                    ar,                  { Real parts of the coefficients}
                    ai: VEC;             { Imaginary parts, coefficients }
                Var x0r,                 { Real part of the root ........}
                    x0i:Double):Integer; { Imaginary part of the root ...}
 {====================================================================*
 *                                                                    *
 *  bauroot computes one root of the polynomial P of degree n-iu:     *
 *                                                 n-iu               *
 *      P(x) = a[iu] + a[iu+1] * x + ... + a[n] * x                   *
 *                                                                    *
 *  with complex  coefficients a[i], i=iu, ..., n.                    *
 *  This program uses Newton's method on the function P(x) / P'(x).   *
 *  The iteration is stabilized using spiralization and extrapolation.*
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;                                               *
 *               Maximal degree of the polynomial ( >= 1 )            *
 *      iu       int iu;                                              *
 *               Index for the constant term of the polynomial,       *
 *               n-iu is the degree of the polynomial with            *
 *               coefficients a[iu], ..., a[n]                        *
 *      ar, ai   REAL   ar[], ai[];                                   *
 *               Real and imaginary parts of the coefficients         *
 *                                                                    *
 *   Output parameter:                                                *
 *   ================                                                 *
 *      x0r,x0i  REAL   *x0r, x0i;                                    *
 *               Real and imaginary part of the computed root         *
 *                                                                    *
 *   Return value :                                                   *
 *   =============                                                    *
 *      = 0      all is ok                                            *
 *      = 1      Division by zero                                     *
 *      = 2      Iteration number ITERMAX exceeeded                   *
 *      = 3      Improper input                                       *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Functions used:                                                  *
 *   ==============                                                   *
 *     chorner():   Complex Horner scheme                             *
 *     comabs():    magnitude of a complex number                     *
 *     comdiv():    Complex division                                  *
 *     quadsolv():  solves quadratic equations                        *
 *                                                                    *
 *   Constants used : TRUE, FALSE, ITERMAX,                           *
 *   ===============  QR, QI, MACH_EPS, EPS, EPSROOT, BETA            *
 *                                                                    *
 *====================================================================}
Label fin,10;
Var
    rc, result, endit, iter, i: Integer;

    xoldr, xoldi, xnewr, xnewi, h, h1, h2, h3, h4, dzmax, dzmin,
       dxr, dxi, tempr, tempi, abs_pold, abs_pnew, abs_p1new, temp,
       ss, u, v, bdze, pr, pi, p1r, p1i, p2r, p2i: Double;
Begin

  result:=2; endit:=0; iter:=0; i:=0;
  xoldr:=0.0; xoldi:=0.0; dxr:=0.0; dxi:=0.0; bdze:=0.0;

  if n < 1 then
  begin
    bauroot:=3;
    goto fin
  end;

  if n - iu = 1 then                        { Polynomial of degree 1 }
  begin
    quadsolv (0.0, 0.0, ar[n], ai[n], ar[n-1], ai[n-1], x0r, x0i);
    bauroot:=0;
    goto fin
  end;

  if n - iu = 2 then                        { Polynomial of degree 2 }
  begin
    quadsolv (ar[n],ai[n], ar[n-1],ai[n-1], ar[n-2],ai[n-2], x0r,x0i);
    bauroot:=0;
    goto fin
  end;

  xnewr := x0r;  xnewi := x0i;
  endit := 0;

  chorner (n, iu, ar, ai, xnewr, xnewi,     { Evaluate polynomial   }
           pr, pi, p1r, p1i, p2r, p2i, ss);

  Inc(iter);

  abs_pnew := comabs (pr, pi);
  if abs_pnew < EPS then
  begin
    bauroot:=0;                             { Starting value is a   }
    goto fin                                { good approximation    }
  end;

  abs_pold := abs_pnew;
  dzmin := BETA * (EPSROOT + comabs (xnewr, xnewi));

  while iter < ITERMAX do                       { Bauhuber-Iteration }
  begin
    abs_p1new := comabs (p1r, p1i);

    if abs_pnew > abs_pold then             { Spiralization step     }
    begin
      i := 0;                                { dx = dx * q            }
      Inc(iter);
      temp := dxr;

      dxr := QR * dxr - QI * dxi;
      dxi := QR * dxi + QI * temp
    end
    else
    begin
      dzmax := 1.0 + comabs(xnewr, xnewi);
      h1 := p1r * p1r - p1i * p1i - pr * p2r + pi * p2i;
      h2 := 2.0 * p1r * p1i - pr * p2i - pi * p2r;
      if (abs_p1new > 10.0 * ss) and (comabs(h1, h2) > 100.0 * ss * ss) then
      begin
        Inc(i);
        if i > 2 then i := 2;
        tempr := pr * p1r - pi * p1i;
        tempi := pr * p1i + pi * p1r;

        rc := comdiv(-tempr, -tempi, h1, h2, dxr, dxi);
        if rc <> 0 then
        begin
          bauroot := 1;
          goto fin
        end;

        if comabs (dxr, dxi) > dzmax then
        begin
          temp := dzmax / comabs (dxr, dxi);       { Newton step     }
          dxr := dxr * temp;
          dxi := dxi * temp;
          i := 0
        end;

        if (i = 2) and (comabs(dxr, dxi) < dzmin / EPSROOT) and (comabs(dxr, dxi) > 0.0) then
        begin
          i := 0;                            { Extrapolation step    }
          rc := comdiv(xnewr - xoldr, xnewi - xoldi, dxr, dxi, h3, h4);
          if rc <> 0 then
          begin
            bauroot := 1;
            goto fin
          end;

          h3 := h3 + 1.0;
          h1 := h3 * h3 - h4 * h4;
          h2 := 2.0 * h3 * h4;
          rc := comdiv (dxr, dxi, h1, h2, h3, h4);
          if rc <> 0 then
          begin
            bauroot := 1;
            goto fin
          end;

          if comabs(h3, h4) < 50.0 * dzmin then
          begin
            dxr := dxr + h3;
            dxi := dxi + h4
          end
        end;

        xoldr := xnewr;
        xoldi := xnewi;
        abs_pold := abs_pnew
      end
      else
      begin
        i := 0;                            { Close to a saddle point }
        h := dzmax / abs_pnew;
        dxr := h * pr;
        dxi := h * pi;

        xoldr := xnewr;
        xoldi := xnewi;
        abs_pold := abs_pnew;

        Repeat

          chorner (n, iu, ar, ai, xnewr+dxr, xnewi+dxi, u, v, h, h1, h2, h3, h4);
          Inc(iter);

          dxr := dxr + dxr;
          dxi := dxi + dxi;                   { dx = dx * 2.0         }

        Until ABS(comabs(u,v)/abs_pnew-1.0) >= EPSROOT
      end
    end;

    if endit<>0 then
    begin
      if comabs(dxr, dxi) < 0.1 * bdze then
      begin
        xnewr := xnewr + dxr;
        xnewi := xnewi + dxi
      end;

      result := 0;
      goto 10                                 { stop iteration   }
    end
    else
    begin
      xnewr := xoldr + dxr;
      xnewi := xoldi + dxi;
      dzmin := BETA * (EPSROOT + comabs (xnewr, xnewi));

      chorner (n, iu, ar, ai, xnewr, xnewi, pr, pi, p1r, p1i, p2r, p2i, ss);

      Inc(iter);
      abs_pnew := comabs( pr, pi);

      if abs_pnew = 0.0 then
      begin
        result := 0;
        goto 10
      end;

      if (comabs(dxr, dxi) < dzmin) or (abs_pnew < EPS) then
      begin
        endit := 1;
        bdze := comabs(dxr, dxi)
      end
    end
  end;  {End Bauhuber iteration}

10: x0r := xnewr;
    x0i := xnewi;

  bauroot:=result;
fin: End; {bauroot}


  Procedure polydiv (n,              { maximal degree ...............}
                     iu: Integer;    { minimal degree ...............}
                 Var ar,             { Real parts of coefficients ...}
                     ai: VEC;        { Imaginary parts, coefficients }
                     x0r,            { Real part of x ...............}
                     x0i: Double);   { Imaginary part of x ..........}
 {====================================================================*
 *                                                                    *
 *  polydiv computes the coefficients of the polynomial Q, with       *
 *  P(x) = Q(x) * (x - x0), where x0 is a computed root of P.         *
 *  Both P and Q and x0 may be complex.                               *
 *                                                                    *
 *====================================================================*
 *                                                                    *
 *   Input parameters:                                                *
 *   ================                                                 *
 *      n        int n;                                               *
 *               Highest degree of the polynomial ( >= 1 )            *
 *      ar, ai   REAL   ar[], ai[];                                   *
 *               Real and imaginary parts of the coefficienten of P   *
 *               of degree n-iu, with a[iu], ..., a[n]                *
 *      x0r,x0i  REAL   x0r, x0i;                                     *
 *               Real and imaginary parts of the root x0              *
 *                                                                    *
 *   Output parameter:                                                *
 *   ================                                                 *
 *      ar, ai   REAL   ar[], ai[];                                   *
 *               Real and imaginary parts of the coefficients         *
 *               ar[iu+1],..,ar[n] of the remainder polynomial Q      *
 *                                                                    *
 *====================================================================}
 Var
    i: Integer;
    temp: Double;
 Begin
   for i := n - 1 Downto iu+1 do
   begin
     temp := ar[i+1];
     ar[i] := ar[i] + temp * x0r - ai[i+1] * x0i;
     ai[i] := ai[i] + ai[i+1] * x0r + temp * x0i
   end
 End;

END.

{ ------------------------- END fbauhube.pas  ------------------------ }