rootfind

package
v0.8.0 Latest Latest
Warning

This package is not in the latest version of its module.

 Go to latest Published: Jul 19, 2026 License: MIT Imports: 6 Imported by: 0

Documentation

Overview

Package rootfind implements polynomial and nonlinear root finding in pure Go.

The package is organized around two polynomial representations and a family of solvers that operate on them or on arbitrary scalar functions.

Polynomial representation

Polynomials are stored as coefficient slices in ascending order of power: the slice p represents the polynomial 

p(x) = p[0] + p[1]*x + p[2]*x^2 + ... + p[n]*x^n

so that p[i] is the coefficient of x^i and the last element is the leading coefficient. Poly holds real coefficients ([]float64) and CPoly holds complex coefficients ([]complex128). This ordering makes indexing by power trivial and matches Horner evaluation from the high-order end.

Global polynomial solvers

Several methods find all roots of a polynomial at once. DurandKerner (the Weierstrass method) and AberthEhrlich are simultaneous iterations that converge to every complex root in parallel; Bairstow extracts real quadratic factors to recover complex-conjugate pairs using only real arithmetic; and CompanionEigenvalues finds roots as the eigenvalues of the companion matrix via the Francis double-shift QR algorithm. PolyRoots is a high-level convenience that returns all roots, and RealRoots filters those with negligible imaginary part.

Real-root theory

SturmSequence and the counting routines built on it give the exact number of distinct real roots in any interval, from which IsolateRoots produces disjoint bracketing intervals that are then refined to full precision. The classical sign-based theorems are provided too: DescartesRuleOfSigns, BudanFourierCount, and a collection of a-priori root-magnitude bounds such as CauchyBound, LagrangeBound, and FujiwaraBound.

Scalar solvers

For a general continuous function of one variable the package offers bracketing methods (Bisection, FalsePosition, Brent, Ridders), open methods (Secant, Steffensen, FixedPoint), and derivative-based methods (Newton, Halley). Laguerre is a globally very reliable polynomial solver. Each returns a Result recording the located root, the residual, the iteration count, and whether convergence was achieved.

Multiplicity and deflation

gcd(p, p') isolates repeated factors: SquareFree returns the squarefree part, SquareFreeFactorization performs Yun's algorithm to group roots by multiplicity, and Multiplicity measures the order of a specific root. Roots may be removed one at a time with Poly.DeflateReal or CPoly.Deflate.

All algorithms use only the Go standard library and perform no allocation beyond what their results require.

Index

Examples

Constants

View Source 
const DefaultMaxIter = 200

DefaultMaxIter is the default iteration budget used by scalar solvers when the caller passes a non-positive maximum.

View Source 
const DefaultRootTol = 1e-14

DefaultRootTol is the default convergence tolerance for the simultaneous polynomial root solvers.

View Source 
const DefaultTol = 1e-12

DefaultTol is the default absolute tolerance used by scalar solvers when the caller passes a non-positive tolerance.

Variables

View Source 
var ErrBadInput = errors.New("rootfind: invalid input")

ErrBadInput is returned for malformed arguments that are not covered by a more specific sentinel error, such as a negative tolerance or iteration count.

View Source 
var ErrDegreeTooLow = errors.New("rootfind: polynomial degree too low for this operation")

ErrDegreeTooLow is returned when a routine requires a polynomial of at least a certain degree (for example a nonconstant polynomial) and it is not met.

View Source 
var ErrEmptyInterval = errors.New("rootfind: invalid or empty interval")

ErrEmptyInterval is returned by interval routines when the supplied interval is degenerate or has its endpoints reversed.

View Source 
var ErrNoBracket = errors.New("rootfind: endpoints do not bracket a root")

ErrNoBracket is returned by a bracketing solver when the supplied endpoints do not straddle a root, i.e. the function has the same sign at both ends.

View Source 
var ErrNoConvergence = errors.New("rootfind: iteration did not converge")

ErrNoConvergence is returned by an iterative solver that failed to reach the requested tolerance within its iteration budget.

View Source 
var ErrZeroDerivative = errors.New("rootfind: derivative vanished during iteration")

ErrZeroDerivative is returned by a derivative-based solver when the derivative vanishes at the current iterate, preventing a well-defined update step.

View Source 
var ErrZeroPolynomial = errors.New("rootfind: operation undefined for the zero polynomial")

ErrZeroPolynomial is returned when an operation that is undefined on the zero polynomial (such as division or making a polynomial monic) receives one.

Functions

func AberthEhrlich 

func AberthEhrlich(c CPoly, tol float64, maxIter int) ([]complex128, int, error)

AberthEhrlich finds all complex roots of a polynomial simultaneously using the Aberth-Ehrlich iteration, a third-order-per-step refinement of Durand-Kerner that uses the logarithmic derivative p'/p. It is typically faster and more robust than Durand-Kerner, converging in far fewer iterations for high-degree polynomials.

func AllRootsReal 

func AllRootsReal(p Poly, tol float64) bool

AllRootsReal reports whether every root of p is real, i.e. the count of distinct real roots accounts for the full degree once multiplicity is included.

func AnnulusBounds 

func AnnulusBounds(p Poly) (lo, hi float64)

AnnulusBounds returns a lower and an upper bound r0 <= |z| <= r1 that enclose every nonzero root z of p in a complex annulus. When 0 is a root the lower bound is 0.

func Bairstow 

func Bairstow(p Poly, u0, v0, tol float64, maxIter int) (u, v float64, iters int, err error)

Bairstow extracts a real quadratic factor x^2 + u*x + v from the real polynomial p, starting from initial guesses u0, v0, using Bairstow's method. It returns the coefficients u and v of the converged factor, from which the corresponding pair of (possibly complex-conjugate) roots follows via QuadraticRoots. Bairstow's method finds complex roots of a real polynomial using only real arithmetic.

func BairstowRoots 

func BairstowRoots(p Poly, tol float64, maxIter int) ([]complex128, error)

BairstowRoots finds all roots of the real polynomial p by repeatedly applying Bairstow to peel off quadratic factors and deflating, dropping to a linear solve when a single degree remains. It returns every root (real roots have zero imaginary part) using only real arithmetic internally. Roots are polished against the original polynomial.

func BracketOutward 

func BracketOutward(f Func, a, b, factor float64, maxIter int) (lo, hi float64, err error)

BracketOutward searches outward from the interval [a, b] by geometrically expanding it until f changes sign across the endpoints, returning the widened bracket. It is a convenience for seeding bracketing solvers. It returns ErrNoBracket if no sign change is found within maxIter expansions.

func BudanFourierCount 

func BudanFourierCount(p Poly, a, b float64) int

BudanFourierCount returns V(a) - V(b), where V is the sign-variation count of the Fourier sequence of p. By the Budan-Fourier theorem this value is an upper bound, differing from the true count of real roots of p in (a, b] (with multiplicity) by a nonnegative even integer, and has the same parity as that count. The endpoints must satisfy a < b.

func BudanFourierUpperBound 

func BudanFourierUpperBound(p Poly) int

BudanFourierUpperBound returns an upper bound on the number of real roots of p (with multiplicity) in the open positive axis (0, +inf), via the Fourier sequence evaluated at 0 and at a value beyond all roots.

func CPolyRoots 

func CPolyRoots(c CPoly) ([]complex128, error)

CPolyRoots returns all complex roots of the complex polynomial c using the Aberth-Ehrlich method, sorted by real then imaginary part.

func CauchyBound 

func CauchyBound(p Poly) float64

CauchyBound returns Cauchy's bound: every (real or complex) root z of p satisfies |z| <= 1 + max_{i<n} |a_i / a_n|, where a_n is the leading coefficient. It returns 0 for constant polynomials.

func CoeffSignVariations 

func CoeffSignVariations(p Poly) int

CoeffSignVariations returns the number of sign changes among the coefficients of p taken in order of ascending power (zeros skipped). By Descartes' rule of signs this is an upper bound, of the correct parity, on the number of positive real roots of p counted with multiplicity.

func CompanionEigenvalues 

func CompanionEigenvalues(p Poly) ([]complex128, error)

CompanionEigenvalues returns all roots of the polynomial p as the eigenvalues of its companion matrix, computed with the Francis double-shift QR algorithm on the real Hessenberg form (the classic EISPACK hqr routine). This is a robust, allocation-light way to obtain every complex root without complex arithmetic in the iteration. Roots are returned sorted by real then imaginary part.

func CompanionMatrix 

func CompanionMatrix(p Poly) ([][]float64, error)

CompanionMatrix returns the companion matrix of the polynomial p as a dense row-major n-by-n slice, where n is the degree of p. The characteristic polynomial of the returned matrix equals p made monic, so its eigenvalues are exactly the roots of p. The matrix is in upper-Hessenberg form with a unit subdiagonal. It returns ErrDegreeTooLow for constant polynomials.

func CountComplexRoots 

func CountComplexRoots(c CPoly) int

CountComplexRoots returns the number of roots of c counted with multiplicity, which by the fundamental theorem of algebra equals its degree.

func CountRealRoots 

func CountRealRoots(p Poly) int

CountRealRoots returns the total number of distinct real roots of p (each counted once regardless of multiplicity) using Sturm's theorem over the whole real line.

func CountRealRootsInInterval 

func CountRealRootsInInterval(p Poly, a, b float64) int

CountRealRootsInInterval returns the number of distinct real roots of p in the half-open interval (a, b].

func DescartesNegativeBound 

func DescartesNegativeBound(p Poly) int

DescartesNegativeBound returns the Descartes upper bound on the number of negative real roots of p, obtained by applying the rule of signs to p(-x).

func DescartesPositiveBound 

func DescartesPositiveBound(p Poly) int

DescartesPositiveBound returns the Descartes upper bound on the number of positive real roots of p (counted with multiplicity): the sign-variation count of its coefficient sequence. The true number of positive roots equals this value minus a nonnegative even integer.

func DescartesRuleOfSigns 

func DescartesRuleOfSigns(p Poly) (positive, negative int)

DescartesRuleOfSigns returns both Descartes bounds at once: an upper bound on the number of positive real roots and an upper bound on the number of negative real roots, each of the correct parity and counted with multiplicity.

Example 
// x^2 - 1 has one sign change -> one positive, one negative root.
p := NewPoly(-1, 0, 1)
pos, neg := DescartesRuleOfSigns(p)
fmt.Println(pos, neg)

Output:
1 1

func DiscriminantCubic 

func DiscriminantCubic(a, b, c, d float64) float64

DiscriminantCubic returns the discriminant of the cubic a*x^3+b*x^2+c*x+d, 18abcd - 4b^3d + b^2c^2 - 4ac^3 - 27a^2d^2. A positive value means three distinct real roots, zero a repeated root, and negative one real and two complex-conjugate roots.

func DiscriminantQuadratic 

func DiscriminantQuadratic(a, b, c float64) float64

DiscriminantQuadratic returns the discriminant b^2 - 4ac of a*x^2 + b*x + c. Its sign classifies the roots: positive means two real roots, zero a repeated real root, negative a complex-conjugate pair.

func DistinctRealRoots 

func DistinctRealRoots(p Poly, tol float64) []float64

DistinctRealRoots returns the distinct real roots of p, each listed once, refined to tolerance tol. It is SturmRealRoots under a descriptive name and discards multiplicity information.

func DurandKerner 

func DurandKerner(c CPoly, tol float64, maxIter int) ([]complex128, int, error)

DurandKerner finds all n complex roots of a degree-n polynomial simultaneously using the Durand-Kerner (Weierstrass) iteration. Every estimate is refined by 

z_k <- z_k - p(z_k) / prod_{j!=k} (z_k - z_j)

The method converges from the standard spiral initialization for essentially all polynomials. It returns the roots, the number of iterations used, and an error when the iteration budget is exhausted before convergence.

func DurandKernerWithInit 

func DurandKernerWithInit(c CPoly, init []complex128, tol float64, maxIter int) ([]complex128, int, error)

DurandKernerWithInit runs the Durand-Kerner iteration from a caller-supplied set of initial guesses, one per root. This is useful when good approximate roots are already known, for example when polishing the output of another method. The number of guesses must equal the degree of c.

func FindBrackets 

func FindBrackets(f Func, a, b float64, n int) [][2]float64

FindBrackets scans the interval [a, b] on a uniform grid of n subintervals and returns every subinterval across which f changes sign. It is a simple way to locate multiple roots of a function before refining each with a bracketing solver.

func FujiwaraBound 

func FujiwaraBound(p Poly) float64

FujiwaraBound returns Fujiwara's bound on the root moduli: |z| <= 2 * max_i ( |a_{n-i}/a_n| ^ (1/i) ), which is one of the sharpest simple upper bounds available.

func HalleyComplex 

func HalleyComplex(c CPoly, x0 complex128, tol float64, maxIter int) (complex128, int, error)

HalleyComplex finds a root of the complex polynomial c from x0 using Halley's cubically convergent method in complex arithmetic.

func Horner 

func Horner(coeffs []float64, x float64) float64

Horner evaluates the polynomial given by ascending-order real coefficients at x using Horner's method, without constructing a Poly.

func HornerComplex 

func HornerComplex(coeffs []complex128, x complex128) complex128

HornerComplex evaluates the polynomial given by ascending-order complex coefficients at x using Horner's method.

func IsolateRoots 

func IsolateRoots(p Poly) [][2]float64

IsolateRoots returns a set of disjoint intervals, each containing exactly one distinct real root of p, covering all real roots. Bisection driven by Sturm root counts subdivides an interval known to contain every real root until each piece is isolating. The returned intervals are sorted in increasing order.

func IsolateRootsInterval 

func IsolateRootsInterval(p Poly, a, b float64) [][2]float64

IsolateRootsInterval returns disjoint isolating intervals for the distinct real roots of p that lie in (a, b], using recursive Sturm bisection. Each returned interval contains exactly one root.

func KojimaBound 

func KojimaBound(p Poly) float64

KojimaBound returns Kojima's bound on the root moduli, formed from the ratios of consecutive coefficients: |z| <= 2 * max_i q_i, where q_1 = |a_{n-1}/a_n| and q_i = |a_{n-i}/a_{n-i+1}| for i>1 (terms with a zero denominator are skipped).

func LagrangeBound 

func LagrangeBound(p Poly) float64

LagrangeBound returns Lagrange's bound: |z| <= max(1, sum_{i<n} |a_i/a_n|) for every root z of p. It is often tighter than CauchyBound.

func Laguerre 

func Laguerre(c CPoly, x0 complex128, tol float64, maxIter int) (complex128, int, error)

Laguerre finds a single root of the polynomial c starting from x0 using Laguerre's method, which has cubic convergence for simple roots and is famously reliable: it converges to some root from almost any starting point, even for complex roots of a real polynomial. It operates in complex arithmetic throughout.

func LaguerreRoots 

func LaguerreRoots(c CPoly, tol float64, maxIter int) ([]complex128, error)

LaguerreRoots finds all roots of the polynomial c by repeated application of Laguerre followed by deflation: each root found is divided out and the search continues on the reduced polynomial. Roots are polished on the original polynomial to undo deflation error. This is the classic reliable driver used by numerical libraries for arbitrary complex polynomials.

func LowerRootBound 

func LowerRootBound(p Poly) float64

LowerRootBound returns a positive lower bound on the moduli of the nonzero roots of p: no nonzero root satisfies |z| < LowerRootBound(p). It is obtained by applying CauchyBound to the reversed polynomial and taking the reciprocal. It returns 0 when 0 is a root or p is constant.

func MaxResidual 

func MaxResidual(c CPoly, roots []complex128) float64

MaxResidual returns the largest modulus |c(z)| over the given candidate roots, a scale-free measure of how well a computed root set satisfies the polynomial.

func Multiplicity 

func Multiplicity(p Poly, r, tol float64) int

Multiplicity returns the multiplicity of the real value r as a root of p, i.e. the largest k with (x-r)^k dividing p. It uses the derivative test: r has multiplicity k when p, p', ... , p^(k-1) all vanish at r but p^(k) does not. Values are compared against a scaled tolerance tol. A return of 0 means r is not a root.

func MultiplicityComplex 

func MultiplicityComplex(c CPoly, r complex128, tol float64) int

MultiplicityComplex returns the multiplicity of the complex value r as a root of the complex polynomial c, via the complex derivative test with tolerance tol.

func NewtonComplex 

func NewtonComplex(c CPoly, x0 complex128, tol float64, maxIter int) (complex128, int, error)

NewtonComplex finds a root of the complex polynomial c from the start x0 using Newton's method in complex arithmetic. It is the complex analogue of Newton and converges quadratically near a simple root.

func PolishComplexRoots 

func PolishComplexRoots(c CPoly, roots []complex128, tol float64, steps int) []complex128

PolishComplexRoots refines each approximate root in roots with a few Newton steps on the complex polynomial c, improving accuracy after a global solve or deflation. The input slice is not modified; a refined copy is returned.

func PolyRoots 

func PolyRoots(p Poly) ([]complex128, error)

PolyRoots returns all complex roots of the real polynomial p using the Aberth-Ehrlich method, the recommended general-purpose driver. The result is sorted by real part then imaginary part. It returns ErrDegreeTooLow for constant polynomials.

Example 
// x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3)
p := NewPoly(-6, 11, -6, 1)
roots, _ := PolyRoots(p)
for _, r := range roots {
	fmt.Printf("%.4f\n", real(r))
}

Output:
1.0000
2.0000
3.0000

func QuadraticRoots 

func QuadraticRoots(a, b, c float64) (complex128, complex128)

QuadraticRoots returns the two roots of the quadratic a*x^2 + b*x + c using a numerically stable formula that avoids catastrophic cancellation. The roots are returned as complex numbers; real roots have zero imaginary part.

func RealRootInterval 

func RealRootInterval(p Poly) (lo, hi float64)

RealRootInterval returns a symmetric interval [-b, b] that contains all real roots of p, where b is the Lagrange bound. Every real root x satisfies -b <= x <= b.

func RealRoots 

func RealRoots(p Poly, imagTol float64) ([]float64, error)

RealRoots returns the real roots of p, that is the roots whose imaginary part is negligible relative to imagTol. Each is refined by a Newton polish on the real polynomial. The result is sorted in ascending order.

func SeparateRoots 

func SeparateRoots(roots []complex128, imagTol float64) (reals []float64, complexes []complex128)

SeparateRoots partitions a list of complex roots into those that are effectively real (|imag| <= imagTol) and those that are genuinely complex. The real parts of the real roots are returned sorted in ascending order.

func SignChange 

func SignChange(f Func, a, b float64) bool

SignChange reports whether f(a) and f(b) have strictly opposite signs, which guarantees a root of a continuous f in (a, b) by the intermediate value theorem. A zero at an endpoint counts as a sign change.

func SignVariations 

func SignVariations(vals []float64) int

SignVariations returns the number of sign changes in the sequence of values, skipping zeros. It is the basic count underlying Descartes' rule of signs and the Budan-Fourier theorem.

func SolveCubic 

func SolveCubic(a, b, c, d float64) ([]complex128, error)

SolveCubic returns the three roots of a*x^3 + b*x^2 + c*x + d as complex numbers, using Cardano's method with complex cube roots so the formula is numerically valid in all cases. Real roots have (near) zero imaginary part. It returns ErrDegreeTooLow when a is zero.

func SolveCubicReal 

func SolveCubicReal(a, b, c, d float64) ([]float64, error)

SolveCubicReal returns the real roots of the cubic a*x^3+b*x^2+c*x+d in ascending order, extracted from SolveCubic by keeping roots whose imaginary part is negligible.

func SolveLinear 

func SolveLinear(a, b float64) (float64, error)

SolveLinear returns the root of the linear equation a*x + b = 0. It returns ErrDegreeTooLow when a is zero.

func SolveQuadraticReal 

func SolveQuadraticReal(a, b, c float64) []float64

SolveQuadraticReal returns the real roots of a*x^2 + b*x + c in ascending order: two values when the discriminant is positive, one (repeated) when it is zero, and none when it is negative.

func SortComplex 

func SortComplex(z []complex128) []complex128

SortComplex returns a copy of z sorted by real part then imaginary part.

func SturmCountRoots 

func SturmCountRoots(seq []Poly, a, b float64) int

SturmCountRoots returns the number of distinct real roots of the polynomial underlying the Sturm sequence seq that lie in the half-open interval (a, b], computed as V(a) - V(b). The endpoints must satisfy a < b.

func SturmRealRoots 

func SturmRealRoots(p Poly, tol float64) []float64

SturmRealRoots returns all distinct real roots of p refined to tolerance tol. It isolates each root with IsolateRoots and refines it with Sturm-count bisection, so it is correct even in the presence of repeated (including even-multiplicity) roots. The roots are returned in increasing order.

Example 
// (x+2)(x-1)(x-3): three distinct real roots.
p := FromRoots(-2, 1, 3)
for _, r := range SturmRealRoots(p, 1e-12) {
	fmt.Printf("%.4f\n", r)
}

Output:
-2.0000
1.0000
3.0000

func SturmRefine 

func SturmRefine(seq []Poly, a, b, tol float64) float64

SturmRefine refines an isolating interval [a, b] known to contain exactly one distinct real root of p down to width tol, using Sturm-count bisection. Unlike sign-based bisection it works for roots of even multiplicity, where p does not change sign. It returns the interval midpoint.

func SturmVariations 

func SturmVariations(seq []Poly, x float64) int

SturmVariations returns the number of sign changes in the Sturm sequence seq evaluated at x, ignoring terms that evaluate to zero. This is the function V(x) of Sturm's theorem.

func SturmVariationsAtNegInf 

func SturmVariationsAtNegInf(seq []Poly) int

SturmVariationsAtNegInf returns the sign-variation count of the Sturm sequence as x -> -infinity, determined from the leading coefficients and degrees. It is used to count roots on unbounded left intervals.

func SturmVariationsAtPosInf 

func SturmVariationsAtPosInf(seq []Poly) int

SturmVariationsAtPosInf returns the sign-variation count of the Sturm sequence as x -> +infinity, determined from the signs of the leading coefficients.

func TotalRealRoots 

func TotalRealRoots(p Poly, tol float64) int

TotalRealRoots returns the number of real roots of p counted with multiplicity, summing the multiplicities reported by RealRootsWithMultiplicity.

Types

type CPoly 

type CPoly []complex128

CPoly is a complex polynomial stored as coefficients in ascending order of power: c[i] is the coefficient of x^i and c[len(c)-1] is the leading coefficient.

func CFromRoots 

func CFromRoots(roots ...complex128) CPoly

CFromRoots builds the monic complex polynomial with the given roots, the product (x - r0)(x - r1)... .

func NewCPoly 

func NewCPoly(coeffs ...complex128) CPoly

NewCPoly returns a CPoly built from the given ascending-order coefficients. The arguments are copied.

func (CPoly) Add 

func (c CPoly) Add(q CPoly) CPoly

Add returns the sum c+q.

func (CPoly) At 

func (c CPoly) At(x complex128) complex128

At is an alias for Eval, evaluating c(x).

func (CPoly) Clone 

func (c CPoly) Clone() CPoly

Clone returns an independent copy of c.

func (CPoly) Coeff 

func (c CPoly) Coeff(i int) complex128

Coeff returns the coefficient of x^i, or 0 when i is out of range.

func (CPoly) Deflate 

func (c CPoly) Deflate(r complex128) (quo CPoly, remainder complex128)

Deflate divides c by the linear factor (x - r), returning the quotient and the remainder c(r) via synthetic division. When r is a root the remainder is (near) zero.

func (CPoly) Degree 

func (c CPoly) Degree() int

Degree returns the degree of c, the largest index with a nonzero coefficient, or -1 for the zero polynomial.

func (CPoly) Derivative 

func (c CPoly) Derivative() CPoly

Derivative returns the formal derivative c'(x).

func (CPoly) Eval 

func (c CPoly) Eval(x complex128) complex128

Eval evaluates c(x) using Horner's method.

func (CPoly) EvalDeriv 

func (c CPoly) EvalDeriv(x complex128) (val, deriv complex128)

EvalDeriv evaluates c(x) and c'(x) together in a single Horner sweep.

func (CPoly) EvalDeriv2 

func (c CPoly) EvalDeriv2(x complex128) (val, d1, d2 complex128)

EvalDeriv2 evaluates c, c', and c” at x in one Horner sweep.

func (CPoly) IsMonic 

func (c CPoly) IsMonic(tol float64) bool

IsMonic reports whether the leading coefficient of c equals 1 within tol.

func (CPoly) IsZero 

func (c CPoly) IsZero() bool

IsZero reports whether c is the zero polynomial.

func (CPoly) LeadingCoeff 

func (c CPoly) LeadingCoeff() complex128

LeadingCoeff returns the leading coefficient, or 0 for the zero polynomial.

func (CPoly) Monic 

func (c CPoly) Monic() (CPoly, error)

Monic returns c divided by its leading coefficient. It returns ErrZeroPolynomial when c is zero.

func (CPoly) Mul 

func (c CPoly) Mul(q CPoly) CPoly

Mul returns the product c*q by convolution.

func (CPoly) Neg 

func (c CPoly) Neg() CPoly

Neg returns -c.

func (CPoly) NumTerms 

func (c CPoly) NumTerms() int

NumTerms returns the number of nonzero coefficients of c.

func (CPoly) Pow 

func (c CPoly) Pow(k int) CPoly

Pow returns c raised to the nonnegative integer power k by repeated squaring.

func (CPoly) Scale 

func (c CPoly) Scale(s complex128) CPoly

Scale returns the polynomial s*c.

func (CPoly) String 

func (c CPoly) String() string

String renders c in descending-power notation with parenthesized complex coefficients, for example "(1+2i)x^2 + (3+0i)".

func (CPoly) Sub 

func (c CPoly) Sub(q CPoly) CPoly

Sub returns the difference c-q.

func (CPoly) ToReal 

func (c CPoly) ToReal() Poly

ToReal returns the real polynomial formed by taking the real part of each coefficient of c, discarding imaginary parts.

func (CPoly) Trim 

func (c CPoly) Trim() CPoly

Trim returns c with trailing zero coefficients removed; the underlying array is shared with c.

type ComplexRootMultiplicity 

type ComplexRootMultiplicity struct {
	// Root is the location of the root in the complex plane.
	Root complex128
	// Multiplicity is how many times Root occurs.
	Multiplicity int
}

ComplexRootMultiplicity pairs a complex root with its multiplicity.

func GroupComplexRoots 

func GroupComplexRoots(roots []complex128, tol float64) []ComplexRootMultiplicity

GroupComplexRoots clusters a list of complex roots so that roots within distance tol of one another are merged into a single representative (their centroid) carrying a multiplicity equal to the cluster size. This recovers multiplicity from the output of a numerical global solver, where a root of multiplicity m appears as a tight cluster of m nearby approximate roots.

type Func 

type Func func(x float64) float64

Func is a real scalar function of one real variable.

type Poly 

type Poly []float64

Poly is a real polynomial stored as coefficients in ascending order of power: p[i] is the coefficient of x^i, so p[len(p)-1] is the leading coefficient. The zero polynomial is represented by an empty or all-zero slice.

func DeflateRoots 

func DeflateRoots(p Poly, roots []float64) Poly

DeflateRoots divides p successively by (x - r) for each r in roots, returning the final quotient. It is the composition of repeated Poly.DeflateReal and removes the listed roots from p.

func FourierSequence 

func FourierSequence(p Poly) []Poly

FourierSequence returns the Fourier (derivative) sequence of p: p, p', p”, ... , p^(n), a chain of length deg(p)+1 whose sign variations feed the Budan-Fourier theorem.

func FromRoots 

func FromRoots(roots ...float64) Poly

FromRoots builds the monic real polynomial whose roots are exactly the given values, i.e. the product (x - r0)(x - r1)... .

func FromRootsWithLead 

func FromRootsWithLead(lead float64, roots ...float64) Poly

FromRootsWithLead builds the real polynomial lead*(x-r0)(x-r1)... with the given leading coefficient.

func NewPoly 

func NewPoly(coeffs ...float64) Poly

NewPoly returns a Poly built from the given ascending-order coefficients. The arguments are copied, so the caller may reuse the slice afterwards.

func PolyFromDesc 

func PolyFromDesc(coeffs ...float64) Poly

PolyFromDesc returns a Poly from coefficients given in descending order of power, i.e. highest-degree coefficient first, as polynomials are usually written. It is the inverse of Poly.CoeffsDesc.

func SquareFree 

func SquareFree(p Poly) (Poly, error)

SquareFree returns the squarefree part of p, that is p divided by gcd(p, p'). The result has the same roots as p but each with multiplicity one. It returns ErrZeroPolynomial for the zero polynomial and a constant for a nonzero constant input.

func SturmSequence 

func SturmSequence(p Poly) []Poly

SturmSequence returns the canonical Sturm chain of the polynomial p: 

p0 = p, p1 = p', p_{k+1} = -(p_{k-1} mod p_k)

continuing until a constant is reached. The sign-variation counts of this sequence, via SturmVariations, give the exact number of distinct real roots of p in any interval by Sturm's theorem. The zero polynomial yields an empty sequence.

func (Poly) Add 

func (p Poly) Add(q Poly) Poly

Add returns the sum p+q.

func (Poly) At 

func (p Poly) At(x float64) float64

At is an alias for Eval, evaluating p(x).

func (Poly) Clone 

func (p Poly) Clone() Poly

Clone returns an independent copy of p.

func (Poly) Coeff 

func (p Poly) Coeff(i int) float64

Coeff returns the coefficient of x^i, or 0 when i is out of range.

func (Poly) CoeffsDesc 

func (p Poly) CoeffsDesc() []float64

CoeffsDesc returns the coefficients of p in descending order of power, highest-degree first, trimmed of leading zeros. The zero polynomial yields a single 0.

func (Poly) Compose 

func (p Poly) Compose(q Poly) Poly

Compose returns the polynomial p(q(x)) formed by substituting q into p, using Horner's method over polynomial arithmetic.

func (Poly) DeflateReal 

func (p Poly) DeflateReal(r float64) (quo Poly, remainder float64)

DeflateReal divides p by the linear factor (x - r), returning the quotient and the remainder p(r). When r is an exact root the remainder is (near) zero. The division is performed by synthetic (Horner) division.

func (Poly) Degree 

func (p Poly) Degree() int

Degree returns the degree of p, that is the largest index whose coefficient is nonzero. The zero polynomial has degree -1 by convention.

func (Poly) Derivative 

func (p Poly) Derivative() Poly

Derivative returns the formal derivative p'(x).

func (Poly) DivMod 

func (p Poly) DivMod(d Poly) (q, r Poly, err error)

DivMod divides p by d and returns the quotient q and remainder r satisfying p = q*d + r with deg(r) < deg(d). It returns ErrZeroPolynomial when d is zero.

func (Poly) Equal 

func (p Poly) Equal(q Poly, tol float64) bool

Equal reports whether p and q are equal as polynomials, that is every coefficient agrees within tol after ignoring trailing zeros.

func (Poly) Eval 

func (p Poly) Eval(x float64) float64

Eval evaluates p(x) using Horner's method, which is both fast and numerically stable.

func (Poly) EvalComplex 

func (p Poly) EvalComplex(x complex128) complex128

EvalComplex evaluates p at a complex argument using Horner's method.

func (Poly) EvalDeriv 

func (p Poly) EvalDeriv(x float64) (val, deriv float64)

EvalDeriv evaluates p(x) and p'(x) simultaneously with a single Horner sweep, returning the value and the first derivative.

func (Poly) EvalDeriv2 

func (p Poly) EvalDeriv2(x float64) (val, d1, d2 float64)

EvalDeriv2 evaluates p, p', and p” at x in a single Horner sweep.

func (Poly) GCD 

func (p Poly) GCD(q Poly) Poly

GCD returns a greatest common divisor of p and q using the Euclidean algorithm, normalized to be monic. The GCD of two zero polynomials is zero.

func (Poly) InfNorm 

func (p Poly) InfNorm() float64

InfNorm returns the largest absolute coefficient of p (the sup norm of its coefficient vector).

func (Poly) Integral 

func (p Poly) Integral(c float64) Poly

Integral returns an antiderivative of p whose constant term equals c.

func (Poly) IsConstant 

func (p Poly) IsConstant() bool

IsConstant reports whether p has degree 0 or is the zero polynomial.

func (Poly) IsMonic 

func (p Poly) IsMonic(tol float64) bool

IsMonic reports whether the leading coefficient of p equals 1 within tol.

func (Poly) IsZero 

func (p Poly) IsZero() bool

IsZero reports whether p is the zero polynomial, i.e. every coefficient is 0.

func (Poly) L1Norm 

func (p Poly) L1Norm() float64

L1Norm returns the sum of absolute values of the coefficients of p.

func (Poly) L2Norm 

func (p Poly) L2Norm() float64

L2Norm returns the Euclidean norm of the coefficient vector of p.

func (Poly) LeadingCoeff 

func (p Poly) LeadingCoeff() float64

LeadingCoeff returns the coefficient of the highest-degree term, or 0 for the zero polynomial.

func (Poly) Monic 

func (p Poly) Monic() (Poly, error)

Monic returns p divided by its leading coefficient, so the result has leading coefficient 1. It returns ErrZeroPolynomial when p is zero.

func (Poly) Mul 

func (p Poly) Mul(q Poly) Poly

Mul returns the product p*q computed by convolution of the coefficient sequences.

func (Poly) Neg 

func (p Poly) Neg() Poly

Neg returns -p.

func (Poly) NumTerms 

func (p Poly) NumTerms() int

NumTerms returns the number of nonzero coefficients of p.

func (Poly) Pow 

func (p Poly) Pow(k int) Poly

Pow returns p raised to the nonnegative integer power k by repeated squaring. Pow(0) is the constant polynomial 1.

func (Poly) ProductOfRoots 

func (p Poly) ProductOfRoots() float64

ProductOfRoots returns the product of all roots of p counted with multiplicity, which by Vieta's formulas equals (-1)^n * a_0/a_n.

func (Poly) Quo 

func (p Poly) Quo(d Poly) (Poly, error)

Quo returns the quotient of p divided by d, discarding the remainder.

func (Poly) ReflectX 

func (p Poly) ReflectX() Poly

ReflectX returns the polynomial p(-x), whose positive roots are the negatives of the negative roots of p.

func (Poly) Rem 

func (p Poly) Rem(d Poly) (Poly, error)

Rem returns the remainder of p divided by d.

func (Poly) Reverse 

func (p Poly) Reverse() Poly

Reverse returns the reversal (reciprocal) polynomial x^n * p(1/x), whose coefficient sequence is that of p reversed. Its nonzero roots are the reciprocals of the nonzero roots of p.

func (Poly) Scale 

func (p Poly) Scale(s float64) Poly

Scale returns the polynomial s*p.

func (Poly) ShiftScale 

func (p Poly) ShiftScale(a, b float64) Poly

ShiftScale returns the polynomial p(a*x + b), the composition of p with the affine map x |-> a*x + b.

func (Poly) String 

func (p Poly) String() string

String renders p in conventional descending-power notation, for example "2x^3 - x + 5". The zero polynomial renders as "0".

func (Poly) Sub 

func (p Poly) Sub(q Poly) Poly

Sub returns the difference p-q.

func (Poly) SumOfRoots 

func (p Poly) SumOfRoots() float64

SumOfRoots returns the sum of all roots of p counted with multiplicity, which by Vieta's formulas equals -a_{n-1}/a_n. It returns 0 for constant polynomials.

func (Poly) ToComplex 

func (p Poly) ToComplex() CPoly

ToComplex converts a real polynomial to its complex counterpart.

func (Poly) Trim 

func (p Poly) Trim() Poly

Trim returns p with trailing (high-order) zero coefficients removed. The result always has length equal to Degree()+1, or length zero for the zero polynomial. The underlying array is shared with p.

type Result 

type Result struct {
	// Root is the best estimate of the root that was found.
	Root float64
	// Value is the function value at Root (the residual).
	Value float64
	// Iterations is the number of iterations actually performed.
	Iterations int
	// Converged reports whether the requested tolerance was reached.
	Converged bool
}

Result records the outcome of a scalar root-finding iteration.

func Bisection 

func Bisection(f Func, a, b, tol float64, maxIter int) (Result, error)

Bisection finds a root of f in the bracketing interval [a, b] by repeated halving. It requires f(a) and f(b) to have opposite signs and returns ErrNoBracket otherwise. Bisection is unconditionally convergent and halves the bracket each step.

func Brent 

func Brent(f Func, a, b, tol float64, maxIter int) (Result, error)

Brent finds a root of f in [a, b] using Brent's method, which combines bisection, the secant method, and inverse quadratic interpolation. It is the recommended general-purpose bracketing solver: robust like bisection yet usually superlinear. It requires f(a) and f(b) to have opposite signs.

Example 
// Solve cos(x) = x on [0, 1].
f := func(x float64) float64 { return math.Cos(x) - x }
res, _ := Brent(f, 0, 1, 1e-12, 100)
fmt.Printf("%.10f\n", res.Root)

Output:
0.7390851332

func FalsePosition 

func FalsePosition(f Func, a, b, tol float64, maxIter int) (Result, error)

FalsePosition (regula falsi) finds a root in [a, b] by interpolating a secant line through the bracket endpoints and keeping the subinterval that still brackets the root. It requires an initial sign change.

func FixedPoint 

func FixedPoint(g Func, x0, tol float64, maxIter int) (Result, error)

FixedPoint iterates x <- g(x) from x0 to find a fixed point of g, which is a root of f(x) = g(x) - x. Convergence requires |g'| < 1 near the fixed point.

func Halley 

func Halley(f, df, d2f Func, x0, tol float64, maxIter int) (Result, error)

Halley finds a root of f using Halley's third-order method, given the first derivative df and second derivative d2f. It converges cubically near a simple root, faster than Newton's method at the cost of a second derivative.

func Illinois 

func Illinois(f Func, a, b, tol float64, maxIter int) (Result, error)

Illinois is the Illinois variant of the false-position method: when an endpoint is retained across consecutive iterations its function value is halved, which cures the slow one-sided convergence of plain regula falsi and restores superlinear behaviour while keeping the guaranteed bracket.

func Newton 

func Newton(f, df Func, x0, tol float64, maxIter int) (Result, error)

Newton finds a root of f using Newton's method, given the derivative df. It converges quadratically near a simple root but may diverge from a poor start or where df vanishes, in which case ErrZeroDerivative is returned.

func Ridders 

func Ridders(f Func, a, b, tol float64, maxIter int) (Result, error)

Ridders finds a root of f in the bracket [a, b] using Ridders' method, which applies an exponential correction to the false-position estimate and converges quadratically while always maintaining a bracket. It requires a sign change.

func Secant 

func Secant(f Func, x0, x1, tol float64, maxIter int) (Result, error)

Secant finds a root of f from two initial guesses x0 and x1 using the secant iteration, which approximates the derivative by a finite difference. It converges superlinearly (order ~1.618) but is not guaranteed to bracket.

func Steffensen 

func Steffensen(f Func, x0, tol float64, maxIter int) (Result, error)

Steffensen finds a root of f using Steffensen's method, a derivative-free iteration that achieves quadratic convergence by estimating the derivative from f(x) and f(x+f(x)). It needs only a single starting point.

type RootMultiplicity 

type RootMultiplicity struct {
	// Root is the location of the root.
	Root float64
	// Multiplicity is how many times Root occurs as a root of p.
	Multiplicity int
}

RootMultiplicity pairs a real root value with its multiplicity in p.

func RealRootsWithMultiplicity 

func RealRootsWithMultiplicity(p Poly, tol float64) []RootMultiplicity

RealRootsWithMultiplicity returns the distinct real roots of p together with their multiplicities. Distinct roots are located by SturmRealRoots and each multiplicity is measured by the derivative test in Multiplicity. The result is sorted by root location.

type SquareFreeFactor 

type SquareFreeFactor struct {
	// Factor is a monic squarefree polynomial whose roots all have the same
	// multiplicity in the original polynomial.
	Factor Poly
	// Multiplicity is the common multiplicity of the roots of Factor.
	Multiplicity int
}

SquareFreeFactor is one factor of a squarefree factorization: the polynomial Factor collects exactly the roots that occur with the given Multiplicity.

func SquareFreeFactorization 

func SquareFreeFactorization(p Poly) ([]SquareFreeFactor, error)

SquareFreeFactorization returns the squarefree factorization of p using Yun's algorithm: it decomposes p (up to its leading constant) as a product prod_i a_i(x)^i where each a_i is monic and squarefree and the a_i are pairwise coprime. Only factors of positive degree are returned, sorted by increasing multiplicity. This groups the roots of p by their multiplicity exactly, using only gcd computations.

Source Files

View all Source files

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL
go.dev uses cookies from Google to deliver and enhance the quality of its services and to analyze traffic. Learn more.