ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation.
pm_polynomial Module Reference

This module contains procedures and generic interfaces for performing various mathematical operations involving polynomials.
More...

Data Types

type  eigen_type
 This is a concrete derived type whose instances are exclusively used to signify the use of the Eigenvalue method of root-finding.
More...
 
interface  getPolyAdd
 Generate and return the vector of coefficients of the polynomial resulting from the addition of a polynomial to another polynomial of arbitrary degrees. More...
 
interface  getPolyDiff
 Generate and return the vector of coefficients of the polynomial resulting from the \(k\)th-order differentiation of a univariate polynomial of arbitrary degree. More...
 
interface  getPolyMul
 Generate and return the vector of coefficients of the polynomial resulting from the multiplication of a polynomial with another polynomial of arbitrary degrees. More...
 
interface  getPolyRoot
 Generate and return the roots of a polynomial of arbitrary degree specified by its coefficients coef.
More...
 
interface  getPolyStr
 Generate and return a string containing the polynomial expression corresponding to the input polynomial coefficients. More...
 
interface  getPolySub
 Generate and return the vector of coefficients of the polynomial resulting from the subtraction of a polynomial to another polynomial of arbitrary degrees. More...
 
interface  getPolyVal
 Generate and return the value of the polynomial of arbitrary degree whose coefficients are specified by the user in the order of increasing power.
More...
 
type  jenkins_type
 This is a concrete derived type whose instances are exclusively used to signify the use of Jenkins-Traub method of root-finding.
More...
 
type  laguerre_type
 This is a concrete derived type whose instances are exclusively used to signify the use of Laguerre method of root-finding.
More...
 
type  method_type
 This is an abstract derived type for constructing concrete derived types to distinguish various procedure signatures that require root-finding methods (e.g., Eigenvalue, Jenkins, Laguerre, ...).
More...
 
interface  setPolyAdd
 Return the vector of coefficients of the polynomial resulting from the addition of a polynomial to another polynomial of arbitrary degrees. More...
 
interface  setPolyDiff
 Return the vector of coefficients of the polynomial resulting from the \(k\)th-order differentiation of a univariate polynomial of arbitrary degree. More...
 
interface  setPolyDiv
 Return the quotient and remainder of dividing a polynomial with another polynomial of arbitrary degrees. More...
 
interface  setPolyMul
 Return the vector of coefficients of the polynomial resulting from the multiplication of a polynomial with another polynomial of arbitrary degrees. More...
 
interface  setPolyRoot
 Return the roots of a polynomial of arbitrary degree specified by its coefficients coef.
More...
 
interface  setPolyRootPolished
 Return the polished (refined) root of a polynomial of arbitrary degree specified by its coefficients coef.
More...
 
interface  setPolySub
 Return the vector of coefficients of the polynomial resulting from the subtraction of a polynomial to another polynomial of arbitrary degrees. More...
 
type  sgl_type
 This is a concrete derived type whose instances are exclusively used to signify the use of Skowron-Gould method of root-finding.
More...
 

Variables

character(*, SK), parameter MODULE_NAME = "@pm_polynomial"
 
type(sgl_type), parameter sgl = sgl_type()
 This is a scalar parameter object of type sgl_type that is exclusively used to signify the use of Skowron-Gould method of root-finding within an interface of a procedure of the ParaMonte library.
More...
 
type(eigen_type), parameter eigen = eigen_type()
 This is a scalar parameter object of type eigen_type that is exclusively used to signify the use of Eigenvalue method of root-finding within an interface of a procedure of the ParaMonte library.
More...
 
type(jenkins_type), parameter jenkins = jenkins_type()
 This is a scalar parameter object of type jenkins_type that is exclusively used to signify the use of Jenkins-Traub method of root-finding within an interface of a procedure of the ParaMonte library.
More...
 
type(laguerre_type), parameter laguerre = laguerre_type()
 This is a scalar parameter object of type laguerre_type that is exclusively used to signify the use of Laguerre method of root-finding within an interface of a procedure of the ParaMonte library.
More...
 

Detailed Description

This module contains procedures and generic interfaces for performing various mathematical operations involving polynomials.

Specifically, this module contains generic interfaces for the following computational polynomial tasks:

  1. Evaluation of polynomial functions using the Horner method.
  2. Evaluation of the roots of polynomials of arbitrary degrees.
  3. Evaluation of the addition of two polynomials of arbitrary degrees.
  4. Evaluation of the subtraction of two polynomials of arbitrary degrees.
  5. Evaluation of the differentiation of polynomials of arbitrary degrees.
  6. converting univariate polynomial coefficients to readable univariate polynomial expressions in string format for better display.

Introduction

In mathematics, a polynomial is an expression consisting of indeterminates (also called variables) and coefficients, that involves only the operations of addition, subtraction, multiplication, and positive-integer powers of variables.
An example of a polynomial of a single indeterminate \(x\) is \(x^2 − 4x + 7\).
An example with three indeterminates is \(x^3 + 2xyz2 − yz + 1\).
Polynomials appear in many areas of mathematics and science.
For example, they are used to form polynomial equations, which encode a wide range of problems.
They are used to define polynomial functions, which appear in settings ranging from basic chemistry and physics to economics and social science.
They are used in calculus and numerical analysis to approximate other functions.
In advanced mathematics, polynomials are used to construct polynomial rings and algebraic varieties, which are central concepts in algebra and algebraic geometry.

Definition

A polynomial expression is an expression that can be built from constants and symbols called variables or indeterminates by means of addition, multiplication and exponentiation to a non-negative integer power.
The constants are generally numbers, but may be any expression that do not involve the indeterminates, and represent mathematical objects that can be added and multiplied.
Two polynomial expressions are considered as defining the same polynomial if they may be transformed, one to the other, by applying the usual properties of commutativity, associativity and distributivity of addition and multiplication.
For example \((x-1)(x-2)\) and \(x^{2}-3x+2\) are two polynomial expressions that represent the same polynomial.
Therefore, one has the equality \((x-1)(x-2) = x^{2}-3x+2\).
A polynomial in a single indeterminate x can always be written (or rewritten) in the form

\begin{equation} a_{n}x^{n}+a_{n-1}x^{n-1}+\dotsb +a_{2}x^{2}+a_{1}x+a_{0} ~, \end{equation}

where \(a_{0}, \ldots , a_{n}\) are constants that are called the coefficients of the polynomial, and \(x\) is the indeterminate.
The word indeterminate means that \(x\) represents no particular value, although any value may be substituted for it.
The mapping that associates the result of this substitution to the substituted value is a function, called a polynomial function.
This can be expressed more concisely by using summation notation:

\begin{equation} \sum_{k = 0}^{n} a_{k}x^{k} ~. \end{equation}

In other words, a polynomial can either be zero or can be written as the sum of a finite number of non-zero terms.
Each term consists of the product of a number - called the coefficient of the term – and a finite number of indeterminates, raised to non-negative integer powers.

Addition

Polynomials can be added using the associative law of addition (grouping all their terms together into a single sum), possibly followed by reordering (using the commutative law) and combining of like terms.
For example, if \(P = 3x^{2} - 2x + 5xy - 2\) and \(Q = -3x^{2} + 3x + 4y^{2} + 8\), then the sum,

\begin{equation} P + Q = 3x^{2} - 2x + 5xy - 2 - 3x^{2} + 3x + 4y^{2} + 8 ~, \end{equation}

can be reordered and regrouped as,

\begin{equation} P + Q = (3x^{2} - 3x^{2}) + (-2x + 3x) + 5xy + 4y^{2} + (8 - 2) ~, \end{equation}

and then simplified to,

\begin{equation} P + Q = x + 5xy + 4y^{2} + 6 ~. \end{equation}

When polynomials are added together, the result is another polynomial.
In summary, to add two polynomials, simply add the corresponding coefficients of terms of equal exponent in the two polynomials.
This results in a new set of polynomial coefficients which represent the polynomial resulting from the sum of the two polynomials.

Note
The degree of the resulting polynomial from the addition of two other polynomials is always the maximum of the degrees of the two polynomials added.

Subtraction

Subtraction of polynomials is similar to their additions.
To subtract two polynomials, simply subtract the corresponding coefficients of terms of equal exponent in the two polynomials.
This results in a new set of polynomial coefficients which represent the polynomial resulting from the subtraction of the two polynomials.

Note
The degree of the resulting polynomial from the subtraction of two other polynomials is always the maximum of the degrees of the two polynomials subtracted.

Multiplication

Polynomials can also be multiplied.
To expand the product of two polynomials into a sum of terms, the distributive law is repeatedly applied, which results in each term of one polynomial being multiplied by every term of the other.
For example, if

\begin{equation} \begin{aligned} \color{Red} P & \color{Red} {= 2x + 3y + 5} \\ \color{Blue} Q & \color{Blue} {= 2x + 5y + xy + 1} \end{aligned} \end{equation}

then,

\begin{equation} \begin{array}{rccrcrcrcr} {\color {Red}{P}}{\color {Blue}{Q}}&{=}&&({\color {Red}{2x}}\cdot {\color {Blue}{2x}})& + &({\color {Red}{2x}}\cdot {\color {Blue}{5y}})& + &({\color {Red}{2x}}\cdot {\color {Blue}{xy}})& + &({\color {Red}{2x}}\cdot {\color {Blue}{1}}) \\ && + &({\color {Red}{3y}}\cdot {\color {Blue}{2x}})&+&({\color {Red}{3y}}\cdot {\color {Blue}{5y}})&+&({\color {Red}{3y}}\cdot {\color {Blue}{xy}})&+&({\color {Red}{3y}}\cdot {\color {Blue}{1}}) \\ && + &({\color {Red}{5}}\cdot {\color {Blue}{2x}})&+&({\color {Red}{5}}\cdot {\color {Blue}{5y}})&+&({\color {Red}{5}}\cdot {\color {Blue}{xy}})&+&({\color {Red}{5}}\cdot {\color {Blue}{1}}) \end{array} \end{equation}

Carrying out the multiplication in each term produces,

\begin{equation} \begin{array}{rccrcrcrcr} PQ &=&& 4x^{2}&+&10xy&+&2x^{2}y&+&2x \\ &&+&6xy&+&15y^{2}&+&3xy^{2}&+&3y\\ &&+&10x&+&25y&+&5xy&+&5. \end{array} \end{equation}

Combining similar terms yields,

\begin{equation} \begin{array}{rcccrcrcrcr} PQ &=&& 4x^{2}&+&(10xy+6xy+5xy)&+&2x^{2}y&+&(2x+10x) \\ &&+&15y^{2}&+&3xy^{2}&+&(3y+25y)&+&5 \end{array} \end{equation}

which can be simplified to,

\begin{equation} PQ = 4x^{2}+21xy+2x^{2}y+12x+15y^{2}+3xy^{2}+28y+5 ~. \end{equation}

As in the example, the product of polynomials is always a polynomial.

In summary, To expand the product of two polynomials into a sum of terms, the distributive law is repeatedly applied.
This results in each term of one polynomial being multiplied by every term of the other.

Note
The degree of the resulting polynomial from the multiplication of two other polynomials is always the sum of the degrees of the two multiplicands plus one.

Division

The division of one polynomial by another is not typically a polynomial.
Instead, such ratios are a more general family of objects, called rational fractions, rational expressions, or rational functions, depending on context.
This is analogous to the fact that the ratio of two integers is a rational number, not necessarily an integer.
For example, the fraction \(1 / (x^2 + 1)\) is not a polynomial, and it cannot be written as a finite sum of powers of the variable \(x\).
For polynomials in one variable, there is a notion of Euclidean division of polynomials, generalizing the Euclidean division of integers.
This notion of the division \(a(x) / b(x)\) results in two polynomials, a quotient \(q(x)\) and a remainder \(r(x)\), such that \(a = b q + r\) and \(\ms{degree}(r) < \ms{degree}(b)\).
The quotient and remainder may be computed by any of several algorithms, including polynomial long division and synthetic division.
When the denominator \(b(x)\) is monic and linear, that is, \(b(x) = x − c\) for some constant \(c\), then the polynomial remainder theorem asserts that the remainder of the division of \(a(x)\) by \(b(x)\) is the evaluation \(a(c)\).
In this case, the quotient may be computed by the Ruffini rule, a special case of synthetic division.

Division algorithm

The methodology employed in this module relies on polynomial long division.
Polynomial long division is an algorithm for dividing a polynomial by another polynomial of the same or lower degree.
It is a generalized version of the familiar arithmetic technique called long division.
Another abbreviated method is polynomial short division (Blomqvist method).
The method of polynomial long division implements the Euclidean division of polynomials.
Starting from the polynomial \(\ms{dividend}\) and a non-zero polynomial \(\ms{divisor}\), it outputs a \(\ms{Quotient}\) polynomial and a \(\ms{Remainder}\) polynomial such that,

\begin{equation} \ms{dividend} = \ms{divisor} \times \ms{Quotient} + \ms{Remainder} ~. \end{equation}

By definition, the resulting degree of \(\ms{Remainder}\) is lower than the degree of \(\ms{divisor}\).
The resulting polynomials \(\ms{Quotient}\) and \(\ms{Remainder}\) are unique and do not depend on the derivation methodology.
The result \(\ms{Remainder} = 0\) occurs if and only if the polynomial \(\ms{dividend}\) has \(\ms{divisor}\) as a factor.

Note
The polynomial division is a means for testing whether one polynomial has another as a factor, and, if it does, for factoring it out.
For example, if \(r\) is a root of the polynomial \(\ms{dividend}\) is known, it can be factored out by dividing \(\ms{dividend}\) by \(\ms{divisor} = (x – r)\).

Polynomial long division

Derivative

Calculating derivatives and integrals of polynomials is particularly simple, compared to other kinds of functions.
The derivative of the polynomial,

\begin{equation} P = a_{n}x^{n}+a_{n-1}x^{n-1}+\dots +a_{2}x^{2}+a_{1}x+a_{0}=\sum _{i=0}^{n}a_{i}x^{i} ~, \end{equation}

with respect to \(x\) is the polynomial,

\begin{equation} na_{n}x^{n-1}+(n-1)a_{n-1}x^{n-2}+\dots +2a_{2}x+a_{1}=\sum _{i=1}^{n}ia_{i}x^{i-1} ~. \end{equation}

Succinctly, given a polynomial of degree \(\ms{degree}\),

\begin{equation} P(x) = \sum_{i ~=~ 0}^{i ~=~ \ms{degree}} ~ c_i x^i ~, \end{equation}

where \(c_i\) is the \(i\)th coefficient of the polynomial, the \(k\)th-order derivative of the polynomial can be computed as,

\begin{equation} \frac{\mathrm{d}^k P(x)}{\mathrm{d}x^k} = \begin{cases} \sum_{i ~=~ k}^{i ~=~ \ms{degree}} ~ \big( \prod_{j = i - k + 1}^{j = i} \big) c_i x^{(i - k)} ~,~& k < \ms{degree} \\ 0 ~,~& k \geq \ms{degree} \end{cases} \end{equation}

Note
The degree of the resulting polynomial from the \(k\)th-order differentiation of a polynomial of degree \(\ms{degree}\) is by definition \(\max(0, \ms{degree} - k)\).

Integral

Similarly, the general antiderivative (or indefinite integral) of \(P\) is,

\begin{equation} {\frac {a_{n}x^{n+1}}{n+1}}+{\frac {a_{n-1}x^{n}}{n}}+\dots +{\frac {a_{2}x^{3}}{3}}+{\frac {a_{1}x^{2}}{2}}+a_{0}x+c=c+\sum _{i=0}^{n}{\frac {a_{i}x^{i+1}}{i+1}} ~, \end{equation}

where \(c\) is an arbitrary constant.
For example, antiderivatives of \(x^2 + 1\) have the form \(\frac{1}{3}x^3 + x + c\).

Polynomial functions and root-finding

A polynomial equation, also called an algebraic equation, is an equation of the form,

\begin{equation} a_{n}x^{n}+a_{n-1}x^{n-1}+\dotsb +a_{2}x^{2}+a_{1}x+a_{0}=0 ~. \end{equation}

For example,

\begin{equation} 3x^{2}+4x-5 = 0 ~, \end{equation}

is a polynomial equation.

A root of a nonzero univariate polynomial \(P\) is a value a of \(x\) such that \(P(a) = 0\).
In other words, a root of \(P\) is a solution of the polynomial equation \(P(x) = 0\) or a zero of the polynomial function defined by \(P\).
In the case of the zero polynomial, every number is a zero of the corresponding function, and the concept of root is rarely considered.
A number \(a\) is a root of a polynomial \(P\) if and only if the linear polynomial \(x − a\) divides \(P\), that is if there is another polynomial \(Q\) such that \(P = (x − a) Q\).
It may happen that a power (greater than 1) of \(x − a\) divides \(P\); in this case, \(a\) is a multiple root of \(P\), and otherwise \(a\) is a simple root of \(P\).
If \(P\) is a nonzero polynomial, there is a highest power \(m\) such that \((x − a)m\) divides \(P\), which is called the multiplicity of \(a\) as a root of \(P\).
The number of roots of a nonzero polynomial \(P\), counted with their respective multiplicities, cannot exceed the degree of \(P\), and equals this degree if all complex roots are considered (this is a consequence of the fundamental theorem of algebra).
The coefficients of a polynomial and its roots are related by the Vieta formulas.

Some polynomials, such as \(x^2 + 1\), do not have any roots among the real numbers.
If, however, the set of accepted solutions is expanded to the complex numbers, every non-constant polynomial has at least one root;
This is the fundamental theorem of algebra.
By successively dividing out factors \(x − a\), one sees that any polynomial with complex coefficients can be written as a constant (its leading coefficient) times a product of such polynomial factors of degree \(1\);
Consequently, the number of (complex) roots counted with their multiplicities is exactly equal to the degree of a polynomial.

Polynomial root-finding: Eigenvalue method

Given the real or complex coefficients \(c_i\) of a polynomial of degree \(n > 0\),

\begin{equation} f(x) = \sum_{i = 0}^{n} ~ c_i x^{i} ~, \end{equation}

the procedures of this module compute the \(n\) (potentially complex) roots of the polynomial \(f(x)\).
The roots are always returned as complex numbers as polynomials with real coefficients can also have complex roots.
The algorithms of this module rely on finding the eigenvalues of a matrix \(A\) that are the roots of its characteristic polynomial \(f(x) = \mathrm{det}[A − xI]\).
It can be shown that a monic polynomial \(f(x)\) is the characteristic polynomial of the special \(n\times n\) square Frobenius Companion Matrix,

\begin{equation} C(f) = \begin{bmatrix} 0 & 0 & \dots & 0 & -c_{0} \\ 1 & 0 & \dots & 0 & -c_{1} \\ 0 & 1 & \dots & 0 & -c_{2} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & -c_{n-1} \end{bmatrix} \end{equation}

Although the eigenvalue root-finding method can be slower than other polynomial root-finding methods, it is generally a more robust technique than others.

See also
pm_polynomial
pm_polynomial
MATH77
Root-Finding Algorithms
S. Goedecker, Remark on algorithms to find roots of polynomials, SIAM J. on Scientific Computing 15, 5 (Sept. 1994) 1058-1063.
B. S. Garbow, J. M. Boyle, J. J. Dongarra, and C. B. Moler, Matrix Eigensystem Routines | EISPACK Guide Extension, Lecture Notes in Computer Science 51, Springer Verlag, Berlin (1977) 343 pages.
B. T. Smith, J. M. Boyle, B. S. Garbow, Y. Ikebe, V. C. Klema, and C. B. Moler, Matrix Eigensystem Routines | EISPACK Guide, Lecture Notes in Computer Science 6, Springer Verlag, Berlin (1974) 387 pages.

Polynomial root-finding: Jenkins–Traub method

Given the real or complex coefficients \(c_i\) of a polynomial of degree \(n > 0\),

\begin{equation} f(x) = \sum_{i = 0}^{n} ~ c_i x^{i} ~, \end{equation}

the root-finding procedures of this module also compute the \(n\) (potentially complex) roots of the polynomial \(f(x)\) using the Jenkins–Traub method.
The roots are always returned as complex numbers as polynomials with real coefficients can also have complex roots.
The algorithms of this module rely on a three-stage algorithm for finding roots of polynomials with real or complex coefficients as outlined in Jenkins and Traub (1970).

Note
Although the Jenkins–Traub method is a popular fairly robust method of polynomial root-finding, it can still fail to find roots of certain highly unstable polynomials, some of which are exemplified in the original paper of Jenkins and Traub (1970).
In general, if the polynomial coefficients have highly variable number of significant digits, that is, they are of vastly different magnitudes (some very small and some very large), the Jenkins–Traub method might be prone to failure.
This can happen when, for example, the difference in the number of significant digits of smallest and largest coefficients is comparable to the precision of the real or complex kind used to represent the coefficients.
In such a case, using a higher-precision real or complex kind may resolve the instability.

Which polynomial root-finding method should I use?
If you have a polynomial of highly varying coefficients, then Eigenvalue method is likely going to be the more reliable approach.
The Jenkins-Traub is also considered a relatively reliable fast Sure-Fire technique for finding the roots of polynomials.

Remarks
The algorithms of this module that take polynomial coefficients of type real are expected to be generally up to four times faster than the corresponding algorithms that accept coefficients of type complex.
The extent of the accuracy of above claim is yet to be tested for the specific implementations of this module.
See also
pm_mathRoot
Root-Finding Algorithms
Jenkins and Traub, 1970, A Three-Stage Algorithm for Real Polynomials Using Quadratic Iteration, SIAM.
Jenkins, ALGORITHM 493 - Zeros of a Real Polynomial [C2], ACM Transactions on Mathematical Software, Vol 1, No 2, June 1975.

Polynomial root-finding: Laguerre method

In numerical analysis, the Laguerre method is a root-finding algorithm tailored to polynomials.
In other words, the Laguerre method can be used to numerically solve the equation \(p(x) = 0\) for a given polynomial \(p(x)\).
One of the most useful properties of this method is that it is, from extensive empirical study, very close to being a sure-fire method, meaning that it is almost guaranteed to always converge to some root of the polynomial, no matter what initial guess is chosen.
However, there are more efficient methods with which it is guaranteed to find all roots or all real roots.
This method is named in honor of the French mathematician Edmond Laguerre.

The algorithm of the Laguerre method to find one root of a polynomial \(p(x)\) of degree \(n\) is:

  1. Choose an initial guess \(x_0\).
  2. For \(k = 0, 1, 2, \ldots\),
    1. If \(p(x_{k})\) is very small, exit the loop.
    2. Calculate \(G = {\frac {p'(x_{k})}{p(x_{k})}}\).
    3. Calculate \(H = G^{2} - {\frac {p''(x_{k})}{p(x_{k})}}\).
    4. Calculate \(a = {\frac{n}{G\pm{\sqrt{(n - 1)(nH - G^{2})}}}}\), where the sign is chosen to give the denominator with the larger absolute value, to avoid catastrophic cancellation.
    5. Set \(x_{{k+1}}=x_{k} - a\).
  3. Repeat until a is small enough or if the maximum number of iterations has been reached.

If a root has been found, the corresponding linear factor can be removed from \(p\).
This deflation step reduces the degree of the polynomial by one, so that eventually, approximations for all roots of \(p\) can be found.
Note however that deflation can lead to approximate factors that differ significantly from the corresponding exact factors.
This error is least if the roots are found in the order of increasing magnitude.

Benchmarks:


Benchmark :: The effects of method on runtime efficiency

The following program compares the runtime performance of setPolyRoot using different polynomial root finding algorithms.
1! Test the performance of `setPolyRoot()` with and without the selection `control` argument.
2program benchmark
3
4 use pm_bench, only: bench_type
5 use pm_distUnif, only: setUnifRand
7 use pm_kind, only: SK, IK, LK, RKH, RK, RKG => RKD
8 use iso_fortran_env, only: error_unit
9
10 implicit none
11
12 integer(IK) , parameter :: degmax = 9_IK
13 integer(IK) , allocatable :: rootCount(:)
14 integer(IK) :: ibench
15 integer(IK) :: ideg
16 integer(IK) :: degree
17 integer(IK) :: fileUnit
18 integer(IK) :: rootUnit
19 real(RKG) :: dumsum = 0._RKG
20 complex(RKG) :: coef(0:2**degmax)
21 complex(RKG) :: root(2**degmax)
22 type(bench_type) , allocatable :: bench(:)
23
24 bench = [ bench_type(name = SK_"Eigen", exec = setPolyRootEigen, overhead = setOverhead) &
25 , bench_type(name = SK_"Jenkins", exec = setPolyRootJenkins, overhead = setOverhead) &
26 , bench_type(name = SK_"Laguerre", exec = setPolyRootLaguerre, overhead = setOverhead) &
27 , bench_type(name = SK_"SGL", exec = setPolyRootSGL, overhead = setOverhead) &
28 ]
29 allocate(rootCount(size(bench)))
30
31 write(*,"(*(g0,:,' '))")
32 write(*,"(*(g0,:,' '))") "Benchmarking setPolyRoot()"
33 write(*,"(*(g0,:,' '))")
34
35 open(newunit = fileUnit, file = "main.out", status = "replace")
36 open(newunit = rootUnit, file = "root.count", status = "replace")
37
38 write(fileUnit, "(*(g0,:,','))") "Polynomial Degree", (bench(ibench)%name, ibench = 1, size(bench))
39 write(rootUnit, "(*(g0,:,','))") "Polynomial Degree", (bench(ibench)%name, ibench = 1, size(bench))
40 loopOverArraySize: do ideg = 0, degmax
41
42 degree = 2**ideg
43 call setUnifRand(coef(0 : degree), (1._RKG, 1._RKG), (2._RKG, 2._RKG))
44 write(*,"(*(g0,:,' '))") "Benchmarking with coef size", degree
45 do ibench = 1, size(bench)
46 bench(ibench)%timing = bench(ibench)%getTiming(minsec = 0.07_RK)
47 end do
48 write(rootUnit,"(*(g0,:,','))") degree, rootCount
49 write(fileUnit,"(*(g0,:,','))") degree, (bench(ibench)%timing%mean, ibench = 1, size(bench))
50
51 end do loopOverArraySize
52 write(*,"(*(g0,:,' '))") dumsum
53 write(*,"(*(g0,:,' '))")
54
55 close(fileUnit)
56
57contains
58
59 !%%%%%%%%%%%%%%%%%%%%
60 ! procedure wrappers.
61 !%%%%%%%%%%%%%%%%%%%%
62
63 subroutine setOverhead()
64 call getDummy()
65 end subroutine
66
67 subroutine getDummy()
68 dumsum = dumsum + rootCount(ibench)
69 end subroutine
70
71 subroutine setPolyRootEigen()
72 use pm_polynomial, only: setPolyRoot, method => eigen
73 call setPolyRoot(root(1 : degree), rootCount(ibench), coef(0 : degree), method)
74 call getDummy()
75 end subroutine
76
77 subroutine setPolyRootJenkins()
78 use pm_polynomial, only: setPolyRoot, method => jenkins
79 call setPolyRoot(root(1 : degree), rootCount(ibench), coef(0 : degree), method)
80 call getDummy()
81 end subroutine
82
83 subroutine setPolyRootLaguerre()
84 use pm_polynomial, only: setPolyRoot, method => laguerre
85 call setPolyRoot(root(1 : degree), rootCount(ibench), coef(0 : degree), method)
86 call getDummy()
87 end subroutine
88
89 subroutine setPolyRootSGL()
90 use pm_polynomial, only: setPolyRoot, method => sgl
91 call setPolyRoot(root(1 : degree), rootCount(ibench), coef(0 : degree), method)
92 rootCount(ibench) = root(1)%re
93 call getDummy()
94 end subroutine
95
96end program benchmark
Allocate or resize (shrink or expand) an input allocatable scalar string or array of rank 1....
Generate and return an object of type timing_type containing the benchmark timing information and sta...
Definition: pm_bench.F90:574
Return a uniform random scalar or contiguous array of arbitrary rank of randomly uniformly distribute...
Return the roots of a polynomial of arbitrary degree specified by its coefficients coef.
This module contains procedures and generic interfaces for resizing allocatable arrays of various typ...
This module contains abstract interfaces and types that facilitate benchmarking of different procedur...
Definition: pm_bench.F90:41
This module contains classes and procedures for computing various statistical quantities related to t...
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
integer, parameter RK
The default real kind in the ParaMonte library: real64 in Fortran, c_double in C-Fortran Interoperati...
Definition: pm_kind.F90:543
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540
integer, parameter RKD
The double precision real kind in Fortran mode. On most platforms, this is an 64-bit real kind.
Definition: pm_kind.F90:568
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539
integer, parameter RKH
The scalar integer constant of intrinsic default kind, representing the highest-precision real kind t...
Definition: pm_kind.F90:858
This module contains procedures and generic interfaces for performing various mathematical operations...
type(eigen_type), parameter eigen
This is a scalar parameter object of type eigen_type that is exclusively used to signify the use of E...
type(jenkins_type), parameter jenkins
This is a scalar parameter object of type jenkins_type that is exclusively used to signify the use of...
type(laguerre_type), parameter laguerre
This is a scalar parameter object of type laguerre_type that is exclusively used to signify the use o...
type(sgl_type), parameter sgl
This is a scalar parameter object of type sgl_type that is exclusively used to signify the use of Sko...
This is the class for creating benchmark and performance-profiling objects.
Definition: pm_bench.F90:386


Example Unix compile command via Intel ifort compiler

1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe


Example Windows Batch compile command via Intel ifort compiler

1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe


Example Unix / MinGW compile command via GNU gfortran compiler

1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe


Postprocessing of the benchmark output

1#!/usr/bin/env python
2
3import matplotlib.pyplot as plt
4import pandas as pd
5import numpy as np
6
7import os
8dirname = os.path.basename(os.getcwd())
9
10fontsize = 14
11
12df = pd.read_csv("main.out", delimiter = ",")
13colnames = list(df.columns.values)
14
15
18
19ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
20ax = plt.subplot()
21
22for colname in colnames[1:]:
23 plt.plot( df[colnames[0]].values
24 , df[colname].values
25 , linewidth = 2
26 )
27
28plt.xticks(fontsize = fontsize)
29plt.yticks(fontsize = fontsize)
30ax.set_xlabel(colnames[0], fontsize = fontsize)
31ax.set_ylabel("Runtime [ seconds ]", fontsize = fontsize)
32ax.set_title(" vs. ".join(colnames[1:])+"\nLower is better.", fontsize = fontsize)
33ax.set_xscale("log")
34ax.set_yscale("log")
35plt.minorticks_on()
36plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
37ax.tick_params(axis = "y", which = "minor")
38ax.tick_params(axis = "x", which = "minor")
39ax.legend ( colnames[1:]
40 #, loc='center left'
41 #, bbox_to_anchor=(1, 0.5)
42 , fontsize = fontsize
43 )
44
45plt.tight_layout()
46plt.savefig("benchmark." + dirname + ".runtime.png")
47
48
51
52ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
53ax = plt.subplot()
54
55plt.plot( df[colnames[0]].values
56 , np.ones(len(df[colnames[0]].values))
57 , linestyle = "--"
58 #, color = "black"
59 , linewidth = 2
60 )
61for colname in colnames[2:]:
62 plt.plot( df[colnames[0]].values
63 , df[colname].values / df[colnames[1]].values
64 , linewidth = 2
65 )
66
67plt.xticks(fontsize = fontsize)
68plt.yticks(fontsize = fontsize)
69ax.set_xlabel(colnames[0], fontsize = fontsize)
70ax.set_ylabel("Runtime compared to {}".format(colnames[1]), fontsize = fontsize)
71ax.set_title("Runtime Ratio Comparison. Lower means faster.\nLower than 1 means faster than {}.".format(colnames[1]), fontsize = fontsize)
72ax.set_xscale("log")
73ax.set_yscale("log")
74plt.minorticks_on()
75plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
76ax.tick_params(axis = "y", which = "minor")
77ax.tick_params(axis = "x", which = "minor")
78ax.legend ( colnames[1:]
79 #, bbox_to_anchor = (1, 0.5)
80 #, loc = "center left"
81 , fontsize = fontsize
82 )
83
84plt.tight_layout()
85plt.savefig("benchmark." + dirname + ".runtime.ratio.png")
86
87
90
91df = pd.read_csv("root.count", delimiter = ",")
92colnames = list(df.columns.values)
93
94ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
95ax = plt.subplot()
96
97plt.plot( range(1, df[colnames[0]].values[-1] * 2)
98 , range(1, df[colnames[0]].values[-1] * 2)
99 , linestyle = "--"
100 , color = "black"
101 , linewidth = 2
102 )
103for colname in colnames[1:-1]:
104 plt.plot( df[colnames[0]].values
105 , df[colname].values
106 , linewidth = 2
107 )
108
109plt.xticks(fontsize = fontsize)
110plt.yticks(fontsize = fontsize)
111ax.set_xlabel(colnames[0], fontsize = fontsize)
112ax.set_ylabel("Number of Roots Found", fontsize = fontsize)
113ax.set_title(" vs. ".join(colnames[1:])+"\nHigher is better.", fontsize = fontsize)
114ax.set_xscale("log")
115ax.set_yscale("log")
116plt.minorticks_on()
117plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
118ax.tick_params(axis = "y", which = "minor")
119ax.tick_params(axis = "x", which = "minor")
120ax.legend ( ["Equality"] + list(colnames[1:])
121 #, loc='center left'
122 #, bbox_to_anchor=(1, 0.5)
123 , fontsize = fontsize
124 )
125
126plt.tight_layout()
127plt.savefig("benchmark." + dirname + ".root.count.png")


Visualization of the benchmark output


Benchmark moral

  1. Among all root finding algorithms, jenkins_type appears to be the fastest.
  2. The eigen_type method also tends to offer a comparably good performance.
  3. Unlike the above two, laguerre_type algorithm tends to significantly trail behind both in performance and reliability in finding all roots of the polynomial.
Test:
test_pm_polynomial
Todo:
High Priority: A generic interface setPolyCoef(coef, root) must be added that computes polynomial coefficients from its roots using recursive FFT-based polynomial multiplications.
See the commented-out generic interface setPolyCoef within this module as the starting point.


Final Remarks


If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

  1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
  2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

      Copyright 2012 Jan Skowron & Andrew Gould

      Licensed under the Apache License, Version 2.0 (the "License");
      you may not use this file except in compliance with the License.
      You may obtain a copy of the License at

          http://www.apache.org/licenses/LICENSE-2.0

      Unless required by applicable law or agreed to in writing, software
      distributed under the License is distributed on an "AS IS" BASIS,
      WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
      See the License for the specific language governing permissions and
      limitations under the License.

      -------------------------------------------------------------------!

      The authors also make this file available under the terms of
      GNU Lesser General Public License version 2 or any later version.
      (text of the LGPL licence 2 in NOTICE file)

      -------------------------------------------------------------------!

      A custom in the scientific comunity is (regardless of the licence
      you chose to use or distribute this software under)
      that if this code was important in the scientific process or
      for the results of your scientific work, we kindly ask you for the
      appropriate citation of the Paper (Skowron & Gould 2012), and
      we would be greatful if you pass the information about
      the proper citation to anyone whom you redistribute this software to.
Author:
Fatemeh Bagheri, Tuesday 08:49 PM, August 10, 2021, Dallas, TX

Variable Documentation

◆ eigen

type(eigen_type), parameter pm_polynomial::eigen = eigen_type()

This is a scalar parameter object of type eigen_type that is exclusively used to signify the use of Eigenvalue method of root-finding within an interface of a procedure of the ParaMonte library.

See the root-finding section in the documentation of pm_polynomial for more details about this root-finding method.
For example usage, see the documentation of the target procedure requiring this object.

See also
sgl
eigen
jenkins
laguerre
sgl_type
eigen_type
jenkins_type
laguerre_type
method_type


Final Remarks


If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

  1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
  2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Author:
Amir Shahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas Austin

Definition at line 3337 of file pm_polynomial.F90.

◆ jenkins

type(jenkins_type), parameter pm_polynomial::jenkins = jenkins_type()

This is a scalar parameter object of type jenkins_type that is exclusively used to signify the use of Jenkins-Traub method of root-finding within an interface of a procedure of the ParaMonte library.

See the root-finding section in the documentation of pm_polynomial for more details about this root-finding method.
For example usage, see the documentation of the target procedure requiring this object.

See also
sgl
eigen
jenkins
laguerre
sgl_type
eigen_type
jenkins_type
laguerre_type
method_type


Final Remarks


If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

  1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
  2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Author:
Amir Shahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas Austin

Definition at line 3399 of file pm_polynomial.F90.

◆ laguerre

type(laguerre_type), parameter pm_polynomial::laguerre = laguerre_type()

This is a scalar parameter object of type laguerre_type that is exclusively used to signify the use of Laguerre method of root-finding within an interface of a procedure of the ParaMonte library.

See the root-finding section in the documentation of pm_polynomial for more details about this root-finding method.
For example usage, see the documentation of the target procedure requiring this object.

See also
sgl
eigen
jenkins
laguerre
sgl_type
eigen_type
jenkins_type
laguerre_type
method_type


Final Remarks


If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

  1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
  2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Author:
Amir Shahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas Austin

Definition at line 3461 of file pm_polynomial.F90.

◆ MODULE_NAME

character(*, SK), parameter pm_polynomial::MODULE_NAME = "@pm_polynomial"

Definition at line 443 of file pm_polynomial.F90.

◆ sgl

type(sgl_type), parameter pm_polynomial::sgl = sgl_type()

This is a scalar parameter object of type sgl_type that is exclusively used to signify the use of Skowron-Gould method of root-finding within an interface of a procedure of the ParaMonte library.

See the root-finding section in the documentation of pm_polynomial for more details about this root-finding method.
For example usage, see the documentation of the target procedure requiring this object.

See also
sgl
eigen
jenkins
laguerre
sgl_type
eigen_type
jenkins_type
laguerre_type
method_type


Final Remarks


If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

  1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
  2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Author:
Fatemeh Bagheri, August 28, 2024, 6:12 PM, NASA Goddard Space Flight Center, Washington, D.C.

Definition at line 3275 of file pm_polynomial.F90.