polypp

This project provides four C++ classes (rpoly.hpp, cpoly.hpp, rpolyvp.hpp and cpolyvp.hpp) for solving complex and real polynomials of any degree in multiprecision. The rpolyvp.hpp and cpolyvp.hpp enforce a variable precision (i.e. adaptive) strategy to estimate the roots of a polynomial with a user supplied precision. Quadratic and cubic solver are based on Numerical Recipe book [1]. Quartic solvers are implemented according to Ref. [2]. Roots of polynomials of higher degree are found by using Ehrlich–Aberth algorithm as discussed in Refs. [3-5]. Stopping criterion and accurate calculation of correction term in Aberth method have been implemented as suggested in [6].

References

[1] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes - The Art of Scientific
    Computing (3rd ed.), Cambridge University Press, Cambridge, UK (2007).
[2] A. G. Orellana and C. De Michele, ACM Transactions on Mathematical Software, 46, No. 2, Article 20 (2020),
    doi: https://doi.org/10.1145/3386241.
[3] D. A. Bini, Numerical Algorithms 13, 179-200 (1996).
[4] D. A. Bini and G. Fiorentino, Numerical Algorithms 23, 127–173 (2000).
[5] D. A. Bini and L. Robol, J. Comp. and Appl. Mat. 272, 276-292 (2014).
[6] T. R. Cameron, Numerical Algorithms, 82, 1065–1084 (2019), doi: https://doi.org/10.1007/s11075-018-0641-9
[7] J. L. Aurentz et al., Siam J. Matrix Anal. Appl., 36, 3 942–973 (2015).

In addition to header files you will find some .cpp files with examples on how to use this class, a tool for a statistical analysis of solver accuracy (statanalysis) and another one for testing its performance (timingtest). Multiprecision is implemented through boost multiprecision libraries (https://www.boost.org/doc/libs/1_73_0/libs/multiprecision/doc/html/index.html) and you need to have both boost and gmp (https://gmplib.org/) installed. Boost and gmp are conveniently provided by boost and gmp homebrew packages (https://brew.sh/). Note that homebrew not only supports Mac OSX but also Linux and Windows (see https://docs.brew.sh/Homebrew-on-Linux). To build the sources you need a c++ compiler which complies with c++17. It is strongly recommended to install g++ from homebrew. The package is called gcc and it will provide the g++ executable.

The class itself can be used without boost and gmp and it does not require gcc from homebrew. This means that the poly_real example can be straightforwardly compiled by the command:

make poly_real

without installing anything else.

If you installed boost, gmp and g++ through homebrew (https://brew.sh), the Makefile should work out of the box for compiling all tools and examples, i.e. by issuing the command:

make all

the following executables will be created:

poly_real: example of usage of rpoly.hpp class for finding the roots of real polynomials.

poly_mp: example of usage of rpoly.hpp to find the roots of real polynomials in multiple precision (using boost).

poly_cmplx: example of usage of cpoly.hpp to find the roots of complex polynomials in multiple precision (using boost).

statanalysis: it performs some statistical analyses inspired by the ones which can be found in Ref. [2]. The syntax is the following:

./statanalysis <trials> <output> <sample> <degree>

where:

trials: it is the number of roots (samples A-E) or coefficients (sample F) to generate

output: every trials save the the probability distribution function P(eps_rel) in the file named PE.dat and the cumulative distribution function F(eps_rel) in the file named cumPE.dat.

sample: it is an integer between 0 and 4 which specifies the sample to generate (cf. Table 4 in Ref. [2] with 0={sample A}, 1={sample B}, 2={sample C}, 3={sample D}, 4={sample E} and 5={sample F}).

degree: degree of the polynomials.

timingtest:

this program performs a timing test of the polynomial solver. The syntax is the following:

time ./timingtest <trials>

where

is the number of quartic equations to solve.

By default a polynomial of degree 1000 with random real coefficients between -0.5 and 0.5 is used. The degree of the polynomial can be changed through the macro NDEG defined in timingtest.cpp.

accuracytest:

this program carries out several accuracy tests which can be found in Ref. [7]. The syntax is the following:

./accuracytest <test_number>

where

<test_number> is a number between 1 and 18 which identifies the test to perform (see the sources to recognize the different tests)

timingtest_vp and accuracytest_vp are the same as timingtest and accuracytest respectively, except that they use the variable precision strategy provided by the cpolyvp.hpp (or rpolyvp.hpp) class to guarantee a given output precision.

Parallelization

Parallel calculation of roots is implemeneted through OpenMP (https://www.openmp.org/). To enable the parallelization gnu gcc must be used and in the Makefile replace

PARALLEL=0

with

PARALLEL=1

then issue the commands:

> make clean
> make

and to execute the timing test in parallel, do:

> export OMP_NUM_THREADS=4
> time ./timingtest 1

where you can play with the number of threads by changing the value of the environment variable OMP_NUM_THREADES and the degree of the polynomial used for the timing test by changing the macro NDEG in the file timingtest.cpp.

For polynomials of degree less than 5 the present solver is much faster than tools like MPSolve (ver. 3.1.8, https://numpi.dm.unipi.it/software/mpsolve), while it has a comparable efficiency for higher degree polynomials. Differently from MPSolve, the present algorithms does not perform (yet) any cluster analysis to boost the convergence in case of multiple/clustered roots,

Copyright © 2025 Cristiano De Michele

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than Cristiano De Michele must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from cristiano.demichele@uniroma1.it.

Citing

If you use polypp in your research, please consider giving proper attribution by citing the following publication:

• A.G. Orellana and C. De Michele, Algorithm 1010: Boosting Efficiency in Solving Quartic Equations with No Compromise in Accuracy, ACM Trans. Math. Softw. 46, 20:1-20:28 (2020) https://doi.org/10.1145/3386241

Enjoy!