Parallel Bairstow root-finding methods for real polynomials. More...

#include <utility>
#include <vector>
#include "matrix2.hpp"
#include "vector2.hpp"

Include dependency graph for rootfinding.hpp:

This graph shows which files directly or indirectly include this file:

Typedefs
using	Vec2 = ginger::Vector2< double >

using	Mat2 = ginger::Matrix2< Vec2 >

Functions
auto	initial_guess (std::vector< double > coeffs) -> std::vector< Vec2 >
	Initial guess for the parallel Bairstow method.

auto	pbairstow_even_st (const std::vector< double > &coeffs, std::vector< Vec2 > &vrs, const Options &options) -> std::pair< unsigned int, bool >
	Multi-threading Bairstow's method (even degree only)

auto	pbairstow_even_mt (const std::vector< double > &coeffs, std::vector< Vec2 > &vrs, const Options &options) -> std::pair< unsigned int, bool >

auto	pbairstow_even (const std::vector< double > &coeffs, std::vector< Vec2 > &vrs, const Options &options) -> std::pair< unsigned int, bool >

auto	horner (std::vector< double > &coeffs1, std::size_t degree, const Vec2 &vr) -> Vec2
	Horner's rule.

auto	suppress_old (Vec2 &vA, Vec2 &vA1, const Vec2 &vri, const Vec2 &vrj) -> void
	Zero suppression step in Bairstow's method (scalar arithmetic)

auto	suppress (Vec2 &vA, Vec2 &vA1, const Vec2 &vri, const Vec2 &vrj) -> void
	Zero suppression step in Bairstow's method (matrix variant)

auto	suppress2 (Vec2 &vA, Vec2 &vA1, const Vec2 &vri, const Vec2 &vrj) -> void
	Zero suppression step in Bairstow's method (variant 2)

auto	makeadjoint (const Vec2 &vr, const Vec2 &vp) -> Mat2
	Create adjoint matrix from two vectors.

auto	delta_scalar (const Vec2 &vA, const Vec2 &vr, const Vec2 &vp) -> Vec2
	Calculate Newton correction delta (scalar arithmetic, no Matrix2 temporaries)

auto	delta (const Vec2 &vA, const Vec2 &vr, const Vec2 &vp) -> Vec2
	Calculate Newton correction delta (matrix-based)

auto	horner_eval (std::vector< double > coeffs1, std::size_t degree, const double &z) -> double

auto	poly_from_quadratic_factors (const std::vector< Vec2 > &vrs) -> std::vector< double >
	Reconstruct a monic polynomial from its quadratic factors.

Detailed Description

Parallel Bairstow root-finding methods for real polynomials.

Typedef Documentation

◆ Mat2

using Mat2 = ginger::Matrix2<Vec2>

◆ Vec2

using Vec2 = ginger::Vector2<double>

Function Documentation

◆ delta()

auto delta	(	const Vec2 &	vA,
		const Vec2 &	vr,
		const Vec2 &	vp
	)		-> Vec2

inline

Calculate Newton correction delta (matrix-based)

Uses the adjoint matrix to compute the correction step in Bairstow's method:

\[ \begin{bmatrix} \Delta r \\ \Delta q \end{bmatrix} = -\frac{\operatorname{adj}(J)}{\det(J)} \, vA \]

Parameters

[in]	vA	A vector of type Vec2.
[in]	vr	A vector representing the direction of rotation.
[in]	vp	The parameter `vp` is a `Vec2` object that is passed by rvalue reference.

Returns: a Vec2 object.

◆ delta_scalar()

auto delta_scalar	(	const Vec2 &	vA,
		const Vec2 &	vr,
		const Vec2 &	vp
	)		-> Vec2

inline

Calculate Newton correction delta (scalar arithmetic, no Matrix2 temporaries)

Pure scalar computation of the Newton correction step in Bairstow's method. Avoids Matrix2 and Vector2 template instantiations in the hot loop.

Adj = [[s, -p], [-p*q, p*r + s]], det = s*(p*r + s) - p*p*q result = (Adj * vA) / det

Parameters

[in]	vA	Residual vector (A, B)
[in]	vr	Current root estimate (r, q)
[in]	vp	Suppression vector (p, s)

Returns: Vec2 Correction delta

◆ horner()

auto horner	(	std::vector< double > &	coeffs1,
		std::size_t	degree,
		const Vec2 &	vr
	)		-> Vec2

extern

Horner's rule.

Horner's rule is a method for evaluating a polynomial at a given point \(x\). It rewrites the polynomial as nested multiplication:

\[ P(x) = a_0 + x(a_1 + x(a_2 + \cdots + x(a_{n-1} + x a_n)\cdots)) \]

This allows evaluation using \(n\) multiplications and \(n\) additions.

Parameters

[in,out]	coeffs1	coeffs1 is a reference to a vector of doubles. It is used to store the coefficients of a polynomial.
[in]	degree	The parameter `degree` represents the size of the vector `coeffs1`. It indicates the number of elements in the vector `coeffs1`.
[in]	vr	vr is a Vec2 object, which represents a 2D vector. It has two components, vr.x() and vr.y(), which are used in the calculations inside the horner function.

Returns: a Vec2 object.

◆ horner_eval()

auto horner_eval	(	std::vector< double >	coeffs1,
		std::size_t	degree,
		const double &	z
	)		-> double

inline

The function horner_eval evaluates a polynomial using Horner's method.

Evaluates:

\[ P(z) = a_0 + a_1 z + a_2 z^2 + \cdots + a_n z^n \]

using the recurrence \(b_0 = a_0, \; b_{k+1} = a_{k+1} + b_k z\), returning \(b_n\).

Parameters

[in,out]	coeffs1	A vector of coefficients for a polynomial, where the coefficient at index i corresponds to the term with degree i.
[in]	degree	The degree parameter represents the degree of the polynomial. It indicates the highest power of the variable in the polynomial equation.
[in]	z	The parameter `z` is a constant value that is used as the input to the polynomial function being evaluated.

Returns: a double value.

◆ initial_guess()

auto initial_guess ( std::vector< double > coeffs ) -> std::vector< Vec2 >

extern

Initial guess for the parallel Bairstow method.

The initial_guess function calculates the initial values for the parallel Bairstow method for finding the roots of a real polynomial.

The initial quadratic factors are:

\[ x^2 - r_k x - q_k, \quad r_k = 2R\cos(2\pi\phi_2(k)), \quad q_k = -R^2, \quad k = 0,\dots,\lfloor n/2\rfloor \]

where \( R \) is estimated from the polynomial coefficients.

Parameters

[in] coeffs coeffs is a vector of doubles that represents the coefficients of a polynomial. The vector is passed by value.

Returns: The function initial_guess returns a vector of Vec2 objects.

◆ makeadjoint()

auto makeadjoint	(	const Vec2 &	vr,
		const Vec2 &	vp
	)		-> Mat2

inline

Create adjoint matrix from two vectors.

Computes the adjugate matrix of the Jacobian in Bairstow's method:

\[ \operatorname{adj}(J) = \begin{bmatrix} s & -p \cdot r_y \\ -p & p \cdot r_x + s \end{bmatrix} \]

where \( (p, s) = vp \) and \( (r_x, r_y) = vr \).

Parameters

[in]	vr	A constant reference to a Vec2 object, representing the vector vr.
[in]	vp	vp is a vector with two components, vp.x() and vp.y().

Returns: a Mat2 object.

◆ pbairstow_even()

auto pbairstow_even	(	const std::vector< double > &	coeffs,
		std::vector< Vec2 > &	vrs,
		const Options &	options
	)		-> std::pair<unsigned int, bool>

inline

◆ pbairstow_even_mt()

auto pbairstow_even_mt	(	const std::vector< double > &	coeffs,
		std::vector< Vec2 > &	vrs,
		const Options &	options
	)		-> std::pair< unsigned int, bool >

extern

◆ pbairstow_even_st()

auto pbairstow_even_st	(	const std::vector< double > &	coeffs,
		std::vector< Vec2 > &	vrs,
		const Options &	options
	)		-> std::pair< unsigned int, bool >

extern

Multi-threading Bairstow's method (even degree only)

The pbairstow_even function implements Bairstow's method for finding the roots of a real polynomial with an even degree using multi-threading.

Each thread handles one quadratic factor \( x^2 - r_i x - q_i \), applying:

\[ \begin{bmatrix} \Delta r_i \\ \Delta q_i \end{bmatrix} = -J_i^{-1} \begin{bmatrix} P_i \\ Q_i \end{bmatrix} \]

where \( P_i, Q_i \) are the remainders from synthetic division.

Parameters

[in]	coeffs	The `coeffs` parameter is a vector representing the coefficients of the polynomial. Each element of the vector corresponds to the coefficient of a term in the polynomial, starting from the highest degree term and ending with the constant term. For example, if the polynomial is `3x^2 + 2x + 1`, the coefficients vector would be `{3, 2, 1}`.
[in,out]	vrs	`vrs` is a vector of iterates, which represents the initial guesses for the roots of the polynomial. Bairstow's method will update these iterates iteratively until the desired tolerance is reached or the maximum number of iterations is reached.
[in]	options	The `options` parameter is an object of type `Options` which contains the maximum number of iterations (`max_iters`) and the tolerance (`tolerance`). These options are used to control the convergence criteria for Bairstow's method.

Returns: The function pbairstow_even returns a std::pair<unsigned int, bool>. The first element of the pair represents the number of iterations performed, and the second element represents whether the method converged to a solution within the specified tolerance.

◆ poly_from_quadratic_factors()

auto poly_from_quadratic_factors ( const std::vector< Vec2 > & vrs ) -> std::vector< double >

extern

Reconstruct a monic polynomial from its quadratic factors.

Given the quadratic factors found by Bairstow's method (each representing \(x^2 - r_i x - q_i\)), multiply them together:

\[ P(x) = \prod_{i=1}^{n/2} (x^2 - r_i x - q_i) = x^n + a_{n-1}x^{n-1} + \cdots + a_0 \]

To get the original polynomial, multiply the result by the original leading coefficient.

Parameters

[in] vrs Vector of quadratic factors, each as a Vec2 with x() = r, y() = q

Returns: std::vector<double> Monic polynomial coefficients (highest degree first)

◆ suppress()

auto suppress	(	Vec2 &	vA,
		Vec2 &	vA1,
		const Vec2 &	vri,
		const Vec2 &	vrj
	)		-> void

extern

Zero suppression step in Bairstow's method (matrix variant)

Finds coefficients of the linear remainder of a deflated polynomial without explicitly constructing the deflated polynomial, avoiding complex arithmetic within iterations.

Uses adjoint-matrix formulation.

Parameters

[in,out]	vA	First remainder coefficient
[in,out]	vA1	Second remainder coefficient
[in]	vri	First known factor
[in]	vrj	Second known factor

◆ suppress2()

auto suppress2	(	Vec2 &	vA,
		Vec2 &	vA1,
		const Vec2 &	vri,
		const Vec2 &	vrj
	)		-> void

extern

Zero suppression step in Bairstow's method (variant 2)

Alternative formulation of the zero suppression technique using adjugate matrix-based update:

\[ \begin{bmatrix} vA \\ vA_1 \end{bmatrix}^+ = \begin{bmatrix} vA \\ vA_1 \end{bmatrix} - \frac{\operatorname{adj}(J)}{\det(J)} \, vA \]

Parameters

[in,out]	vA	First remainder coefficient
[in,out]	vA1	Second remainder coefficient
[in]	vri	First known factor
[in]	vrj	Second known factor

◆ suppress_old()

auto suppress_old	(	Vec2 &	vA,
		Vec2 &	vA1,
		const Vec2 &	vri,
		const Vec2 &	vrj
	)		-> void

extern

Zero suppression step in Bairstow's method (scalar arithmetic)

In-place scalar variant that avoids Matrix2/Vector2 temporaries. Directly computes the adjusted remainder coefficients using Cramer's rule on the 2×2 Jacobian system.

Parameters

[in,out]	vA	First remainder coefficient
[in,out]	vA1	Second remainder coefficient
[in]	vri	First known factor
[in]	vrj	Second known factor

Typedefs

Functions

Detailed Description

Typedef Documentation

◆ Mat2

◆ Vec2

Function Documentation

◆ delta()

◆ delta_scalar()

◆ horner()

◆ horner_eval()

◆ initial_guess()

◆ makeadjoint()

◆ pbairstow_even()

◆ pbairstow_even_mt()

◆ pbairstow_even_st()

◆ poly_from_quadratic_factors()

◆ suppress()

◆ suppress2()

◆ suppress_old()