polyroots-fortran: Polynomial Roots with Modern Fortran

Description

A modern Fortran library for finding the roots of polynomials.

Methods

Many of the methods are from legacy libraries. They have been extensively modified and refactored into Modern Fortran.

Method name	Polynomial type	Coefficients	Roots	Reference
cpoly	General	complex	complex	Jenkins & Traub (1972)
rpoly	General	real	complex	Jenkins & Traub (1975)
rpzero	General	real	complex	SLATEC
cpzero	General	complex	complex	SLATEC
rpqr79	General	real	complex	SLATEC
cpqr79	General	complex	complex	SLATEC
dqtcrt	Quartic	real	complex	NSWC Library
dcbcrt	Cubic	real	complex	NSWC Library
dqdcrt	Quadratic	real	complex	NSWC Library
quadpl	Quadratic	real	complex	NSWC Library
dpolz	General	real	complex	MATH77 Library
cpolz	General	complex	complex	MATH77 Library
polyroots	General	real	complex	LAPACK
cpolyroots	General	complex	complex	LAPACK
rroots_chebyshev_cubic	Cubic	real	complex	Lebedev (1991)
qr_algeq_solver	General	real	complex	Edelman & Murakami (1995)
polzeros	General	complex	complex	Bini (1996)
cmplx_roots_gen	General	complex	complex	Skowron & Gould (2012)
fpml	General	complex	complex	Cameron (2019)

The library also includes some utility routines:

• newton_root_polish
• sort_roots

Example

An example of using polyroots to compute the zeros for a 5th order real polynomial $$P(x) = x^5 + 2x^4 + 3x^3 + 4x^2 + 5x + 6$$

program example

use iso_fortran_env
use polyroots_module, wp => polyroots_module_rk

implicit none

integer,parameter :: degree = 5 !! polynomial degree
real(wp),dimension(degree+1) :: p = [1,2,3,4,5,6] !! coefficients

integer :: i !! counter
integer :: istatus !! status code
real(wp),dimension(degree) :: zr !! real components of roots
real(wp),dimension(degree) :: zi !! imaginary components of roots

call polyroots(degree, p, zr, zi, istatus)

write(*,'(A,1x,I3)') 'istatus: ', istatus
write(*, '(*(a22,1x))') 'real part', 'imaginary part'
do i = 1, degree
    write(*,'(*(e22.15,1x))') zr(i), zi(i)
end do

end program example

The result is:

istatus:    0
             real part         imaginary part
 0.551685463458982E+00  0.125334886027721E+01
 0.551685463458982E+00 -0.125334886027721E+01
-0.149179798813990E+01  0.000000000000000E+00
-0.805786469389031E+00  0.122290471337441E+01
-0.805786469389031E+00 -0.122290471337441E+01

Compiling

A fpm.toml file is provided for compiling polyroots-fortran with the Fortran Package Manager. For example, to build:

fpm build --profile release

By default, the library is built with double precision (real64) real values. Explicitly specifying the real kind can be done using the following processor flags:

Preprocessor flag	Kind	Number of bytes
REAL32	real(kind=real32)	4
REAL64	real(kind=real64)	8
REAL128	real(kind=real128)	16

For example, to build a single precision version of the library, use:

fpm build --profile release --flag "-DREAL32"

All routines, except for polyroots are available for any of the three real kinds. polyroots is not available for real128 kinds since there is no corresponding LAPACK eigenvalue solver.

To run the unit tests:

fpm test

To use polyroots-fortran within your fpm project, add the following to your fpm.toml file:

[dependencies]
polyroots-fortran = { git="https://github.com/jacobwilliams/polyroots-fortran.git" }

or, to use a specific version:

[dependencies]
polyroots-fortran = { git="https://github.com/jacobwilliams/polyroots-fortran.git", tag = "1.2.0"  }

To generate the documentation using ford, run: ford ford.md

Documentation

The latest API documentation for the master branch can be found here. This was generated from the source code using FORD.

License

The polyroots-fortran source code and related files and documentation are distributed under a permissive free software license (BSD-style).

See also

• Roots-Fortran

Similar libraries in other programming languages

• R: polyroot
• MATLAB: roots
• C: GSL - Polynomials, MPSolve
• Julia: PolynomialRoots.jl, FastPolynomialRoots.jl, Polynomials.jl
• Python: numpy.polynomial.polynomial

Other references and codes

• GAMS Class F1a.
• eiscor - eigensolvers based on unitary core transformations containing the AMVW method from the work of Aurentz et al. (2015), Fast and Backward Stable Computation of Roots of Polynomials (an earlier version can be picked up from the website of Ran Vandebril, one of the co-authors of that paper).
• PA03 HSL Archive code for computing all the roots of a cubic polynomial
• PA05 HSL Archive code for computing all the roots of a quartic polynomial
• PA16, PA17 HSL codes for computing zeros of polynomials using method of Madsen and Reid
• Various codes from Alan Miller
• A solver using the companion matrix and LAPACK
• Root-finding algorithms: Roots of Polynomials | Wikipedia
• Polynomial Roots | Wolfram MathWorld
• What is a Companion Matrix | Nick Higham
• 19 Dubious Ways to Compute the Zeros of a Polynomial | Cleve's Corner
• New Progress in Polynomial Root-finding | Victor Y. Pan

See also

• Code coverage statistics [codecov.io]