Polynomial Zeros in Quadruple Precision

I host an archive of Alan Miller’s Fortran Software, a large collection of routines from a wide range of fields. My goal is simply to make sure these algorithms remain easily accessible, not to actively maintain or enhance them. Occasionally though, users of these algorithms contact me with improvements that others may find useful.

In particular, N. Shamsundar wrote to me to regarding an improvement to rpoly.f90, a Fortran implementation of the Jenkins-Traub algorithm for finding zeros of polynomials. If you use this algorithm and your compiler supports quadruple precision, then you may be interested in making the changes described in his comments below.

I am writing about a suggested enhancement to rpoly.f90. This routine returns the zeros of a polynomial with real coefficients. Run as is, for the test examples it produces these results:

EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,...,10.
 Real part           Imaginary part
   1.00000000250611       0.00000000000000
   2.00002883224649       0.00000000000000
   2.99991349226559       0.00000000000000
   4.00008651702092       0.00000000000000
   4.99997115596089       0.00000000000000
 EXAMPLE 2. POLYNOMIAL WITH ZEROS 1,2,...,10.
 Real part           Imaginary part
   1.00000000236858       0.00000000000000
   2.00002195897130       0.00000000000000
   3.00046637606344       0.00000000000000
   3.99858534914590       0.00000000000000
   5.01275188243750       0.00000000000000
   5.94435077376986       0.00000000000000
   7.13492931078591       0.00000000000000
   7.85730317747578       0.00000000000000
   9.06783881447278       0.00000000000000
   9.98375235450897       0.00000000000000

 Now try case where 1 is an obvious root
  Real part           Imaginary part
 -0.283274635857482E-05   1.00000686290706
 -0.283274635857482E-05  -1.00000686290706
  0.283307420786060E-05  0.999993137276023
  0.283307420786060E-05 -0.999993137276023
  0.999999999344302       0.00000000000000

For Example 2 the polynomial is (x–1).(x–2)…(x–10), and the results are somewhat inaccurate. I believe that the cause is the lower precision of the IEEE floating point arithmetic that is implemented in modern microprocessors as compared to the BURROUGHS B6700 mentioned in the original F77 program of Jenkins and Traub, http://www.netlib.org/toms/493.

The following minor changes enable the results to be improved considerably if the compiler supports quadruple precision, as do Gfortran, Intel and Absoft.

1. In Subroutine RPOLY, at line–73, change the value assigned to ETA to 1D-20.

2. In Subroutine REALIT, line–500 of the file, change the type of variable pv to REAL*16 or kind=selected_real_kind(30,300), whichever is supported by the compiler. I suggest REAL*16; if the user’s compiler objects to it, the user can change it to suit.

After these modifications, the results become:

 EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,...,10.
 Real part           Imaginary part
   1.00000000000000       0.00000000000000
   2.00000000000000       0.00000000000000
   3.00000000000000       0.00000000000000
   4.00000000000000       0.00000000000000
   5.00000000000000       0.00000000000000
 EXAMPLE 1. POLYNOMIAL WITH ZEROS 1,2,...,10.
 Real part           Imaginary part
   1.00000000000000       0.00000000000000
   2.00000000000000       0.00000000000000
   3.00000000000000       0.00000000000000
   4.00000000000000       0.00000000000000
   5.00000000000013       0.00000000000000
   5.99999999999862       0.00000000000000
   7.00000000000477       0.00000000000000
   7.99999999999255       0.00000000000000
   9.00000000000543       0.00000000000000
   9.99999999999849       0.00000000000000

 Now try case where 1 is an obvious root
  Real part           Imaginary part
 -0.333797545686623E-09   1.00000000007661
 -0.333797545686623E-09  -1.00000000007661
  0.333797544801067E-09  0.999999999923389
  0.333797544801067E-09 -0.999999999923389
   1.00000000000000       0.00000000000000

I should note that there is also a quadruple precision version of this algorithm in the archive, q_rpoly.f90, but that version makes use of software-implemented quadruple precision via the module in quad.f90, which is slower than native quadruple precision.