|
ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation. |
|
|
ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation. |
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... | |
Variables | |
| character(*, SK), parameter | MODULE_NAME = "@pm_polynomial" |
| 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... | |
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:
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.
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.
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.
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.
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.
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.
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.
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}
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\).
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.
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.
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).
real or complex kind used to represent the coefficients.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.
real are expected to be generally up to four times faster than the corresponding algorithms that accept coefficients of type complex.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:
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.
setPolyCoef(coef, root) must be added that computes polynomial coefficients from its roots using recursive FFT-based polynomial multiplications.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.
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.
| 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.
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.
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.
Definition at line 3200 of file pm_polynomial.F90.
| 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.
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.
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.
Definition at line 3258 of file pm_polynomial.F90.
| 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.
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.
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.
Definition at line 3316 of file pm_polynomial.F90.
| character(*, SK), parameter pm_polynomial::MODULE_NAME = "@pm_polynomial" |
Definition at line 389 of file pm_polynomial.F90.
1.9.3