Conference rusure::math

    I can send you a couple of routines that will help... for real
    polynomials, I've found it best to put the polynomial into a
    Frobenius matrix (also called a companion matrix) and then use
    an efficient unsymmetric eigenvalue solver to solve for the roots.

    An alternative is to use a procedure written up by Jenkens and
    Traub for factoring polynomials directly.  Interestingly, that
    routine is really a reformulation of a matrix eigenvalue solver; it
    isn't really any faster than the matrix routine, but doesn't
    require space for the matrix.  The speed seems like a contradiction, but
    what happens with the eigenvalue solver, is that it is improving
    all the eigenvalues at once while it finds the next one it's working
    on, so by the time its iterated the first few, the remaining ones
    are already pretty close.  The polynomial rootfinder deflates the
    polynomial by each root and has to start searching from scratch for
    each one.

    Either routine could be coded in basic (yuk) with little trouble.

    I've written other rootfinders, but these are the most reliable
    general purpose methods I've found.

    - Jim

    I have a hard copy of the FORTRAN sources to the Jenkins/Traub 
    polynomial root finders.  I think the one for complex coefficients is
    available in the ACM library of numerical algorithms and the FORTRAN
    for real coefficients is available in some magazine article.  I might
    be able to find the references if someone is interested.
    
    I would like to have a soft copy of these fortran routines.  I am not
    eager to type this in by hand (especially since FORTRAN does not
    require you to declare your variables its very easy to make typos).
    
    Also, the package MAPLE has interesting algorithm for finding
    polynomial roots:  if you give it rational coffecients (ie,
    coefficients that can be expressed exactly either by symbolic means or
    specifying an integer numerator and denominator) it will give you the
    roots of the polynomial, *exactly*.
    
    Does anyone have any references on this algorithm?
    
                                Sieg
    
    

    I have Jenkins' PHD thesis, as well as some papers from the Numer.
    Math and SIAM journal on numerical analysis.  The thesis is useful
    because it contains ALGOL procedures for both the complex and real
    polynomials, as well as a routine to put a posteriori error bounds
    on the roots knowing the finite precision arithmetic being used.

    References are

	J. F. Traub

	A Class of Globally Convergent Iteration Functions for the
	Solution of Polynomial Equations, Math. Comp. V 20, 1966, pp 113-138

	and

	The Calculations of Zeroes of Polynomials and Analytic
	Functions,  AMS Proceedings Symposium of Applied Mathematics,
	Vol 19, 1967 pp 138-152

	[ Provides background on the iteration functions ]

	M. A. Jenkins, J. F. Traub

	A Three Stage Variable-Shift Iteration for Polynomial
	Zeroes and its Relation to Generalized Rayleigh Iteration,
	Numer. Math. 14, 1970, pp 252-263

	[ Describes the complex algorithm ]

	M. A. Jenkins, J. F. Traub

	A Three Stage algorithm for Real Polynomials using Quadratic
	Iterations, SIAM J. Numer Anal, V 7, 1970 pp 545-566

	[ Describes real algorithm ]

    - Jim

C
C  Include file for CPZERO, Jenkins-Traub complex polynomial rootfinder
C
	INTEGER*4 MAXDEG

	PARAMETER (MAXDEG = 99)

	INTEGER*4 NN

	REAL*8 SR, SI, TR, TI, PVR, PVI, ARE, MRE, ETA, INFIN,
	1 PR(MAXDEG), PI(MAXDEG), HR(MAXDEG), HI(MAXDEG),
	2 QPR(MAXDEG), QPI(MAXDEG), QHR(MAXDEG), QHI(MAXDEG),
	3 SHR(MAXDEG), SHI(MAXDEG)

	COMMON /CPZERO/ PR, PI, HR, HI, QPR, QPI, QHR, QHI, SHR, SHI,
	1  SR, SI, TR, TI, PVR, PVI, ARE, MRE, ETA, INFIN, NN




	SUBROUTINE CPZERO(OPR, OPI, NDEG, ZEROR, ZEROI, FAIL)
C
C  ACM Algorithm 419
C
C  Find all zeros of a complex polynomial P(Z) by 3 stage variable
C  shift Jenkins-Traub iteration.  Zeroes are found roughly in
C  order of increasing modulus, and the polynomial is deflated after
C  each zero is found.
C
C  OPR(K),
C  OPI(K) =     Real and imaginary components of coefficients in
C		order of decreasing powers of Z
C
C  NDEG =       Degree of polynomial
C
C  ZEROR(K),
C  ZEROI(K) =	Real and imaginary components of zeroes, roughly in
C		increasing order of modulus
C
C  FAIL	=       TRUE if major and minor passes failed to converge
C
	IMPLICIT NONE

	INCLUDE 'CPZERO.INC'

	INTEGER*4 NDEG, CNT1, CNT2, I
	REAL*8 OPR(1), OPI(1), ZEROR(1), ZEROI(1)
	REAL*8 XX, YY, COSR, SINR, SMALNO, BASE, TT, ZR, ZI, BND,
	1 CMOD, SCALE, CAUCHY, DSQRT
	LOGICAL FAIL, CONV
C
C  Initialize machine specific constants
C
	CALL MCON(ETA, INFIN, SMALNO, BASE)

	ARE = ETA
	MRE = 2.0*DSQRT(2.0D0)*ETA

	XX = 0.70710678
	YY = -XX
	COSR = -0.069756474
	SINR =  0.99756405
	FAIL = .FALSE.
	NN = NDEG+1
C
C  Algorithm fails if leading coefficient is zero
C
	IF (OPR(1).NE.0.0.OR.OPI(1).NE.0.0) GOTO 10
	FAIL = .TRUE.
	RETURN
C
C  Remove zeros at the origin if any
C
10	IF (OPR(NN).NE.0.0.OR.OPI(NN).NE.0.0) GOTO 20
	ZEROR(NDEG-NN+2) = 0.0
	ZEROI(NDEG-NN+2) = 0.0
	NN = NN-1
	GOTO 10
C
C  Make copy of the coefficients
C
20	DO 30 I = 1,NN
	PR(I) = OPR(I)
	PI(I) = OPI(I)
	SHR(I) = CMOD(PR(I),PI(I))
30	CONTINUE
C
C  Scale the polynomial
C
	BND = SCALE(NN, SHR, ETA, INFIN, SMALNO, BASE)
	IF (BND.EQ.1.0) GOTO 40
	DO 35 I = 1,NN
	PR(I) = BND*PR(I)
	PI(I) = BND*PI(I)
35	CONTINUE
C
C  Start the algorithm for one zero
C
40	IF (NN.GT.2) GOTO 50
	IF (NN.EQ.1) RETURN
C
C  Calculate final zero and return
C
	CALL CDIVID(-PR(2), -PI(2), PR(1), PI(1),
	1 ZEROR(NDEG), ZEROI(NDEG))
	RETURN
C
C  Calculate BND, a lower bound on the modulus of the zeroes using
C  Cauchy's estimate
C
50	DO 60 I = 1,NN
	SHR(I) = CMOD(PR(I), PI(I))
60	CONTINUE
	BND = CAUCHY(NN, SHR, SHI)
C
C  Outer loop to control two major passes with different sequences of shifts
C
	DO 100 CNT1 = 1,2
C
C  First stage calculation with 5 no shift iterations
C
	CALL NOSHFT(5)
C
C  Inner loop to select a shift
C
	DO 90 CNT2 = 1,9
C
C  Shift is chosen with modulus BND and amplitude rotated by 94 degrees
C  from previous shift
C
	TT = COSR*XX-SINR*YY
	YY = SINR*XX+COSR*YY
	XX = TT
	SR = BND*XX
	SI = BND*YY
C
C  Second stage calculation with fixed shift
C
	CALL FXSHFT(10*CNT2, ZR, ZI, CONV)
	IF (.NOT.CONV) GOTO 80
C
C  The second stage jumps directly to the third stage iteration.
C  If success, the zero is stored and the polynomial is deflated.
C
	ZEROR(NDEG-NN+2) = ZR
	ZEROI(NDEG-NN+2) = ZI
	NN = NN-1
	DO 70 I = 1,NN
	PR(I) = QPR(I)
	PI(I) = QPI(I)
70	CONTINUE
	GOTO 40

80	CONTINUE
C
C  If iteration was unsuccessful, another shift is tried
C
90	CONTINUE
C
C  If 9 shifts fail, the outer loop is repeated with another sequence
C  of shifts
C
100	CONTINUE
C
C  Failed on two major passes; return empty handed
C
	FAIL = .TRUE.
	RETURN

	END

	SUBROUTINE NOSHFT(L1)
C
C  Compute derivative of polynomial as initial H polynomial and
C  compute L1 no-shift H polynomials
C
	IMPLICIT NONE

	INCLUDE 'CPZERO.INC'

	INTEGER*4 L1, I, J, JJ, N

	REAL*8 XNI, T1, T2, CMOD, XN

	N = NN-1
	XN = N
	DO 10 I = 1,N
	XNI = NN-I
	HR(I) = XNI*PR(I)/XN
	HI(I) = XNI*PI(I)/XN
10	CONTINUE

	DO 50 JJ = 1,L1
	IF (CMOD(HR(N), HI(N)).LE.ETA*10.0D0*CMOD(PR(N), PI(N)))
	1 GOTO 30

	CALL CDIVID(-PR(NN), -PI(NN), HR(N), HI(N), TR, TI)

	DO 20 I = 1,N-1
	J = NN-I
	T1 = HR(J-1)
	T2 = HI(J-1)
	HR(J) = TR*T1-TI*T2+PR(J)
	HI(J) = TR*T2+TI*T2+PI(J)
20	CONTINUE
	HR(1) = PR(1)
	HI(1) = PI(1)
	GOTO 50
C
C  If the constant term is essentially zero, shift H coefficients
C
30	DO 40 I = 1,N-1
	J = NN-I
	HR(J) = HR(J-1)
	HI(J) = HI(J-1)
40	CONTINUE
	HR(1) = 0.0
	HI(1) = 0.0

50	CONTINUE

	RETURN

	END

	SUBROUTINE FXSHFT(L2, ZR, ZI, CONV)
C
C  Compute L2 fixed-shift H polynomials and test for convergence
C
C  Initiates a variable shift iteration and return with the approximate
C  zero if successful
C
C  L2 =	      Limit of fixed shift steps
C  ZR, ZI =   Approximate zero if CONV is true
C  CONV =     Flag indicating convergence of stage 3 iteration
C
	IMPLICIT NONE

	INCLUDE 'CPZERO.INC'

	REAL*8 ZR, ZI, OTR, OTI, SVSR, SVSI, CMOD
	LOGICAL CONV, TEST, PASSED, BOOL
	INTEGER*4 L2, N, I, J

	N = NN-1
C
C  Evaluate P at S
C
	CALL POLYEV(NN, SR, SI, PR, PI, QPR, QPI, PVR, PVI)

	TEST = .TRUE.
	PASSED = .FALSE.
C
C  Calculate first T = -P(S)/H(S)
C
	CALL CALCT(BOOL)
C
C  Main loop for one second stage step
C
	DO 50 J = 1,L2
	OTR = TR
	OTI = TI
C
C  Compute next H polynomial and new T
C
	CALL NEXTH(BOOL)
	CALL CALCT(BOOL)
	ZR = SR+TR
	ZI = SI+TI
C
C  Tests for convergence unless stage 3 has failed once or this is the
C  last H polynomial
C
	IF (BOOL.OR..NOT.TEST.OR.J.EQ.L2) GOTO 50
	IF (CMOD(TR-OTR, TI-OTI).GE.0.5*CMOD(ZR, ZI)) GOTO 40
	IF (.NOT.PASSED) GOTO 30
C
C  The weak convergence test has been passed twice; start the third
C  stage iteration after saving the current H polynomial and shift
C 
	DO 10 I = 1,N
	SHR(I) = HR(I)
	SHI(I) = HI(I)
10	CONTINUE
	SVSR = SR
	SVSI = SI

	CALL VRSHFT(10, ZR, ZI, CONV)
	IF (CONV) RETURN
C
C  The iteration failed to converge; turn off testing and restore
C  H, S, PV, and T
C
	TEST = .FALSE.
	DO 20 I=1,N
	HR(I) = SHR(I)
	HI(I) = SHI(I)
20	CONTINUE
	SR = SVSR
	SI = SVSI
	CALL POLYEV(NN, SR, SI, PR, PI, QPR, QPI, PVR, PVI)
	CALL CALCT(BOOL)
	GOTO 50

30	PASSED = .TRUE.
	GOTO 50

40	PASSED = .FALSE.

50	CONTINUE
C
C  Attempt an iteration with final H polynomial from second stage
C
	CALL VRSHFT(10, ZR, ZI, CONV)

	RETURN

	END

	SUBROUTINE VRSHFT(L3, ZR, ZI, CONV)
C
C  Carry out third stage iteration
C
C  L3 =	      Limit of steps in stage 3
C  ZR, ZI =   On entry contains initial iterate
C	      if iteration converges, contains final iterate
C  CONV	=     .TRUE. if success
C
	IMPLICIT NONE

	INCLUDE 'CPZERO.INC'

	INTEGER*4 L3, I, J

	LOGICAL CONV, B, BOOL

	REAL*8 ZR, ZI, MP, MS, OMP, RELSTP, R1, R2, CMOD, DSQRT,
	1 ERREV, TP

	CONV = .FALSE.

	B = .FALSE.

	SR = ZR
	SI = ZI
C
C  Main loop for stage 3
C
	DO 60 I = 1,L3
C
C  Evaluate P at S and test for convergence
C
	CALL POLYEV(NN, SR, SI, PR, PI, QPR, QPI, PVR, PVI)
	MP = CMOD(PVR, PVI)
	MS = CMOD(SR, SI)

	IF (MP.GT.20.0*ERREV(NN, QPR, QPI, MS, MP, ARE, MRE))
	1 GOTO 10
C
C  Success - polynomial value is smaller than a bound on the error
C
	CONV = .TRUE.
	ZR = SR
	ZI = SI
	RETURN

10	IF (I.EQ.1) GOTO 40
	IF (B.OR.MP.LT.OMP.OR.RELSTP.GE.0.5) GOTO 30
C
C  Iteration has stalled, probably due to a cluster of zeros
C  Do 5 fixed shift iterations into the cluster to try and force one
C  zero to dominate
C
	TP = RELSTP
	B = .TRUE.
	IF (RELSTP.LT.ETA) TP = ETA
	R1 = DSQRT(TP)
	R2 = SR*(1.0+R1)-SI*R1
	SI = SR*R1+SI*(1.0+R1)
	SR = R2
	CALL POLYEV(NN, SR, SI, PR, PI, QPR, QPI, PVR, PVI)
	DO 20 J = 1,5
	CALL CALCT(BOOL)
	CALL NEXTH(BOOL)
20	CONTINUE

	OMP = INFIN
	GOTO 50
C
C  Return if polynomial value increases significantly
C
30	IF (MP*0.1.GT.OMP) RETURN
40	OMP = MP
C
C  Calculate next iterate
C
50	CALL CALCT(BOOL)
	CALL NEXTH(BOOL)
	CALL CALCT(BOOL)
	IF (BOOL) GOTO 60
	RELSTP = CMOD(TR, TI)/CMOD(SR, SI)
	SR = SR+TR
	SI = SI+TI
60	CONTINUE

	RETURN
	END

	SUBROUTINE CALCT(BOOL)
C
C  Calculate T = -P(S)/H(S)
C
C  BOOL = True if H(S) is essentially zero
C
	IMPLICIT NONE

	INCLUDE 'CPZERO.INC'

	INTEGER*4 N
	REAL*8 HVR, HVI, CMOD
	LOGICAL BOOL

	N = NN-1
C
C  Evaluate H(S)
C
	CALL POLYEV(N, SR, SI, HR, HI, QHR, QHI, HVR, HVI)

	BOOL = CMOD(HVR, HVI).LE.ARE*10.0*CMOD(HR(N), HI(N))
	IF (BOOL) GOTO 10
	CALL CDIVID(-PVR, -PVI, HVR, HVI, TR, TI)

	RETURN

10	TR = 0.0
	TI = 0.0
	RETURN

	END

	SUBROUTINE NEXTH(BOOL)
C
C  Calculate the next shifted H polynomial
C
	IMPLICIT NONE

	INCLUDE 'CPZERO.INC'

	REAL*8 T1, T2

	INTEGER*4 I, J

	LOGICAL BOOL

	IF (BOOL) GOTO 20

	DO 10 J = 2,NN-1
	T1 = QHR(J-1)
	T2 = QHI(J-1)
	HR(J) = TR*T1-TI*T2+QPR(J)
	HI(J) = TR*T2+TI*T1+QPI(J)
10	CONTINUE
	HR(1) = QPR(1)
	HI(1) = QPI(1)

	RETURN
C
C  If H(S) is zero replace H with QH
C
20	DO 30 J = 2,NN-1
	HR(J) = QHR(J-1)
	HI(J) = QHI(J-1)
30	CONTINUE
	HR(1) = 0.0
	HI(1) = 0.0

	RETURN

	END

	SUBROUTINE POLYEV(NN, SR, SI, PR, PI, QR, QI, PVR, PVI)
C
C  Evaluate a polynomial via Horner recurrance placing partial
C  sums in Q and computed value in PV
C
	IMPLICIT NONE

	INTEGER*4 I, NN
	REAL*8 PR(NN), PI(NN), QR(NN), QI(NN), SR, SI, PVR, PVI, PVT

	QR(1) = PR(1)
	QI(1) = PI(1)
	PVR = QR(1)
	PVI = QI(1)
	DO 10 I = 2,NN
	PVT = PVR*SR-PVI*SI+PR(I)
	PVI = PVR*SI+PVI*SR+PI(I)
	PVR = PVT
	QR(I) = PVR
	QI(I) = PVI
10	CONTINUE

	RETURN

	END

	REAL*8 FUNCTION ERREV(NN, QR, QI, MS, MP, ARE, MRE)
C
C  Bounds the error in evaluating the polynomial by Horner recurrance
C
C  QR, QI =	Partial sums
C  MS =		Modulus of S
C  MP =		Modulus of P(S)
C  ARE, MRE     Relative error bounds on complex addition and multiplication
C
	IMPLICIT NONE

	INTEGER*4 I, NN

	REAL*8 QR(NN), QI(NN), MS, MP, ARE, MRE, E, CMOD

	E = CMOD(QR(1), QI(1))*MRE/(ARE+MRE)

	DO 10 I = 1,NN
	E = E*MS+CMOD(QR(I), QI(I))
10	CONTINUE

	ERREV = E*(ARE+MRE)-MP*MRE

	RETURN

	END

	REAL*8 FUNCTION CAUCHY(NN, PT, Q)
C
C  Compute a lower bound on the modulus of the zeros of a polynomial
C  using Cauchy's estimate.
C
C  PT(K) =	  Modulus of K'th coefficient
C
	IMPLICIT NONE

	INTEGER*4 N, NN, I

	REAL*8 Q(NN), PT(NN), X, XM, F, DX, DF, DABS, DEXP, DLOG, XN

	N = NN-1
	XN = N
	PT(NN) = -PT(NN)
C
C  Compute upper estimate of bound
C
	X = DEXP((DLOG(-PT(NN))-DLOG(PT(1)))/XN)
	IF (PT(N).EQ.0.0) GOTO 20
C
C  If Newton step at origin is better, use it
C
	XM = -PT(NN)/PT(N)
	IF (XN.LT.X) X = XM
C
C  Chop the interval [0,X] until P <= 0
C
20	XM = X*0.1
	F = PT(1)
	DO 30 I = 2,NN
	F = F*XM+PT(I)
30	CONTINUE
	IF (F.LE.0.0) GOTO 40
	X = XM
	GOTO 20
40	DX = X
C
C  Do Newton iteration until X converges to two decimal places
C
50	IF (DABS(DX/X).LE.0.005) GOTO 70
	Q(1) = PT(1)
	DO 60 I = 2,NN
	Q(I) = Q(I-1)*X+PT(I)
60	CONTINUE

	F = Q(NN)
	DF = Q(1)
	DO 65 I=2,N
	DF = DF*X+Q(I)
65	CONTINUE

	DX = F/DF
	X = X-DX
	GOTO 50

70	CAUCHY = X

	RETURN

	END

	REAL*8 FUNCTION SCALE(NN, PT, ETA, INFIN, SMALNO, BASE)
C
C  Return a scale factor to multiply the coefficients of the polynomial
C
C  The scaling is done to avoid overflow and to avoid undetected underflow
C  interfering with the convergence criterion.  The factor is a power of
C  the base.
C
C  PT(K) =	  Modulus of K'th coefficient
C  ETA, INFIN,
C  SMALNO, BASE = Constants describing machine arithmetic
C
	IMPLICIT NONE

	INTEGER*4 NN, I, L

	REAL*8 PT(NN), ETA, INFIN, SMALNO, BASE, HI, LO,
	1 MAX, MIN, X, SC, DSQRT, DLOG
C
C  Find largest and smallest moduli of coefficients
C
	HI = DSQRT(INFIN)
	LO = SMALNO/ETA
	MAX = 0.0
	MIN = INFIN
	DO 10 I = 1,NN
	X = PT(I)
	IF (X.GT.MAX) MAX = X
	IF (X.NE.0.0.AND.X.LT.MIN) MIN = X
10	CONTINUE
C
C  Scale only if there are very large or small components
C
	SCALE = 1.0
	IF (MIN.GE.LO.AND.MAX.LE.HI) RETURN

	X = LO/MIN
	IF (X.GT.1.0) GOTO 20
	SC = 1.0/(DSQRT(MAX)*DSQRT(MIN))
	GOTO 30
20	SC = X
	IF (INFIN/SC.GT.MAX) SC = 1.0
30	L = DLOG(SC)/DLOG(BASE)+0.5
	SCALE = BASE**L

	RETURN

	END

	SUBROUTINE CDIVID(AR, AI, BR, BI, CR, CI)
C
C  Complex division C = A/B avoiding overflow
C
	IMPLICIT NONE

	REAL*8 AR, AI, BR, BI, CR, CI, T, R, D, INFIN, DABS

	IF (BR.NE.0.0.OR.BI.NE.0.0) goto 10
C
C  Division by zero, return infinity
C
	CALL MCON(T, INFIN, T, T)
	CR = INFIN
	CI = INFIN
	RETURN

10	IF (DABS(BR).GE.DABS(BI)) GOTO 20
	R = BR/BI
	D = BI+R*BR
	CR = (AR*R+AI)/D
	CI = (AI*R-AR)/D
	RETURN

20	R = BI/BR
	D = BR+R*BI
	CR = (AR+AI*R)/D
	CI = (AI-AR*R)/D
	RETURN

	END

	REAL*8 FUNCTION CMOD(R, I)
C
C  Modulus of complex floating number avoiding overflow
C
	IMPLICIT NONE

	REAL*8 R, I, AR, AI, DABS, DSQRT

	AR = DABS(R)
	AI = DABS(I)

	IF (AR.GE.AI) GOTO 10
	CMOD = AI*DSQRT(1.0+(AR/AI)**2)
	RETURN

10	IF (AR.LE.AI) GOTO 20
	CMOD = AR*DSQRT(1.0+(AI/AR)**2)
	RETURN

20	CMOD = AR*DSQRT(2.0D0)
	RETURN

	END

	SUBROUTINE MCON(ETA, LARGENO, SMALLNO, BASE)
C
C
C  MCON provides machine constants used in various parts of the program
C
C  ETA =      Maximum relative representation error which can
C	      be represented as the smallest positive floating point
C	      number such that 1.0 + ETA > 1.0
C  LARGENO =  Largest positive floating point number
C  SMALLNO =  Smallest positive floating point number
C  BASE =     Base of floating point system used
C
C  Let T be the number of base digits in each floating point number
C  Then ETA is either 0.5*B**(1-T) or B**(1-T) depending on whether
C  rounding or truncation is used.
C
C  Let M be the largest exponent and N be the smallest exponent.
C  Then LARGENO = (1-BASE**(-T))*BASE**M
C  and SMALLNO = BASE**N
C
C  The values below correspond to VAX D floating numbers
C
	IMPLICIT NONE

	REAL*8 ETA, LARGENO, SMALLNO, BASE
	REAL*8 XL, XS
	INTEGER M, N, T
C
C Hexadecimal representations of large and small numbers:
C
C Large FFFFFF7F FFFFFFFF
C Small 00000080 00000000
C
	BASE = 2.0
	T = 56
	M =  127
	N = -128
	ETA = BASE**-T
	LARGENO = BASE**(126)*(1.0-BASE**-T)*BASE
	SMALLNO = BASE**(-126)/BASE/BASE

	RETURN

	END

    I finally cheated and purchased EUREKA, a Borland software package
    for the PC.  So far it has been serving me well.  Thanks for all
    the help, I did review my fortran and may get this program running
    as well to see what it does.
    
    
    					Thanks again...
    
    
    					Dennis