jacobwilliams / polyroots-fortran Goto Github PK

Modern Fortran library for finding the roots of real and complex polynomial equations

License: Other

Fortran 99.16% Python 0.84%

fortran fortran-package-manager polynomial-roots jenkins-traub polynomials companion-matrix eigenvalues root-finding laguerre-method nonlinear-equations

polyroots-fortran's Introduction

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:

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).

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

polyroots-fortran's People

Contributors

Stargazers

Watchers

Forkers

ivan-pi caguerra cem728

polyroots-fortran's Issues

Build fails on PPC due to missing ieee_arithmetic

[  0%]           polyroots_module.F90
 + mkdir -p build/gfortran_EDF3F128DBFCA123/polyroots-fortran/
 + gfortran -c ././src/polyroots_module.F90  -O3 -funroll-loops -Wimplicit-interface -fPIC -fmax-errors=1 -fcoarray=single -J build/gfortran_EDF3F128DBFCA123 -Ibuild/gfortran_EDF3F128DBFCA123 -o build/gfortran_EDF3F128DBFCA123/polyroots-fortran/src_polyroots_module.F90.o
././src/polyroots_module.F90:24:8:

     use ieee_arithmetic
        1
Fatal Error: Can't open module file 'ieee_arithmetic.mod' for reading at (1): No such file or directory
compilation terminated.
[ 50%]           polyroots_module.F90  done.

././src/polyroots_module.F90:24:8:

     use ieee_arithmetic
        1
Fatal Error: Can't open module file 'ieee_arithmetic.mod' for reading at (1): No such file or directory
compilation terminated.
<ERROR> Compilation failed for object " src_polyroots_module.F90.o "
<ERROR>stopping due to failed compilation

Repository keywords

The following keywords (taken from NIST GAMS) may be appropriate: Jenkins-Traub method, Laguerre method, Nonlinear equations, Polynomial, Polynomial equations

Searching for "Jenkins-Traub" gives mostly variants of the original RPOLY routine in C++, Python, Javascript, Haskell, C, Ada, and C#.

Searching for "laguerre-method", I've found the following relevant results:

FPML (a newer library in Fortran by @trcameron; the author also provides FPML-Comp, a library for polishing the polynomial roots)
polynomial-solver (in C++) and rootpolynomialfinder (in C) which implement Laguerre's method, but look more like homework projects

The other keywords provide much more broad and varied results related also to fitting, evaluation, and other types of polynomial equations I never knew existed...

Tests for polynomial root finders

I have extracted the polynomials from an early issue of ACM TOMS (in case anyone is interested in the Algol code, a copy is available here):

Grau, A. A. (1960). Solution of polynomial equation by Bairstow-Hitchcock method. Communications of the ACM, 3(2), 74-75. https://doi.org/10.1145/366959.366966

I've attached the polynomials in the MPSolve polynomial file format: polyBH.tar.gz. MPSolve can be used to find the roots with arbitrary precision. (Some of the polynomials have nice roots, which could be easily specified analytically.)

The ground work for testing polynomial root solvers are laid out in

Jenkins, M. A., & Traub, J. F. (1975). Principles for testing polynomial zero finding programs. ACM Transactions on Mathematical Software (TOMS), 1(1), 26-34. https://doi.org/10.1145/355626.355632

which also provides a collection of test polynomials.

The following resources also provide or refer to collections of test polynomials:

Lang, M., & Frenzel, B. C. (1994). Polynomial root finding. IEEE Signal Processing Letters, 1(10), 141-143. [PDF]
Malek, F. (1995). Polynomial zerofinding matrix algorithms. University of Ottawa (Canada). [PDF]; (I've done my best to transcribe these polynomials from the thesis in the file getpoly.f.txt, however a few support routines are missing and the code is in terrible shape)
Edelman, A., & Murakami, H. (1995). Polynomial roots from companion matrix eigenvalues. Mathematics of Computation, 64(210), 763-776. https://doi.org/10.1090/S0025-5718-1995-1262279-2
Zeng, Z. (2004). Algorithm 835: MultRoot---a Matlab package for computing polynomial roots and multiplicities. ACM Transactions on Mathematical Software (TOMS), 30(2), 218-236. https://doi.org/10.1145/980175.980189

That's about all I can do for one evening.

Make interfaces consistent

Should refactor all the routines so that:

dummy variable names are the same
all accept the coefficients in the same order (currently, some use increasing powers, other decreasing).

one of tests fails

I do not use fpm so I have just compiled the module and the tests using Ubuntu 22.04 default gfortran v. 11.2.0 and liblapack3 (3.10.0) package (for "-llapack"). I used option "-O2" both for module and tests build.

The 'rpoly_test.f90' error-stops in CASE - 2 part of test, at cpqr79:

cpqr79
real part imaginary part root
no convergence in 30 qr iterations. ierr = 4
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
ERROR STOP ** failure in cpqr79 **

Flocke, N. (2015). Algorithm 954: An accurate and efficient cubic and quartic equation solver for physical applications. ACM Transactions on Mathematical Software (TOMS), 41(4), 1-24. https://doi.org/10.1145/2699468

For quartics, an even newer solver in C is available:

Orellana, A. G., & Michele, C. D. (2020). Algorithm 1010: Boosting Efficiency in Solving Quartic Equations with No Compromise in Accuracy. ACM Transactions on Mathematical Software (TOMS), 46(2), 1-28. https://doi.org/10.1145/3386241

The Zip files don't contain license files, so the ACM Software License Agreement applies.