DECLARE FUNCTION comabs# (ar#, ai#)
DECLARE FUNCTION icomdiv% (ar#, ai#, br#, bi#, cr#, ci#)
DECLARE FUNCTION ibauhub% (real0 AS INTEGER, scale AS INTEGER, n%, ar#(), ai#(), rootr#(), rooti#(), absf#())
DECLARE SUB scpoly (n%, ar#(), ai#(), scal#)
DECLARE SUB chorner (n%, iu%, ar#(), ai#(), xr#, xi#, pr#, pi#, p1r#, p1i#, p2r#, p2i#, rf1#)
DECLARE SUB polydiv (n%, iu%, ar#(), ai#(), x0r#, x0i#)
DECLARE SUB quadsolv (ar#, ai#, br#, bi#, cr#, ci#, tr#, ti#)
DECLARE FUNCTION ibauroot% (n%, iu%, ar#(), ai#(), x0r#, x0i#)
'***************************************************************************
'*               Test program for Bauhuber's method                        *
'* ----------------------------------------------------------------------- *
'*  This program 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.                               *
'* ----------------------------------------------------------------------- *
'* SAMPLE RUN:                                                             *
'* (Find all roots of polynomial of degree 4:                              *
'*  P(x) = 0.000089248 - 0.04368 X + 2.9948 X2 - 6.798 X3 + X4)            *
'*                                                                         *
'*  ----------------------------------------------------------             *
'*  Complex and Real Roots of a Polynomial (Bauhuber's method)             *
'*  ----------------------------------------------------------             *
'*  Polynomial coefficients:                                               *
'*    0     .000089248     0                                               *
'*    1    -.04368         0                                               *
'*    2     2.9948         0                                               *
'*    3    -6.798          0                                               *
'*    4     1              0                                               *
'*  Roots (without scaling)                                                *
'*    No   Real part    Imaginary part   Function value                    *
'*  0        0.002454       0.000000   2.055622946982948D-15               *
'*  1        0.012573       0.000000   2.05556701633865D-15                *
'*  2        0.457319       0.000000   2.07678617766666D-15                *
'*  3        6.325654       0.000000   5.019937198408532D-14               *
'*  Roots (with scaling)                                                   *
'*    No   Real part    Imaginary part   Function value                    *
'*  0        0.002454       0.000000   2.055638722971591D-15               *
'*  1        0.012573       0.000000   2.055715472097546D-15               *
'*  2        0.457319       0.000000   2.205981399958949D-15               *
'*  3        6.325654       0.000000   5.019937198408532D-14               *
'*  ----------------------------------------------------------             *
'*                                                                         *
'* ----------------------------------------------------------------------- *
'* Ref.: "Numerical algorithms with C, By Gisela Engeln-Muellges and       *
'*        Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].                  *
'*                                                                         *
'*                                     Basic Release By J-P Moreau, Paris. *
'*                                              (www.jpmoreau.fr)          *
'***************************************************************************
'Program TBauhube
DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0

  F$ = "######.######"

  CLS
  PRINT
  PRINT "  ----------------------------------------------------------"
  PRINT "  Complex and Real Roots of a Polynomial (Bauhuber's method)"
  PRINT "  ----------------------------------------------------------"

  n = 4      'order of polynomial

  DIM ar(n + 1), ai(n + 1), rootr(n + 1), rooti(n + 1), value(n + 1)
  DIM rc AS INTEGER
  DIM skala AS INTEGER

  'define ar vector
  ar(0) = .000089248#
  ar(1) = -.04368#
  ar(2) = 2.9948#
  ar(3) = -6.798#
  ar(4) = 1#

  'ai vector is null
  FOR i = 0 TO n
    ai(i) = 0#
  NEXT i

  PRINT "  Polynomial coefficients:"
  FOR i = 0 TO n
    PRINT TAB(4); i; TAB(10); ar(i); TAB(25); ai(i)
  NEXT i

  FOR j = 0 TO 1
    skala = j
    rc = ibauhub(0, skala, n, ar(), ai(), rootr(), rooti(), value())
    IF rc = 0 THEN
      IF skala = 0 THEN
	PRINT "  Roots (without scaling)"
      ELSE
	PRINT "  Roots (with scaling)"
      END IF

      PRINT "    No   Real part    Imaginary part   Function value"

      FOR i = 0 TO n - 1
	PRINT TAB(2); i; "  ";
	PRINT USING F$; rootr(i);
	PRINT "  ";
	PRINT USING F$; rooti(i);
	PRINT "  "; ; value(i)
      NEXT i
    ELSE
      PRINT "  *** Error in bauhube, rc="; rc
    END IF
  NEXT j

  PRINT "  ----------------------------------------------------------"
  PRINT
  INPUT "", RR$

END 'of main program


  SUB chorner (n, iu, ar(), ai(), xr, xi, pr, pi, p1r, p1i, p2r, p2i, rf1)
 '*====================================================================*
 '*                                                                    *
 '*  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   Double  ar(), ai();                                  *
 '*               Real and imaginary parts of the coefficients of the  *
 '*               polynomial with the coefficients a(iu), ..., a(n)    *
 '*      xr,xi    Double  xr, xi;                                      *
 '*               Real and imaginary parts of the point of evaluation  *
 '*                                                                    *
 '*   Ausgabeparameter:                                                *
 '*   ================                                                 *
 '*      pr, pi   Double  pr, pi;                                      *
 '*               Real and imaginary part of the polynomial            *
 '*      p1r, p1i Double  p1r, p1i;                                    *
 '*               Real and imaginary parts of the 1st derivative there *
 '*      p2r, p2i Double  p2r, p2i;                                    *
 '*               Real and imaginary parts of the 2nd derivative       *
 '*      rf1      Double  rf1;                                         *
 '*               Error estimate for the first derivative              *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   Functions used:                                                  *
 '*   ==============                                                   *
 '*     comabs():    modulus of a complex number                       *
 '*====================================================================*
 '*                                                                    *
 '*   Constant used:  EPS                                              *
 '*   =============                                                    *
 '*                                                                    *
 '*====================================================================*
  XMACHEPS = .000000000000001# ' Small real number
  EPS = 64# * XMACHEPS         ' Accuracy for functional value

  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 TO iu STEP -1
    temp = pr                         ' Polynomial val0ue (pr,pi)
    pr = pr * xr - pi * xi + ar(i)
    pi = pi * xr + temp * xi + ai(i)
    IF i = iu THEN GOTO 20                         ' 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)
    IF temp < temp1 THEN temp = temp1

    IF (temp > rf1) THEN
      rf1 = temp
      i1 = j - 1
    END IF

    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  j = j - 1
  NEXT i

20 temp = comabs(xr, xi)

  IF temp <> 0# THEN
    rf1 = rf1 * temp ^ i1 * (1# * i1 + 1#)
  ELSE
    rf1 = comabs(p1r, p1i)
  END IF

  rf1 = rf1 * EPS

  p2r = p2r + p2r
  p2i = p2i + p2i

END SUB 'chorner

 ' Complex absolute value ....................
 FUNCTION comabs (ar, ai)
 '*====================================================================*
 '*                                                                    *
 '*  Complex absolute value of   a                                     *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   Input parameters:                                                *
 '*   ================                                                 *
 '*      ar,ai    Double  ar, ai;                                      *
 '*               Real, imaginary parts of  a                          *
 '*                                                                    *
 '*   Return value :                                                   *
 '*   =============                                                    *
 '*      Absolute value of a                                           *
 '*                                                                    *
 '*   Functions used :   SQR, ABS                                      *
 '*   ==============                                                   *
 '*                                                                    *
 '*====================================================================*
   IF ar = 0# AND ai = 0# THEN
     comabs = 0#
     EXIT FUNCTION
   END IF

   ar1 = ABS(ar)
   ai1 = ABS(ai)

   IF ai1 > ar1 THEN                               'Switch  ai1 and ar1
     t = ai1: ar1 = ai1: ai1 = t
   END IF

   IF ai1 = 0# THEN
     comabs = ar1
   ELSE
     comabs = ar1 * SQR(1# + (ai1 / ar1) ^ 2)
   END IF
END FUNCTION 'comabs

  FUNCTION ibauhub (real0 AS INTEGER, scale AS INTEGER, n, ar(), ai(), rootr(), rooti(), absf())
 '*====================================================================*
 '*                                                                    *
 '*  ibauhub 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   Double;                                              *
 '*               Real and imaginary parts of the polynomial           *
 '*               coefficients (ar(0), ..., ar(n))                     *
 '*                                                                    *
 '*   Output parameters:                                               *
 '*   =================                                                *
 '*      rootr    Double;     (Vector of length  n+1 )                 *
 '*               rootr(0),..,rootr(n-1) are the real parts of the n   *
 '*               roots                                                *
 '*      rooti    Double;     (Vector of length n+1 )                  *
 '*               rooti(0),..,rooti(n-1) are the imaginary parts       *
 '*               of the roots                                         *
 '*      absf     Double;                                              *
 '*               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:                                    *
 '*   ============================                                     *
 '*     ibauroot():  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                     *
 '*                                                                    *
 '*====================================================================*
  DIM res AS INTEGER

  scalefak = 1#
  IF n < 1 THEN       'n is to small!
    ibauhub = 1
    EXIT FUNCTION
  END IF

  IF ar(n) = 0# AND ai(n) = 0# THEN          ' Leading coefficient must
    ibauhub = 2                               ' differ from zero
    EXIT FUNCTION
  END IF

  FOR i = 0 TO n                     ' store given coefficients in root
    rootr(i) = ar(i)
    rooti(i) = ai(i)
    IF i < n THEN absf(i) = 0#
  NEXT i

  scalefak = 1#
  IF scale <> 0 THEN                     ' Scale polynomial, if desired
    CALL scpoly(n, rootr(), rooti(), scalefak)
  END IF

  x0r = 0#                                           ' Starting val0ue
  x0i = 0#

  FOR i = 0 TO n - 1
						 ' compute the ith root
    res = ibauroot(n, i, rootr(), rooti(), x0r, x0i)

    rootr(i) = scalefak * x0r                              ' store root
    rooti(i) = scalefak * x0i

    IF res <> 0 THEN
      ibauhub = 3                           ' Iteration maximum reached?
      EXIT FUNCTION
    END IF

    ' Polynomial val0ue of input polynomial at (rootr(i), rooti(i))

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

    absf(i) = comabs(tempr, tempi)                       ' store error

    CALL polydiv(n, i, rootr(), rooti(), x0r, x0i)       ' reduce degree

    IF real0 <> 0 THEN                         ' New starting val0ue
      x0i = -x0i                               ' depending on real x..
    ELSE
      x0r = 0#
      x0i = 0#
    END IF
  NEXT i

  ibauhub = 0                                          ' normal exit

END FUNCTION 'ibauhub

  FUNCTION ibauroot (n, iu, ar(), ai(), x0r, x0i)
 '*====================================================================*
 '*                                                                    *
 '*  ibauroot 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   Double  ar(0:n), ai(0:n);                            *
 '*               Real and imaginary parts of the coefficients         *
 '*                                                                    *
 '*   Output parameter:                                                *
 '*   ================                                                 *
 '*      x0r,x0i  Double  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                     *
 '*     icomdiv():   Complex division                                  *
 '*     quadsolv():  solves quadratic equations                        *
 '*                                                                    *
 '*   Constants used : ITERMAX,                                        *
 '*   ===============  QR, QI, XMACHEPS, EPS, EPSROOT, BETA            *
 '*                                                                    *
 '*====================================================================*

  DIM rc AS INTEGER
  DIM result0 AS INTEGER
  DIM endit AS INTEGER

  ITERMAX = 100               ' Maximal number of function
			      ' evaluations per root
  XMACHEPS = .000000000000001#' Small real number
  EPSROOT = 31622000#         ' SQRT(XMACHEPS)
  EPS = 64# * XMACHEPS        ' Accuracy for functional value
  BETA = 8# * EPS
  QR = .1#                    ' Real and imaginary parts of the
  QI = .9#                    ' spiralization constant

  result0 = 2: endit = 0: iter = 0: i = 0
  xoldr = 0#: xoldi = 0#: dxr = 0#: dxi = 0#: bdze = 0#

  IF n < 1 THEN
    ibauroot = 3
    EXIT FUNCTION
  END IF

  IF n - iu = 1 THEN                        ' Polynomial of degree 1
    CALL quadsolv(0#, 0#, ar(n), ai(n), ar(n - 1), ai(n - 1), x0r, x0i)
    ibauroot = 0
    EXIT FUNCTION
  END IF

  IF n - iu = 2 THEN                       ' Polynomial of degree 2
    CALL quadsolv(ar(n), ai(n), ar(n - 1), ai(n - 1), ar(n - 2), ai(n - 2), x0r, x0i)
    ibauroot = 0
    EXIT FUNCTION

  END IF

  xnewr = x0r:  xnewi = x0i
  endit = 0

  ' Evaluate polynomial
  CALL chorner(n, iu, ar(), ai(), xnewr, xnewi, pr, pi, p1r, p1i, p2r, p2i, ss)

  iter = iter + 1

  abspnew = comabs(pr, pi)
  IF abspnew < EPS THEN
    ibauroot = 0                            ' Starting value is a
    EXIT FUNCTION                           ' good approximation
  END IF

  abspold = abspnew
  dzmin = BETA * (EPSROOT + comabs(xnewr, xnewi))


  WHILE iter < ITERMAX                      ' ibauhuber-Iteration

    absp1new = comabs(p1r, p1i)

    IF abspnew > abspold THEN

      i = 0                                 ' dx = dx * q
      iter = iter + 1
      temp = dxr

      dxr = QR * dxr - QI * dxi
      dxi = QR * dxi + QI * temp
    
    ELSE
      
      dzmax = 1# + comabs(xnewr, xnewi)
      h1 = p1r * p1r - p1i * p1i - pr * p2r + pi * p2i
      h2 = 2# * p1r * p1i - pr * p2i - pi * p2r
       IF absp1new > 10# * ss AND comabs(h1, h2) > 100# * ss * ss THEN
	i = i + 1
	IF i > 2 THEN i = 2
	tempr = pr * p1r - pi * p1i
	tempi = pr * p1i + pi * p1r

	rc = icomdiv(-tempr, -tempi, h1, h2, dxr, dxi)
	
	IF rc <> 0 THEN     'Divide by zero?
	  ibauroot = 1
	  EXIT FUNCTION
	END IF

	IF comabs(dxr, dxi) > dzmax THEN
	  temp = dzmax / comabs(dxr, dxi)          ' Newton step
	  dxr = dxr * temp
	  dxi = dxi * temp
	  i = 0
	END IF

	IF i = 2 AND comabs(dxr, dxi) < dzmin / EPSROOT AND comabs(dxr, dxi) > 0# THEN
	  i = 0                                    ' Extrapolation step
	  rc = icomdiv(xnewr - xoldr, xnewi - xoldi, dxr, dxi, h3, h4)
	  IF (rc <> 0) THEN
	    ibauroot = 1
	    EXIT FUNCTION
	  END IF

	  h3 = h3 + 1#
	  h1 = h3 * h3 - h4 * h4
	  h2 = 2# * h3 * h4
	  rc = icomdiv(dxr, dxi, h1, h2, h3, h4)
	  IF rc <> 0 THEN
	    ibauroot = 1
	    EXIT FUNCTION
	  END IF

	  IF comabs(h3, h4) < 50# * dzmin THEN
	    dxr = dxr + h3
	    dxi = dxi + h4
	  END IF
	END IF

	xoldr = xnewr
	xoldi = xnewi
	abspold = abspnew
      ELSE
	i = 0                               ' Close to a saddle point
	h = dzmax / abspnew
	dxr = h * pr
	dxi = h * pi

	xoldr = xnewr
	xoldi = xnewi
	abspold = abspnew

	WHILE ABS(comabs(u, v) / abspnew - 1#) < EPSROOT

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

	  dxr = dxr + dxr
	  dxi = dxi + dxi                   ' dx = dx * 2

	WEND
      END IF
    END IF

    IF endit <> 0 THEN
      IF comabs(dxr, dxi) < .1# * bdze THEN
	xnewr = xnewr + dxr
	xnewi = xnewi + dxi
      END IF

      result0 = 0
      GOTO 100                                 ' stop iteration
    ELSE
      xnewr = xoldr + dxr
      xnewi = xoldi + dxi
      dzmin = BETA * (EPSROOT + comabs(xnewr, xnewi))
      CALL chorner(n, iu, ar(), ai(), xnewr, xnewi, pr, pi, p1r, p1i, p2r, p2i, ss)
      iter = iter + 1
      abspnew = comabs(pr, pi)

      IF abspnew = 0# THEN
	result0 = 0
	GOTO 100
      END IF

      IF comabs(dxr, dxi) < dzmin OR abspnew < EPS THEN
	endit = 1
	bdze = comabs(dxr, dxi)
      END IF
    END IF
  WEND  'End ibauhuber iteration

100 x0r = xnewr
    x0i = xnewi
  ibauroot = result0
  
END FUNCTION 'ibauroot

  FUNCTION icomdiv (ar, ai, br, bi, cr, ci)
 '*====================================================================*
 '*                                                                    *
 '*  Complex division  c = a / b                                       *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   Input parameters:                                                *
 '*   ================                                                 *
 '*      ar,ai    Double  ar, ai;                                      *
 '*               Real, imaginary parts of numerator                   *
 '*      br,bi    Double  br, bi;                                      *
 '*               Real, imaginary parts of denominator                 *
 '*                                                                    *
 '*   Output parameters:                                               *
 '*   ==================                                               *
 '*      cr,ci    Double  cr, ci;                                      *
 '*               Real , imaginary parts of the quotient               *
 '*                                                                    *
 '*   Return value :                                                   *
 '*   =============                                                    *
 '*      = 0      ok                                                   *
 '*      = 1      division by 0                                        *
 '*                                                                    *
 '*====================================================================*
  IF br = 0# AND bi = 0# THEN
    icomdiv = 1
    EXIT FUNCTION
  END IF

  IF ABS(br) > ABS(bi) THEN
    tmp = bi / br
    br = tmp * bi + br
    cr = (ar + tmp * ai) / br
    ci = (ai - tmp * ar) / br
  ELSE
    tmp = br / bi
    bi = tmp * br + bi
    cr = (tmp * ar + ai) / bi
    ci = (tmp * ai - ar) / bi
  END IF

  icomdiv = 0

END FUNCTION 'icomdiv

  SUB polydiv (n, iu, ar(), ai(), x0r, x0i)
 '*====================================================================*
 '*                                                                    *
 '*  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   Double  ar(), ai();                                  *
 '*               Real and imaginary parts of the coefficienten of P   *
 '*               of degree n-iu, with a(iu), ..., a(n)                *
 '*      x0r,x0i  Double  x0r, x0i;                                    *
 '*               Real and imaginary parts of the root x0              *
 '*                                                                    *
 '*   Output parameter:                                                *
 '*   ================                                                 *
 '*      ar, ai   Double  ar(), ai();                                  *
 '*               Real and imaginary parts of the coefficients         *
 '*               ar(iu+1),..,ar(n) of the remainder polynomial Q      *
 '*                                                                    *
 '*====================================================================*
   FOR i = n - 1 TO iu + 1 STEP -1
     temp = ar(i + 1)
     ar(i) = ar(i) + temp * x0r - ai(i + 1) * x0i
     ai(i) = ai(i) + ai(i + 1) * x0r + temp * x0i
   NEXT i

 END SUB 'polydiv

SUB quadsolv (ar, ai, br, bi, cr, ci, tr, ti)
 '*====================================================================*
 '*                                                                    *
 '*  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             Double  ar, ai;      *
 '*      br, bi   coefficient of t                Double  br, bi;      *
 '*      cr, ci   constant term                   Double  cr, ci;      *
 '*                                                                    *
 '*  Output parameter:                                                 *
 '*  ================                                                  *
 '*      tr, ti   complex solution of minimal magnitude                *
 '*                                               Double  tr, ti;      *
 '*                                                                    *
 '*====================================================================*

  pr = br * br - bi * bi
  pi = 2# * br * bi                         '  p = b * b

  qr0 = ar * cr - ai * ci
  qi0 = ar * ci + ai * cr                   '  q = a * c

  pr = pr - 4# * qr0
  pi = pi - 4# * qi0                        ' p = b * b - 4 * a * c

  h = SQR(pr * pr + pi * pi)                ' q = sqrt (p)

  qr0 = h + pr
  IF qr0 > 0# THEN
    qr0 = SQR(qr0 * .5#)
  ELSE
    qr0 = 0#
  END IF

  qi0 = h - pr
  IF qi0 > 0# THEN
    qi0 = SQR(qi0 * .5#)
  ELSE
    qi0 = 0#
  END IF

  IF pi < 0# THEN qi0 = -qi0

  h = qr0 * br + qi0 * bi      ' p = -b +/- q, choose sign for large}
			       ' magnitude  p
  IF h > 0# THEN
    qr0 = -qr0
    qi0 = -qi0
  END IF

  pr = qr0 - br
  pi = qi0 - bi
  h = pr * pr + pi * pi                       ' t = (2 * c) / p

  IF h = 0# THEN
    tr = 0#
    ti = 0#
  ELSE
    tr = 2# * (cr * pr + ci * pi) / h
    ti = 2# * (ci * pr - cr * pi) / h
  END IF

END SUB 'quadsolv

  SUB scpoly (n, ar(), ai(), scal)
 '*====================================================================*
 '*                                                                    *
 '*  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   Double;                                              *
 '*               Real and imaginary parts of the coefficients         *
 '*               a(0),..,a(n)                                         *
 '*                                                                    *
 '*   Output parameters:                                               *
 '*   =================                                                *
 '*      ar, ai   Double;                                              *
 '*               Real and imaginary parts of the coefficients         *
 '*               a(0),..,a(n) of the scaled polynomial.               *
 '*      scal     Double                                               *
 '*               Scaling factor                                       *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   Functions used:                                                  *
 '*   ==============                                                   *
 '*      comabs:  Modulus of a complex number                          *
 '*      Max   :  Maximum of two real numbers                          *
 '*                                                                    *
 '*====================================================================*

  scal = 0#
			      ' scal =
  p = comabs(ar(n), ai(n))    '              a(i)  1/(n-i)
  FOR i = 0 TO n - 1          '   max{ cabs( ---- )       , i=0..n-1
			      '              a(n)
    IF ar(i) <> 0# OR ai(i) <> 0# THEN

      ai(i) = ai(i) / p
      ar(i) = ar(i) / p

      pot = comabs(ar(i), ai(i)) ^ (1# / (n - i))
      IF scal < pot THEN scal = pot
	
    END IF
  NEXT i

  ar(n) = ar(n) / p                   ' Absolute value of a(n) = 1
  ai(n) = ai(n) / p

  IF scal = 0# THEN scal = 1#

  p = 1#
  FOR i = n - 1 TO 0 STEP -1
    p = p * scal              '                    n-i
    ar(i) = ar(i) / p         ' a(i) = a(i) / (scal    ), i=0..n-1
    ai(i) = ai(i) / p         '
  NEXT i
END SUB 'scpoly