# Source code for naginterfaces.library.zeros

# -*- coding: utf-8 -*-
r"""
Module Summary
--------------
Interfaces for the NAG Mark 27.3 zeros Chapter.

zeros - Zeros of Polynomials

This module is concerned with computing the zeros of a polynomial with real or complex coefficients.

--------
naginterfaces.library.examples.zeros :
This subpackage contains examples for the zeros module.
See also the :ref:library_zeros_ex subsection.

Functionality Index
-------------------

**All zeros of cubic**

complex coefficients: :meth:cubic_complex

real coefficients: :meth:cubic_real

**All zeros of polynomial**

complex coefficients

fast modified Laguerre's method: :meth:poly_complex_fpml

modified Laguerre's method: :meth:poly_complex

real coefficients

fast modified Laguerre's method: :meth:poly_real_fpml

modified Laguerre's method: :meth:poly_real

complex coefficients: :meth:quadratic_complex

real coefficients: :meth:quadratic_real

**All zeros of quartic**

complex coefficients: :meth:quartic_complex

real coefficients: :meth:quartic_real

For full information please refer to the NAG Library document

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02intro.html
"""

[docs]def poly_complex_fpml(a, itmax=30, polish=1):
r"""
poly_complex_fpml finds all the roots of a complex polynomial equation, using a fourth-order convergent modification of Laguerre's method.

.. _c02aa-py2-py-doc:

For full information please refer to the NAG Library document for c02aa

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02aaf.html

.. _c02aa-py2-py-parameters:

**Parameters**
**a** : complex, array-like, shape :math:\left(n+1\right)
:math:\mathrm{a}[\textit{i}] must contain the coefficient of :math:z^{{n-\textit{i}}}, for :math:\textit{i} = 0,1,\ldots,n.

**itmax** : int, optional
The maximum number of iterations to be performed.

**polish** : int, optional
Specifies the polishing technique used to refine root approximations.

:math:\mathrm{polish} = 0

No polishing.

:math:\mathrm{polish} = 1

Single iteration of Newton's method.

:math:\mathrm{polish} = 2

Iterative refinement using the compensated Horner's method.

See also Selecting a Polishing Process <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02aaf.html#fcomments_polish>__.

**Returns**
**z** : complex, ndarray, shape :math:\left(n\right)
:math:\mathrm{z}[\textit{j}-1] holds the approximation of the root :math:z_{\textit{j}}, for :math:\textit{j} = 1,2,\ldots,n.

**berr** : float, ndarray, shape :math:\left(n\right)
:math:\mathrm{berr}[\textit{j}-1] holds the relative backward error, :math:\eta \left(z_{\textit{j}}\right), for :math:\textit{j} = 1,2,\ldots,n.

**cond** : float, ndarray, shape :math:\left(n\right)
:math:\mathrm{cond}[\textit{j}-1] holds the condition number, :math:\kappa \left(z_{\textit{j}}\right), for :math:\textit{j} = 1,2,\ldots,n.

**conv** : int, ndarray, shape :math:\left(n\right)
:math:\mathrm{conv}[\textit{j}-1] indicates the convergence status of the root approximation, :math:z_{\textit{j}}, for :math:\textit{j} = 1,2,\ldots,n.

:math:\mathrm{conv}[j-1]\geq 0

Successfully converged after :math:\mathrm{conv}[j-1] iterations.

:math:\mathrm{conv}[j-1] = -1

Failed to converge after :math:\mathrm{itmax} iterations.

:math:\mathrm{conv}[j-1] = -2

Overflow was encountered.

Note: if :math:\mathrm{polish} = 2, :math:\mathrm{conv} refers to convergence in the compensated polishing process.

.. _c02aa-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, the complex variable :math:\mathrm{a}[0] = \left(0.0, 0.0\right).

(errno :math:2)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n \geq 1.

(errno :math:3)
On entry, :math:\mathrm{itmax} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{itmax}\geq 1.

(errno :math:4)
On entry, :math:\mathrm{polish} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{polish} = 0, :math:1 or :math:2.

(errno :math:5)
Convergence has failed for at least one root approximation.

(errno :math:6)
poly_complex_fpml encountered overflow during at least one root approximation.

.. _c02aa-py2-py-notes:

**Notes**
poly_complex_fpml attempts to find all the roots of the :math:n\ th degree complex polynomial equation

.. math::
P\left(z\right) = a_0z^n+a_1z^{{n-1}}+a_2z^{{n-2}} + \cdots +a_{{n-1}}z+a_n = 0\text{.}

The roots are located using a modified form of Laguerre's method, as implemented by Cameron (2018).
An implicit deflation strategy is employed, which allows for high accuracy even when solving high degree polynomials.
Linear (:math:n = 1) and quadratic (:math:n = 2) problems are obtained via the 'standard' closed formulae.

First, initial estimates of the roots are made using a method proposed by Bini (1996), which selects complex numbers along circles of suitable radii.
Updates to each root approximation :math:z_{\textit{j}}, for :math:\textit{j} = 1,2,\ldots,n, are then made using the iterative formula from Petković et al. (1997):

.. math::
\hat{z}_j = z_j-\frac{{n}}{{G_j\pm \sqrt{\left(n-1\right)\left(nH_j-G_j^2\right)}}}\text{,}

where

.. math::
G_j = \frac{{P^{\prime }\left(z_j\right)}}{{P\left(z_j\right)}}-\sum_{1\text{; }\textit{i}\neq j}^{n}{\frac{{1}}{{\left(z_j-z_{\textit{i}}\right)}}}\hspace{1em}\text{and}\hspace{1em}H_j = -\left(\frac{{P^{\prime }\left(z_j\right)}}{{P\left(z_j\right)}}\right)^{\prime }-\sum_{1\text{; }\textit{i}\neq j}^{n}{\frac{{1}}{{\left(z_j-z_{\textit{i}}\right)^2}}}\text{.}

The nearest :math:\hat{z}_j to the current approximation is used, by selecting the sign that maximizes the denominator of the correction term.
The subtraction of the sum terms when computing :math:G_j and :math:H_j determines the implicit deflation strategy of the modified Laguerre method.

The relative backward error of a root approximation :math:z_j is given by

.. math::
\eta \left(z_j\right) = \frac{{\left\lvert P\left(z_j\right)\right\rvert }}{{\alpha \left(z_j\right)}}\text{,}

where :math:\alpha is the perturbed polynomial

.. math::
\alpha \left(z\right) = \sum_{0}^{n}{\left(3.8\textit{k}+1\right)\left\lvert a_{{n-\textit{k}}}\right\rvert \left\lvert z\right\rvert^{\textit{k}}}\text{.}

A root approximation is deemed to have converged if :math:\eta \left(z_j\right)\leq 2\text{ }\times machine precision, at which point updates of that root cease.
If the stopping criterion holds, then the computed root is the exact root of a polynomial whose coefficients are no more perturbed than the floating-point computation of :math:P\left(z_j\right).

The condition number of each root is also computed, as a measure of sensitivity to changes in the coefficients of the polynomial.
It is given by

.. math::
\kappa \left(z_j\right) = \frac{{\alpha \left(z_j\right)}}{{\left\lvert z_j\right\rvert \left\lvert P^{\prime }\left(z_j\right)\right\rvert }}\text{.}

Root approximations can be further refined with optional 'polishing' processes.
A simple polishing process is provided that carries out a single iteration of Newton's method, which proves quick yet often effective.
Alternatively, a compensated polishing process from Cameron and O'Neill (2019) can be applied.
This iterative method combines the implicit deflation of the modified Laguerre method, with the accuracy of evaluating polynomials and their derivatives using the compensated Horner's method from Graillat et al. (2005).
Compensated polishing yields approximations with a limiting accuracy as if computed in twice the working precision.

It is recommended that you read Selecting a Polishing Process <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02aaf.html#fcomments_polish>__ for advice on selecting an appropriate polishing process.

.. _c02aa-py2-py-references:

**References**
Bini, D A, 1996, Numerical computation of polynomial zeros by means of Aberth's method, Numerical Algorithms (13), 179--200, Springer US

Cameron, T R, 2018, An effective implementation of a modified Laguerre method for the roots of a polynomial, Numerical Algorithms, Springer US, https://doi.org/10.1007/s11075-018-0641-9

Cameron, T R and O'Neill, A, 2019, On a compensated polishing technique for polynomial root solvers, To be published,

Graillat, S, Louvet, N, and Langlois, P, 2005, Compensated Horner scheme, Technical Report, Université de Perpignan Via Domitia

Petković, M, Ilić, S, and Tričković, S, 1997, A family of simultaneous zero-finding methods, Computers & Mathematics with Applications (Volume 34) (10), 49--59, https://doi.org/10.1016/S0898-1221(97)00206-X

Wilkinson, J H, 1959, The evaluation of the zeros of ill-conditioned polynomials. Part I, Numerische Mathematik (Volume 1) (1), 150--166, Springer-Verlag

--------
:meth:naginterfaces.library.examples.zeros.poly_complex_fpml_ex.main
"""
raise NotImplementedError

[docs]def poly_real_fpml(a, itmax=30, polish=1):
r"""
poly_real_fpml finds all the roots of a real polynomial equation, using a fourth-order convergent modification of Laguerre's method.

.. _c02ab-py2-py-doc:

For full information please refer to the NAG Library document for c02ab

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02abf.html

.. _c02ab-py2-py-parameters:

**Parameters**
**a** : float, array-like, shape :math:\left(n+1\right)
:math:\mathrm{a}[\textit{i}] must contain the coefficient of :math:z^{{n-\textit{i}}}, for :math:\textit{i} = 0,1,\ldots,n.

**itmax** : int, optional
The maximum number of iterations to be performed.

**polish** : int, optional
Specifies the polishing technique used to refine root approximations.

:math:\mathrm{polish} = 0

No polishing.

:math:\mathrm{polish} = 1

Single iteration of Newton's method.

:math:\mathrm{polish} = 2

Iterative refinement using the compensated Horner's method.

See also Selecting a Polishing Process <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02abf.html#fcomments_polish>__.

**Returns**
**z** : complex, ndarray, shape :math:\left(n\right)
:math:\mathrm{z}[\textit{j}-1] holds the approximation of the root :math:z_{\textit{j}}, for :math:\textit{j} = 1,2,\ldots,n.

**berr** : float, ndarray, shape :math:\left(n\right)
:math:\mathrm{berr}[\textit{j}-1] holds the relative backward error, :math:\eta \left(z_{\textit{j}}\right), for :math:\textit{j} = 1,2,\ldots,n.

**cond** : float, ndarray, shape :math:\left(n\right)
:math:\mathrm{cond}[\textit{j}-1] holds the condition number, :math:\kappa \left(z_{\textit{j}}\right), for :math:\textit{j} = 1,2,\ldots,n.

**conv** : int, ndarray, shape :math:\left(n\right)
:math:\mathrm{conv}[\textit{j}-1] indicates the convergence status of the root approximation, :math:z_{\textit{j}}, for :math:\textit{j} = 1,2,\ldots,n.

:math:\mathrm{conv}[j-1]\geq 0

Successfully converged after :math:\mathrm{conv}[j-1] iterations.

:math:\mathrm{conv}[j-1] = -1

Failed to converge after :math:\mathrm{itmax} iterations.

:math:\mathrm{conv}[j-1] = -2

Overflow was encountered.

Note: if :math:\mathrm{polish} = 2, :math:\mathrm{conv} refers to convergence in the compensated polishing process.

.. _c02ab-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, the real variable :math:\mathrm{a}[0] = 0.0.

(errno :math:2)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n \geq 1.

(errno :math:3)
On entry, :math:\mathrm{itmax} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{itmax}\geq 1.

(errno :math:4)
On entry, :math:\mathrm{polish} = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:\mathrm{polish} = 0, :math:1 or :math:2.

(errno :math:5)
Convergence has failed for at least one root approximation.

(errno :math:6)
poly_real_fpml encountered overflow during at least one root approximation.

.. _c02ab-py2-py-notes:

**Notes**
poly_real_fpml attempts to find all the roots of the :math:n\ th degree real polynomial equation

.. math::
P\left(z\right) = a_0z^n+a_1z^{{n-1}}+a_2z^{{n-2}} + \cdots +a_{{n-1}}z+a_n = 0\text{.}

The roots are located using a modified form of Laguerre's method, as implemented by Cameron (2018).

poly_real_fpml is a wrapper around the corresponding complex routine :meth:poly_complex_fpml.

The relative backward error of a root approximation :math:z_j is given by

.. math::
\eta \left(z_j\right) = \frac{{\left\lvert P\left(z_j\right)\right\rvert }}{{\alpha \left(z_j\right)}}\text{,}

where :math:\alpha is the perturbed polynomial

.. math::
\alpha \left(z\right) = \sum_{0}^{n}{\left(3.8\textit{k}+1\right)\left\lvert a_{{n-\textit{k}}}\right\rvert \left\lvert z\right\rvert^{\textit{k}}}\text{.}

A root approximation is deemed to have converged if :math:\eta \left(z_j\right)\leq 2\text{ }\times machine precision, at which point updates of that root cease.
If the stopping criterion holds, then the computed root is the exact root of a polynomial whose coefficients are no more perturbed than the floating-point computation of :math:P\left(z_j\right).

The condition number of each root is also computed, as a measure of sensitivity to changes in the coefficients of the polynomial.
It is given by

.. math::
\kappa \left(z_j\right) = \frac{{\alpha \left(z_j\right)}}{{\left\lvert z_j\right\rvert \left\lvert P^{\prime }\left(z_j\right)\right\rvert }}\text{.}

Root approximations can be further refined with optional 'polishing' processes.
A simple polishing process is provided that carries out a single iteration of Newton's method, which proves quick yet often effective.
Alternatively, a compensated polishing process from Cameron and O'Neill (2019) can be applied.
This iterative method combines the implicit deflation of the modified Laguerre method, with the accuracy of evaluating polynomials and their derivatives using the compensated Horner's method from Graillat et al. (2005).
Compensated polishing yields approximations with a limiting accuracy as if computed in twice the working precision.

It is recommended that you read Selecting a Polishing Process <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02abf.html#fcomments_polish>__ for advice on selecting an appropriate polishing process.

.. _c02ab-py2-py-references:

**References**
Bini, D A, 1996, Numerical computation of polynomial zeros by means of Aberth's method, Numerical Algorithms (13), 179--200, Springer US

Cameron, T R, 2018, An effective implementation of a modified Laguerre method for the roots of a polynomial, Numerical Algorithms, Springer US, https://doi.org/10.1007/s11075-018-0641-9

Cameron, T R and O'Neill, A, 2019, On a compensated polishing technique for polynomial root solvers, To be published,

Graillat, S, Louvet, N, and Langlois, P, 2005, Compensated Horner scheme, Technical Report, Université de Perpignan Via Domitia

Petković, M, Ilić, S, and Tričković, S, 1997, A family of simultaneous zero-finding methods, Computers & Mathematics with Applications (Volume 34) (10), 49--59, https://doi.org/10.1016/S0898-1221(97)00206-X

Wilkinson, J H, 1959, The evaluation of the zeros of ill-conditioned polynomials. Part I, Numerische Mathematik (Volume 1) (1), 150--166, Springer-Verlag
"""
raise NotImplementedError

[docs]def poly_complex(a, scal=True):
r"""
poly_complex finds all the roots of a complex polynomial equation, using a variant of Laguerre's method.

.. deprecated:: 27.1.0.0
poly_complex is deprecated.
Please use :meth:poly_complex_fpml instead.
See also the :ref:Replacement Calls <replace> document.

.. _c02af-py2-py-doc:

For full information please refer to the NAG Library document for c02af

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02aff.html

.. _c02af-py2-py-parameters:

**Parameters**
**a** : complex, array-like, shape :math:\left(n+1\right)
:math:\mathrm{Re}\left(\mathrm{a}[\textit{i}]\right) and :math:\mathrm{Im}\left(\mathrm{a}[\textit{i}]\right) must contain the real and imaginary parts of :math:a_{\textit{i}} (i.e., the coefficient of :math:z^{{n-\textit{i}}}), for :math:\textit{i} = 0,1,\ldots,n.

**scal** : bool, optional
Indicates whether or not the polynomial is to be scaled. See Further Comments <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02aff.html#fcomments>__ for advice on when it may be preferable to set :math:\mathrm{scal} = \mathbf{False} and for a description of the scaling strategy.

**Returns**
**z** : complex, ndarray, shape :math:\left(n\right)
The roots.

.. _c02af-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{a}[0] = 0.0.

Constraint: :math:\mathrm{a}[0]\neq 0.0.

(errno :math:1)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n \geq 1.

(errno :math:2)
The iterative procedure has failed to converge. This error is very unlikely to occur. If it does, please contact NAG <https://www.nag.com>__ immediately, as some basic assumption for the arithmetic has been violated.

(errno :math:3)
poly_complex cannot evaluate :math:P\left(z\right) near some of its zeros without underflow. If this message occurs please contact NAG <https://www.nag.com>__.

(errno :math:3)
poly_complex cannot evaluate :math:P\left(z\right) near some of its zeros without overflow. If this message occurs please contact NAG <https://www.nag.com>__.

.. _c02af-py2-py-notes:

**Notes**
In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.

poly_complex attempts to find all the roots of the :math:n\ th degree complex polynomial equation

.. math::
P\left(z\right) = a_0z^n+a_1z^{{n-1}}+a_2z^{{n-2}} + \cdots +a_{{n-1}}z+a_n = 0\text{.}

The roots are located using a modified form of Laguerre's method, originally proposed by Smith (1967).

The method of Laguerre (see Wilkinson (1965)) can be described by the iterative scheme

.. math::
L\left(z_k\right) = z_{{k+1}}-z_k = \frac{{{-n}P\left(z_k\right)}}{{P^{\prime }\left(z_k\right)\pm \sqrt{H\left(z_k\right)}}}\text{,}

where :math:H\left(z_k\right) = \left(n-1\right)\left[\left(n-1\right){\left(P^{\prime }\left(z_k\right)\right)}^2-nP\left(z_k\right)P^{{\prime \prime }}\left(z_k\right)\right] and :math:z_0 is specified.

The sign in the denominator is chosen so that the modulus of the Laguerre step at :math:z_k, viz. :math:\left\lvert L\left(z_k\right)\right\rvert, is as small as possible.
The method can be shown to be cubically convergent for isolated roots (real or complex) and linearly convergent for multiple roots.

The function generates a sequence of iterates :math:z_1,z_2,z_3,\ldots \text{}, such that :math:\left\lvert P\left(z_{{k+1}}\right)\right\rvert < \left\lvert P\left(z_k\right)\right\rvert and ensures that :math:z_{{k+1}}+L\left(z_{{k+1}}\right) 'roughly' lies inside a circular region of radius :math:\left\lvert F\right\rvert about :math:z_k known to contain a zero of :math:P\left(z\right); that is, :math:\left\lvert L\left(z_{{k+1}}\right)\right\rvert \leq \left\lvert F\right\rvert, where :math:F denotes the Fejér bound (see Marden (1966)) at the point :math:z_k.
Following Smith (1967), :math:F is taken to be :math:\mathrm{min}\left(B, {1.445nR}\right), where :math:B is an upper bound for the magnitude of the smallest zero given by

.. math::
B = 1.0001\times \mathrm{min}\left({\sqrt{n}L\left(z_k\right)}, \left\lvert r_1\right\rvert, \left\lvert a_n/a_0\right\rvert^{{1/n}}\right)\text{,}

:math:r_1 is the zero :math:X of smaller magnitude of the quadratic equation

.. math::
\frac{{P^{{\prime \prime }}\left(z_k\right)}}{{2n\left(n-1\right)}}X^2+\frac{{P^{\prime }\left(z_k\right)}}{n}X+\frac{1}{2}P\left(z_k\right) = 0

and the Cauchy lower bound :math:R for the smallest zero is computed (using Newton's Method) as the positive root of the polynomial equation

.. math::
\left\lvert a_0\right\rvert z^n+\left\lvert a_1\right\rvert z^{{n-1}}+\left\lvert a_2\right\rvert z^{{n-2}} + \cdots +\left\lvert a_{{n-1}}\right\rvert z-\left\lvert a_n\right\rvert = 0\text{.}

Starting from the origin, successive iterates are generated according to the rule :math:z_{{k+1}} = z_k+L\left(z_k\right), for :math:k = 1,2,3,\ldots \text{}, and :math:L\left(z_k\right) is 'adjusted' so that :math:\left\lvert P\left(z_{{k+1}}\right)\right\rvert < \left\lvert P\left(z_k\right)\right\rvert and :math:\left\lvert L\left(z_{{k+1}}\right)\right\rvert \leq \left\lvert F\right\rvert.
The iterative procedure terminates if :math:P\left(z_{{k+1}}\right) is smaller in absolute value than the bound on the rounding error in :math:P\left(z_{{k+1}}\right) and the current iterate :math:z_p = z_{{k+1}} is taken to be a zero of :math:P\left(z\right).
The deflated polynomial :math:\tilde{P}\left(z\right) = P\left(z\right)/\left(z-z_p\right) of degree :math:n-1 is then formed, and the above procedure is repeated on the deflated polynomial until :math:n < 3, whereupon the remaining roots are obtained via the 'standard' closed formulae for a linear (:math:n = 1) or quadratic (:math:n = 2) equation.

Note that :meth:quadratic_complex, :meth:cubic_complex and :meth:quartic_complex can be used to obtain the roots of a quadratic, cubic (:math:n = 3) and quartic (:math:n = 4) polynomial, respectively.

.. _c02af-py2-py-references:

**References**
Marden, M, 1966, Geometry of polynomials, Mathematical Surveys (3), American Mathematical Society, Providence, RI

Smith, B T, 1967, ZERPOL: a zero finding algorithm for polynomials using Laguerre's method, Technical Report, Department of Computer Science, University of Toronto, Canada

Thompson, K W, 1991, Error analysis for polynomial solvers, Fortran Journal (Volume 3) (3), 10--13

Wilkinson, J H, 1965, The Algebraic Eigenvalue Problem, Oxford University Press, Oxford
"""
raise NotImplementedError

[docs]def poly_real(a, scal=True):
r"""
poly_real finds all the roots of a real polynomial equation, using a variant of Laguerre's method.

.. deprecated:: 27.1.0.0
poly_real is deprecated.
Please use :meth:poly_real_fpml instead.
See also the :ref:Replacement Calls <replace> document.

.. _c02ag-py2-py-doc:

For full information please refer to the NAG Library document for c02ag

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02agf.html

.. _c02ag-py2-py-parameters:

**Parameters**
**a** : float, array-like, shape :math:\left(n+1\right)
:math:\mathrm{a}[\textit{i}] must contain :math:a_{\textit{i}} (i.e., the coefficient of :math:z^{{n-\textit{i}}}), for :math:\textit{i} = 0,1,\ldots,n.

**scal** : bool, optional
Indicates whether or not the polynomial is to be scaled. See Further Comments <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02agf.html#fcomments>__ for advice on when it may be preferable to set :math:\mathrm{scal} = \mathbf{False} and for a description of the scaling strategy.

**Returns**
**z** : complex, ndarray, shape :math:\left(n\right)
The roots. Complex conjugate pairs of roots are stored in consecutive elements of :math:\mathrm{z}.

.. _c02ag-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{a}[0] = 0.0.

Constraint: :math:\mathrm{a}[0]\neq 0.0.

(errno :math:1)
On entry, :math:n = \langle\mathit{\boldsymbol{value}}\rangle.

Constraint: :math:n \geq 1.

(errno :math:2)
The iterative procedure has failed to converge. This error is very unlikely to occur. If it does, please contact NAG <https://www.nag.com>__ immediately, as some basic assumption for the arithmetic has been violated.

(errno :math:3)
poly_real cannot evaluate :math:P\left(z\right) near some of its zeros without underflow. If this message occurs please contact NAG <https://www.nag.com>__.

(errno :math:3)
poly_real cannot evaluate :math:P\left(z\right) near some of its zeros without overflow. If this message occurs please contact NAG <https://www.nag.com>__.

.. _c02ag-py2-py-notes:

**Notes**
In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.

poly_real attempts to find all the roots of the :math:n\ th degree real polynomial equation

.. math::
P\left(z\right) = a_0z^n+a_1z^{{n-1}}+a_2z^{{n-2}} + \cdots +a_{{n-1}}z+a_n = 0\text{.}

The roots are located using a modified form of Laguerre's method, originally proposed by Smith (1967).

The method of Laguerre (see Wilkinson (1965)) can be described by the iterative scheme

.. math::
L\left(z_k\right) = z_{{k+1}}-z_k = \frac{{{-n}P\left(z_k\right)}}{{P^{\prime }\left(z_k\right)\pm \sqrt{H\left(z_k\right)}}}\text{,}

where :math:H\left(z_k\right) = \left(n-1\right)\left[\left(n-1\right){\left(P^{\prime }\left(z_k\right)\right)}^2-nP\left(z_k\right)P^{{\prime \prime }}\left(z_k\right)\right] and :math:z_0 is specified.

The sign in the denominator is chosen so that the modulus of the Laguerre step at :math:z_k, viz. :math:\left\lvert L\left(z_k\right)\right\rvert, is as small as possible.
The method can be shown to be cubically convergent for isolated roots (real or complex) and linearly convergent for multiple roots.

The function generates a sequence of iterates :math:z_1,z_2,z_3,\ldots \text{}, such that :math:\left\lvert P\left(z_{{k+1}}\right)\right\rvert < \left\lvert P\left(z_k\right)\right\rvert and ensures that :math:z_{{k+1}}+L\left(z_{{k+1}}\right) 'roughly' lies inside a circular region of radius :math:\left\lvert F\right\rvert about :math:z_k known to contain a zero of :math:P\left(z\right); that is, :math:\left\lvert L\left(z_{{k+1}}\right)\right\rvert \leq \left\lvert F\right\rvert, where :math:F denotes the Fejér bound (see Marden (1966)) at the point :math:z_k.
Following Smith (1967), :math:F is taken to be :math:\mathrm{min}\left(B, {1.445nR}\right), where :math:B is an upper bound for the magnitude of the smallest zero given by

.. math::
B = 1.0001\times \mathrm{min}\left({\sqrt{n}L\left(z_k\right)}, \left\lvert r_1\right\rvert, \left\lvert a_n/a_0\right\rvert^{{1/n}}\right)\text{,}

:math:r_1 is the zero :math:X of smaller magnitude of the quadratic equation

.. math::
\frac{{P^{{\prime \prime }}\left(z_k\right)}}{{2n\left(n-1\right)}}X^2+\frac{{P^{\prime }\left(z_k\right)}}{n}X+\frac{1}{2}P\left(z_k\right) = 0

and the Cauchy lower bound :math:R for the smallest zero is computed (using Newton's Method) as the positive root of the polynomial equation

.. math::
\left\lvert a_0\right\rvert z^n+\left\lvert a_1\right\rvert z^{{n-1}}+\left\lvert a_2\right\rvert z^{{n-2}} + \cdots +\left\lvert a_{{n-1}}\right\rvert z-\left\lvert a_n\right\rvert = 0\text{.}

Starting from the origin, successive iterates are generated according to the rule :math:z_{{k+1}} = z_k+L\left(z_k\right), for :math:k = 1,2,3,\ldots \text{}, and :math:L\left(z_k\right) is 'adjusted' so that :math:\left\lvert P\left(z_{{k+1}}\right)\right\rvert < \left\lvert P\left(z_k\right)\right\rvert and :math:\left\lvert L\left(z_{{k+1}}\right)\right\rvert \leq \left\lvert F\right\rvert.
The iterative procedure terminates if :math:P\left(z_{{k+1}}\right) is smaller in absolute value than the bound on the rounding error in :math:P\left(z_{{k+1}}\right) and the current iterate :math:z_p = z_{{k+1}} is taken to be a zero of :math:P\left(z\right) (as is its conjugate :math:\bar{z}_p if :math:z_p is complex).
The deflated polynomial :math:\tilde{P}\left(z\right) = P\left(z\right)/\left(z-z_p\right) of degree :math:n-1 if :math:z_p is real (:math:\tilde{P}\left(z\right) = P\left(z\right)/\left(\left(z-z_p\right)\left(z-\bar{z}_p\right)\right) of degree :math:n-2 if :math:z_p is complex) is then formed, and the above procedure is repeated on the deflated polynomial until :math:n < 3, whereupon the remaining roots are obtained via the 'standard' closed formulae for a linear (:math:n = 1) or quadratic (:math:n = 2) equation.

Note that :meth:quadratic_real, :meth:cubic_real and :meth:quartic_real can be used to obtain the roots of a quadratic, cubic (:math:n = 3) and quartic (:math:n = 4) polynomial, respectively.

.. _c02ag-py2-py-references:

**References**
Marden, M, 1966, Geometry of polynomials, Mathematical Surveys (3), American Mathematical Society, Providence, RI

Smith, B T, 1967, ZERPOL: a zero finding algorithm for polynomials using Laguerre's method, Technical Report, Department of Computer Science, University of Toronto, Canada

Thompson, K W, 1991, Error analysis for polynomial solvers, Fortran Journal (Volume 3) (3), 10--13

Wilkinson, J H, 1965, The Algebraic Eigenvalue Problem, Oxford University Press, Oxford
"""
raise NotImplementedError

r"""
quadratic_complex determines the roots of a quadratic equation with complex coefficients.

.. _c02ah-py2-py-doc:

For full information please refer to the NAG Library document for c02ah

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02ahf.html

.. _c02ah-py2-py-parameters:

**Parameters**
**a** : complex
:math:a, the coefficient of :math:z^2.

**b** : complex
:math:b, the coefficient of :math:z.

**c** : complex
:math:c, the constant coefficient.

**Returns**
**z_small** : complex
The smallest root in magnitude.

**z_large** : complex
The largest root in magnitude.

.. _c02ah-py2-py-errors:

**Warns**
**NagAlgorithmicWarning**
(errno :math:1)
On entry, :math:\mathrm{a} = 0.0.

(errno :math:2)
On entry, :math:\mathrm{a} = 0.0 and :math:\mathrm{b} = 0.0.

(errno :math:3)
On entry, :math:\mathrm{a} = 0.0 and the root :math:{-\mathrm{c}}/\mathrm{b} overflows: :math:\mathrm{Re}\left(\mathrm{a}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Re}\left(\mathrm{c}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Re}\left(\mathrm{b}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Im}\left(\mathrm{a}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Im}\left(\mathrm{c}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Im}\left(\mathrm{b}\right) = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:4)
On entry, :math:\mathrm{c} = 0.0 and the root :math:{-\mathrm{b}}/\mathrm{a} overflows: :math:\mathrm{Re}\left(\mathrm{c}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Re}\left(\mathrm{b}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Re}\left(\mathrm{a}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Im}\left(\mathrm{c}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Im}\left(\mathrm{b}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Im}\left(\mathrm{a}\right) = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:5)
On entry, :math:B = \mathrm{max}\left(\left\lvert \mathrm{Re}\left(\mathrm{b}\right)\right\rvert, \left\lvert \mathrm{Im}\left(\mathrm{b}\right)\right\rvert \right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:A = \mathrm{max}\left(\left\lvert \mathrm{Re}\left(\mathrm{a}\right)\right\rvert, \left\lvert \mathrm{Im}\left(\mathrm{a}\right)\right\rvert \right) = \langle\mathit{\boldsymbol{value}}\rangle and :math:C = \mathrm{max}\left(\left\lvert \mathrm{Re}\left(\mathrm{c}\right)\right\rvert, \left\lvert \mathrm{Im}\left(\mathrm{c}\right)\right\rvert \right) = \langle\mathit{\boldsymbol{value}}\rangle. :math:B is so large that :math:B^2 is indistinguishable from :math:\left(B^2-4\times A\times C\right) and the root :math:{-\mathrm{b}}/\mathrm{a} overflows. :math:\mathrm{Re}\left(\mathrm{b}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Im}\left(\mathrm{b}\right) = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{Re}\left(\mathrm{a}\right) = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{Im}\left(\mathrm{a}\right) = \langle\mathit{\boldsymbol{value}}\rangle.

.. _c02ah-py2-py-notes:

**Notes**
No equivalent traditional C interface for this routine exists in the NAG Library.

quadratic_complex attempts to find the roots of the quadratic equation :math:az^2+bz+c = 0 (where :math:a, :math:b and :math:c are complex coefficients), by carefully evaluating the 'standard' closed formula

.. math::
z = \frac{{-b\pm \sqrt{b^2-4ac}}}{{2a}}\text{.}

It is based on the function CQDRTC from Smith (1967).

**Note:** it is not necessary to scale the coefficients prior to calling the function.

.. _c02ah-py2-py-references:

**References**
Smith, B T, 1967, ZERPOL: a zero finding algorithm for polynomials using Laguerre's method, Technical Report, Department of Computer Science, University of Toronto, Canada
"""
raise NotImplementedError

r"""
quadratic_real determines the roots of a quadratic equation with real coefficients.

.. _c02aj-py2-py-doc:

For full information please refer to the NAG Library document for c02aj

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02ajf.html

.. _c02aj-py2-py-parameters:

**Parameters**
**a** : float
Must contain :math:a, the coefficient of :math:z^2.

**b** : float
Must contain :math:b, the coefficient of :math:z.

**c** : float
Must contain :math:c, the constant coefficient.

**Returns**
**z_small** : complex
The smallest root in magnitude.

**z_large** : complex
The largest root in magnitude.

.. _c02aj-py2-py-errors:

**Warns**
**NagAlgorithmicWarning**
(errno :math:1)
On entry, :math:\mathrm{a} = 0.0.

(errno :math:2)
On entry, :math:\mathrm{a} = 0.0 and :math:\mathrm{b} = 0.0.

(errno :math:3)
On entry, :math:\mathrm{a} = 0.0 and the root :math:{-\mathrm{c}}/\mathrm{b} overflows: :math:\mathrm{a} = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{c} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{b} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:4)
On entry, :math:\mathrm{c} = 0.0 and the root :math:{-\mathrm{b}}/\mathrm{a} overflows: :math:\mathrm{c} = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{b} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{a} = \langle\mathit{\boldsymbol{value}}\rangle.

(errno :math:5)
On entry, :math:\mathrm{b} is so large that :math:\mathrm{b}^2 is indistinguishable from :math:\left(\mathrm{b}^2-4\times \mathrm{a}\times \mathrm{c}\right) and the root :math:{-\mathrm{b}}/\mathrm{a} overflows: :math:\mathrm{b} = \langle\mathit{\boldsymbol{value}}\rangle, :math:\mathrm{a} = \langle\mathit{\boldsymbol{value}}\rangle and :math:\mathrm{c} = \langle\mathit{\boldsymbol{value}}\rangle.

.. _c02aj-py2-py-notes:

**Notes**
No equivalent traditional C interface for this routine exists in the NAG Library.

quadratic_real attempts to find the roots of the quadratic equation :math:az^2+bz+c = 0 (where :math:a, :math:b and :math:c are real coefficients), by carefully evaluating the 'standard' closed formula

.. math::
z = \frac{{-b\pm \sqrt{b^2-4ac}}}{{2a}}\text{.}

It is based on the function QDRTC from Smith (1967).

**Note:** it is not necessary to scale the coefficients prior to calling the function.

.. _c02aj-py2-py-references:

**References**
Smith, B T, 1967, ZERPOL: a zero finding algorithm for polynomials using Laguerre's method, Technical Report, Department of Computer Science, University of Toronto, Canada
"""
raise NotImplementedError

[docs]def cubic_real(u, r, s, t):
r"""
cubic_real determines the roots of a cubic equation with real coefficients.

.. _c02ak-py2-py-doc:

For full information please refer to the NAG Library document for c02ak

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02akf.html

.. _c02ak-py2-py-parameters:

**Parameters**
**u** : float
:math:u, the coefficient of :math:z^3.

**r** : float
:math:r, the coefficient of :math:z^2.

**s** : float
:math:s, the coefficient of :math:z.

**t** : float
:math:t, the constant coefficient.

**Returns**
**zero** : complex, ndarray, shape :math:\left(3\right)
:math:\mathrm{zero}[i-1] contains the :math:i\ th root.

**errest** : float, ndarray, shape :math:\left(3\right)
:math:\mathrm{errest}[i-1] contains an approximate error estimate for the :math:i\ th root.

.. _c02ak-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{u} = 0.0.

Constraint: :math:\mathrm{u}\neq 0.0.

(errno :math:2)
The companion matrix :math:H cannot be formed without overflow.

(errno :math:3)
Failure to converge in :meth:lapackeig.dhseqr <naginterfaces.library.lapackeig.dhseqr>.

.. _c02ak-py2-py-notes:

**Notes**
In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.

cubic_real attempts to find the roots of the cubic equation

.. math::
uz^3+rz^2+sz+t = 0\text{,}

where :math:u, :math:r, :math:s and :math:t are real coefficients with :math:u\neq 0.
The roots are located by finding the eigenvalues of the associated :math:3\times 3 (upper Hessenberg) companion matrix :math:H given by

.. math::
H = \begin{pmatrix}0&0&-t/u\\1&0&-s/u\\0&1&-r/u\end{pmatrix}\text{.}

The eigenvalues are obtained by a call to :meth:lapackeig.dhseqr <naginterfaces.library.lapackeig.dhseqr>.
Further details can be found in Further Comments <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02akf.html#fcomments>__.

To obtain the roots of a quartic equation, :meth:quartic_real can be used.

.. _c02ak-py2-py-references:

**References**
Golub, G H and Van Loan, C F, 1996, Matrix Computations, (3rd Edition), Johns Hopkins University Press, Baltimore
"""
raise NotImplementedError

[docs]def quartic_real(e, a, b, c, d):
r"""
quartic_real determines the roots of a quartic equation with real coefficients.

.. _c02al-py2-py-doc:

For full information please refer to the NAG Library document for c02al

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02alf.html

.. _c02al-py2-py-parameters:

**Parameters**
**e** : float
:math:e, the coefficient of :math:z^4.

**a** : float
:math:a, the coefficient of :math:z^3.

**b** : float
:math:b, the coefficient of :math:z^2.

**c** : float
:math:c, the coefficient of :math:z.

**d** : float
:math:d, the constant coefficient.

**Returns**
**zero** : complex, ndarray, shape :math:\left(4\right)
:math:\mathrm{zero}[i-1] contains the :math:i\ th root.

**errest** : float, ndarray, shape :math:\left(4\right)
:math:\mathrm{errest}[i-1] contains an approximate error estimate for the :math:i\ th root.

.. _c02al-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{e} = 0.0.

Constraint: :math:\mathrm{e}\neq 0.0.

(errno :math:2)
The companion matrix :math:H cannot be formed without overflow.

(errno :math:3)
Failure to converge in :meth:lapackeig.dhseqr <naginterfaces.library.lapackeig.dhseqr>.

.. _c02al-py2-py-notes:

**Notes**
In the NAG Library the traditional C interface for this routine uses a different algorithmic base. Please contact NAG if you have any questions about compatibility.

quartic_real attempts to find the roots of the quartic equation

.. math::
ez^4+az^3+bz^2+cz+d = 0\text{,}

where :math:e, :math:a, :math:b, :math:c and :math:d are real coefficients with :math:e\neq 0.
The roots are located by finding the eigenvalues of the associated :math:4\times 4 (upper Hessenberg) companion matrix :math:H given by

.. math::
H = \begin{pmatrix}0&0&0&-d/e\\1&0&0&-c/e\\0&1&0&-b/e\\0&0&1&-a/e\end{pmatrix}\text{.}

The eigenvalues are obtained by a call to :meth:lapackeig.dhseqr <naginterfaces.library.lapackeig.dhseqr>.
Further details can be found in Further Comments <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02alf.html#fcomments>__.

To obtain the roots of a cubic equation, :meth:cubic_real can be used.

.. _c02al-py2-py-references:

**References**
Golub, G H and Van Loan, C F, 1996, Matrix Computations, (3rd Edition), Johns Hopkins University Press, Baltimore
"""
raise NotImplementedError

[docs]def cubic_complex(u, r, s, t):
r"""
cubic_complex determines the roots of a cubic equation with complex coefficients.

.. _c02am-py2-py-doc:

For full information please refer to the NAG Library document for c02am

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02amf.html

.. _c02am-py2-py-parameters:

**Parameters**
**u** : complex
:math:u, the coefficient of :math:z^3.

**r** : complex
:math:r, the coefficient of :math:z^2.

**s** : complex
:math:s, the coefficient of :math:z.

**t** : complex
:math:t, the constant coefficient.

**Returns**
**zero** : complex, ndarray, shape :math:\left(3\right)
:math:\mathrm{zero}[i-1] contains the :math:i\ th root.

**errest** : float, ndarray, shape :math:\left(3\right)
:math:\mathrm{errest}[i-1] contains an approximate error estimate for the :math:i\ th root.

.. _c02am-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{u} = \left(0.0, 0.0\right).

Constraint: :math:\mathrm{u}\neq \left(0.0, 0.0\right)

(errno :math:2)
The companion matrix :math:H cannot be formed without overflow.

(errno :math:3)
Failure to converge in :meth:lapackeig.zhseqr <naginterfaces.library.lapackeig.zhseqr>.

.. _c02am-py2-py-notes:

**Notes**
No equivalent traditional C interface for this routine exists in the NAG Library.

cubic_complex attempts to find the roots of the cubic equation

.. math::
uz^3+rz^2+sz+t = 0\text{,}

where :math:u, :math:r, :math:s and :math:t are complex coefficients with :math:u\neq 0.
The roots are located by finding the eigenvalues of the associated :math:3\times 3 (upper Hessenberg) companion matrix :math:H given by

.. math::
H = \begin{pmatrix}0&0&-t/u\\1&0&-s/u\\0&1&-r/u\end{pmatrix}\text{.}

The eigenvalues are obtained by a call to :meth:lapackeig.zhseqr <naginterfaces.library.lapackeig.zhseqr>.
Further details can be found in Further Comments <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02amf.html#fcomments>__.

To obtain the roots of a quadratic equation, :meth:quadratic_complex can be used.

.. _c02am-py2-py-references:

**References**
Golub, G H and Van Loan, C F, 1996, Matrix Computations, (3rd Edition), Johns Hopkins University Press, Baltimore
"""
raise NotImplementedError

[docs]def quartic_complex(e, a, b, c, d):
r"""
quartic_complex determines the roots of a quartic equation with complex coefficients.

.. _c02an-py2-py-doc:

For full information please refer to the NAG Library document for c02an

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02anf.html

.. _c02an-py2-py-parameters:

**Parameters**
**e** : complex
:math:e, the coefficient of :math:z^4.

**a** : complex
:math:a, the coefficient of :math:z^3.

**b** : complex
:math:b, the coefficient of :math:z^2.

**c** : complex
:math:c, the coefficient of :math:z.

**d** : complex
:math:d, the constant coefficient.

**Returns**
**zero** : complex, ndarray, shape :math:\left(4\right)
:math:\mathrm{zero}[i-1] contains the :math:i\ th root.

**errest** : float, ndarray, shape :math:\left(4\right)
:math:\mathrm{errest}[i-1] contains an approximate error estimate for the :math:i\ th root.

.. _c02an-py2-py-errors:

**Raises**
**NagValueError**
(errno :math:1)
On entry, :math:\mathrm{e} = \left(0.0, 0.0\right).

Constraint: :math:\mathrm{e}\neq \left(0.0, 0.0\right)

(errno :math:2)
The companion matrix :math:H cannot be formed without overflow.

(errno :math:3)
Failure to converge in :meth:lapackeig.zhseqr <naginterfaces.library.lapackeig.zhseqr>.

.. _c02an-py2-py-notes:

**Notes**
No equivalent traditional C interface for this routine exists in the NAG Library.

quartic_complex attempts to find the roots of the quartic equation

.. math::
ez^4+az^3+bz^2+cz+d = 0\text{,}

where :math:e, :math:a, :math:b, :math:c and :math:d are complex coefficients with :math:e\neq 0.
The roots are located by finding the eigenvalues of the associated :math:4\times 4 (upper Hessenberg) companion matrix :math:H given by

.. math::
H = \begin{pmatrix}0&0&0&-d/e\\1&0&0&-c/e\\0&1&0&-b/e\\0&0&1&-a/e\end{pmatrix}\text{.}

The eigenvalues are obtained by a call to :meth:lapackeig.zhseqr <naginterfaces.library.lapackeig.zhseqr>.
Further details can be found in Further Comments <https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c02/c02anf.html#fcomments>__.

To obtain the roots of a cubic equation, :meth:cubic_complex can be used.

.. _c02an-py2-py-references:

**References**
Golub, G H and Van Loan, C F, 1996, Matrix Computations, (3rd Edition), Johns Hopkins University Press, Baltimore
"""
raise NotImplementedError