Struct Polynomial 

pub struct Polynomial { /* private fields */ }

A multivariate polynomial over rationals. Represented as a sum of terms, sorted by monomial order.

Implementations

impl Polynomial

pub fn refine_roots( &self, var: Var, intervals: &[(BigRational, BigRational)], max_iterations: usize, tolerance: &BigRational, ) -> Vec<BigRational>

Refine all roots in the given intervals using Newton-Raphson.

Takes a list of isolating intervals (from isolate_roots) and refines each root to higher precision using Newton-Raphson iteration.

Arguments

• var - The variable to solve for
• intervals - List of isolating intervals from isolate_roots
• max_iterations - Maximum number of Newton-Raphson iterations per root
• tolerance - Convergence tolerance

Returns

Vector of refined root approximations

Examples

use num_bigint::BigInt;
use num_rational::BigRational;
use oxiz_math::polynomial::Polynomial;

// Solve x^2 - 2 = 0
let p = Polynomial::from_coeffs_int(&[
    (1, &[(0, 2)]),  // x^2
    (-2, &[]),       // -2
]);

let intervals = p.isolate_roots(0);
let tolerance = BigRational::new(BigInt::from(1), BigInt::from(1000000));
let refined_roots = p.refine_roots(0, &intervals, 10, &tolerance);

pub fn taylor_expansion( &self, var: Var, point: &BigRational, degree: u32, ) -> Polynomial

Compute the Taylor series expansion of a univariate polynomial around a point.

For a polynomial p(x), computes the Taylor series: p(x) = Σ_{k=0}^n (p^(k)(a) / k!) * (x - a)^k

where p^(k) is the k-th derivative.

Arguments

• var - The variable to expand around
• point - The point to expand around (typically 0 for Maclaurin series)
• degree - Maximum degree of the Taylor expansion

Returns

A polynomial representing the Taylor series expansion

Examples

use num_bigint::BigInt;
use num_rational::BigRational;
use oxiz_math::polynomial::Polynomial;

// Expand x^2 + 2x + 1 around x = 0
let p = Polynomial::from_coeffs_int(&[
    (1, &[(0, 2)]),  // x^2
    (2, &[(0, 1)]),  // 2x
    (1, &[]),        // 1
]);

let taylor = p.taylor_expansion(0, &BigRational::from_integer(BigInt::from(0)), 3);
// Should be the same polynomial since it's already a polynomial

pub fn maclaurin_expansion(&self, var: Var, degree: u32) -> Polynomial

Compute the Maclaurin series expansion (Taylor series around 0).

This is a convenience method for Taylor expansion around 0.

Arguments

• var - The variable to expand
• degree - Maximum degree of the expansion

Examples

use oxiz_math::polynomial::Polynomial;

// Expand x^3 - 2x around x = 0
let p = Polynomial::from_coeffs_int(&[
    (1, &[(0, 3)]),   // x^3
    (-2, &[(0, 1)]),  // -2x
]);

let maclaurin = p.maclaurin_expansion(0, 4);

pub fn bisection( &self, var: Var, lower: BigRational, upper: BigRational, max_iterations: usize, tolerance: &BigRational, ) -> Option<BigRational>

Find a root using the bisection method.

The bisection method is a robust root-finding algorithm that works by repeatedly bisecting an interval and selecting the subinterval where the function changes sign.

Arguments

• var - The variable to solve for
• lower - Lower bound (must have f(lower) and f(upper) with opposite signs)
• upper - Upper bound (must have f(lower) and f(upper) with opposite signs)
• max_iterations - Maximum number of iterations
• tolerance - Convergence tolerance (stop when interval width < tolerance)

Returns

The approximate root, or None if the initial bounds don’t bracket a root

Examples

use num_bigint::BigInt;
use num_rational::BigRational;
use oxiz_math::polynomial::Polynomial;

// Solve x^2 - 2 = 0 (find sqrt(2))
let p = Polynomial::from_coeffs_int(&[
    (1, &[(0, 2)]),  // x^2
    (-2, &[]),       // -2
]);

let lower = BigRational::from_integer(BigInt::from(1));
let upper = BigRational::from_integer(BigInt::from(2));
let tolerance = BigRational::new(BigInt::from(1), BigInt::from(1000000));

let root = p.bisection(0, lower, upper, 100, &tolerance);
assert!(root.is_some());

pub fn secant( &self, var: Var, x0: BigRational, x1: BigRational, max_iterations: usize, tolerance: &BigRational, ) -> Option<BigRational>

Find a root using the secant method.

The secant method is similar to Newton-Raphson but doesn’t require computing derivatives. It uses two previous points to approximate the derivative.

Arguments

• var - The variable to solve for
• x0 - First initial guess
• x1 - Second initial guess (should be close to x0)
• max_iterations - Maximum number of iterations
• tolerance - Convergence tolerance

Returns

The approximate root, or None if the method fails to converge

Examples

use num_bigint::BigInt;
use num_rational::BigRational;
use oxiz_math::polynomial::Polynomial;

// Solve x^2 - 2 = 0
let p = Polynomial::from_coeffs_int(&[
    (1, &[(0, 2)]),  // x^2
    (-2, &[]),       // -2
]);

let x0 = BigRational::from_integer(BigInt::from(1));
let x1 = BigRational::new(BigInt::from(3), BigInt::from(2)); // 1.5
let tolerance = BigRational::new(BigInt::from(1), BigInt::from(1000000));

let root = p.secant(0, x0, x1, 20, &tolerance);
assert!(root.is_some());

pub fn sign_at(&self, var_signs: &FxHashMap<Var, i8>) -> Option<i8>

Sign of the polynomial given signs of variables. Returns Some(1) for positive, Some(-1) for negative, Some(0) for zero, or None if undetermined.

pub fn is_irreducible(&self, var: Var) -> bool

Check if a univariate polynomial is irreducible over rationals. This is a heuristic test - returns false if definitely reducible, true if likely irreducible (but not guaranteed).

pub fn factor(&self, var: Var) -> Vec<(Polynomial, u32)>

Factor a univariate polynomial into irreducible factors. Returns a vector of (factor, multiplicity) pairs. Each factor is monic and primitive.

pub fn content(&self) -> BigRational

Content of a polynomial: GCD of all coefficients. Returns the rational content.

pub fn compose(&self, var: Var, other: &Polynomial) -> Polynomial

Compose this polynomial with another: compute p(q(x)).

This is an alias for substitute with clearer semantics for composition. If self is p(x) and other is q(x), returns p(q(x)).

pub fn lagrange_interpolate( var: Var, points: &[(BigRational, BigRational)], ) -> Option<Polynomial>

Lagrange polynomial interpolation.

Given a set of points (x_i, y_i), constructs the unique polynomial of minimal degree that passes through all the points.

Arguments

• var - The variable to use for the polynomial
• points - A slice of (x, y) coordinate pairs

Returns

The interpolating polynomial, or None if points is empty or contains duplicate x values

Reference

Standard Lagrange interpolation formula from numerical analysis textbooks

pub fn newton_interpolate( var: Var, points: &[(BigRational, BigRational)], ) -> Option<Polynomial>

Newton polynomial interpolation using divided differences.

Constructs the same interpolating polynomial as Lagrange interpolation, but using Newton’s divided difference form which can be more efficient for incremental construction.

Arguments

• var - The variable to use for the polynomial
• points - A slice of (x, y) coordinate pairs

Returns

The interpolating polynomial, or None if points is empty or contains duplicate x values

pub fn chebyshev_first_kind(var: Var, n: u32) -> Polynomial

Generate the nth Chebyshev polynomial of the first kind T_n(x).

Chebyshev polynomials of the first kind are defined by:

• T_0(x) = 1
• T_1(x) = x
• T_{n+1}(x) = 2x T_n(x) - T_{n-1}(x)

These polynomials are orthogonal on [-1, 1] with respect to the weight function 1/√(1-x²) and are useful for polynomial approximation.

Arguments

• var - The variable to use for the polynomial
• n - The degree of the Chebyshev polynomial

Returns

The nth Chebyshev polynomial T_n(var)

pub fn chebyshev_second_kind(var: Var, n: u32) -> Polynomial

Generate the nth Chebyshev polynomial of the second kind U_n(x).

Chebyshev polynomials of the second kind are defined by:

• U_0(x) = 1
• U_1(x) = 2x
• U_{n+1}(x) = 2x U_n(x) - U_{n-1}(x)

These polynomials are orthogonal on [-1, 1] with respect to the weight function √(1-x²).

Arguments

• var - The variable to use for the polynomial
• n - The degree of the Chebyshev polynomial

Returns

The nth Chebyshev polynomial U_n(var)

pub fn chebyshev_nodes(n: u32) -> Vec<BigRational>

Compute Chebyshev nodes (zeros of Chebyshev polynomial) for use in interpolation.

The Chebyshev nodes minimize the Runge phenomenon in polynomial interpolation. For degree n, returns n+1 nodes in [-1, 1].

The nodes are: x_k = cos((2k+1)π / (2n+2)) for k = 0, 1, …, n

Note: This function returns approximate rational values. For exact values, use algebraic numbers or symbolic computation.

Arguments

• n - The degree (returns n+1 nodes)

Returns

Approximations of the Chebyshev nodes as BigRational values

pub fn legendre(var: Var, n: u32) -> Polynomial

Generate the nth Legendre polynomial P_n(x).

Legendre polynomials are orthogonal polynomials on [-1, 1] with respect to the weight function 1. They are defined by the recurrence:

• P_0(x) = 1
• P_1(x) = x
• (n+1) P_{n+1}(x) = (2n+1) x P_n(x) - n P_{n-1}(x)

These polynomials are useful for Gaussian quadrature and least-squares polynomial approximation.

Arguments

• var - The variable to use for the polynomial
• n - The degree of the Legendre polynomial

Returns

The nth Legendre polynomial P_n(var)

pub fn hermite(var: Var, n: u32) -> Polynomial

Generate the nth Hermite polynomial H_n(x) (physicist’s version).

Hermite polynomials are orthogonal with respect to the weight function e^(-x²). The physicist’s version is defined by:

• H_0(x) = 1
• H_1(x) = 2x
• H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x)

These polynomials are useful in quantum mechanics, probability theory, and numerical analysis.

Arguments

• var - The variable to use for the polynomial
• n - The degree of the Hermite polynomial

Returns

The nth Hermite polynomial H_n(var)

pub fn laguerre(var: Var, n: u32) -> Polynomial

Generate the nth Laguerre polynomial L_n(x).

Laguerre polynomials are orthogonal with respect to the weight function e^(-x) on [0, ∞). They are defined by:

• L_0(x) = 1
• L_1(x) = 1 - x
• (n+1) L_{n+1}(x) = (2n+1-x) L_n(x) - n L_{n-1}(x)

These polynomials are useful in quantum mechanics and numerical analysis.

Arguments

• var - The variable to use for the polynomial
• n - The degree of the Laguerre polynomial

Returns

The nth Laguerre polynomial L_n(var)

impl Polynomial

pub fn neg(&self) -> Polynomial

Negate the polynomial.

pub fn add(&self, other: &Polynomial) -> Polynomial

Add two polynomials.

pub fn sub(&self, other: &Polynomial) -> Polynomial

Subtract two polynomials.

pub fn scale(&self, c: &BigRational) -> Polynomial

Multiply by a scalar.

pub fn mul(&self, other: &Polynomial) -> Polynomial

Multiply two polynomials. Uses Karatsuba algorithm for large univariate polynomials.

pub fn mul_monomial(&self, m: &Monomial) -> Polynomial

Multiply by a monomial.

pub fn pow(&self, k: u32) -> Polynomial

Compute p^k.

impl Polynomial

pub fn zero() -> Self

Create the zero polynomial.

pub fn one() -> Self

Create the one polynomial.

pub fn constant(c: BigRational) -> Self

Create a constant polynomial.

pub fn from_var(var: Var) -> Self

Create a polynomial from a single variable.

pub fn from_var_power(var: Var, power: u32) -> Self

Create a polynomial x^k.

pub fn from_terms( terms: impl IntoIterator<Item = Term>, order: MonomialOrder, ) -> Self

Create a polynomial from terms. Normalizes and combines like terms.

pub fn from_coeffs_int(coeffs: &[(i64, &[(Var, u32)])]) -> Self

Create a polynomial from integer coefficients.

pub fn linear(coeffs: &[(BigRational, Var)], constant: BigRational) -> Self

Create a linear polynomial a1x1 + a2x2 + … + c.

pub fn univariate(var: Var, coeffs: &[BigRational]) -> Self

Create a univariate polynomial from coefficients. coeffs[i] is the coefficient of x^i.

pub fn is_zero(&self) -> bool

Check if the polynomial is zero.

pub fn is_constant(&self) -> bool

Check if the polynomial is a non-zero constant.

pub fn constant_value(&self) -> BigRational

Get the constant value of the polynomial.

Returns the constant coefficient if the polynomial is constant, or zero if the polynomial is zero.

pub fn is_one(&self) -> bool

Check if the polynomial is one.

pub fn is_univariate(&self) -> bool

Check if the polynomial is univariate.

pub fn is_linear(&self) -> bool

Check if the polynomial is linear (all terms have degree <= 1).

pub fn num_terms(&self) -> usize

Get the number of terms.

pub fn terms(&self) -> &[Term]

Get the terms.

pub fn total_degree(&self) -> u32

Get the total degree of the polynomial.

pub fn degree(&self, var: Var) -> u32

Get the degree with respect to a specific variable.

pub fn max_var(&self) -> Var

Get the maximum variable in the polynomial, or NULL_VAR if constant.

pub fn vars(&self) -> Vec<Var>

Get all variables in the polynomial.

pub fn leading_term(&self) -> Option<&Term>

Get the leading term (with respect to monomial order).

pub fn leading_coeff(&self) -> BigRational

Get the leading coefficient.

pub fn leading_monomial(&self) -> Option<&Monomial>

Get the leading monomial.

pub fn constant_term(&self) -> BigRational

Get the constant term.

pub fn univ_coeff(&self, var: Var, k: u32) -> BigRational

Get the coefficient of x^k for a univariate polynomial.

pub fn coeff(&self, var: Var, k: u32) -> Polynomial

Get the coefficient polynomial for x^k. For polynomial p(y_1, …, y_n, x), returns coefficient of x^k.

pub fn leading_coeff_wrt(&self, var: Var) -> Polynomial

Get the leading coefficient with respect to variable x.

impl Polynomial

pub fn derivative(&self, var: Var) -> Polynomial

Compute the derivative with respect to a variable.

pub fn nth_derivative(&self, var: Var, n: u32) -> Polynomial

Compute the nth derivative of the polynomial with respect to a variable.

Arguments

• var - The variable to differentiate with respect to
• n - The order of the derivative (n = 0 returns the polynomial itself)

Examples

use oxiz_math::polynomial::Polynomial;

// f(x) = x^3 = x³
let f = Polynomial::from_coeffs_int(&[(1, &[(0, 3)])]);

// f'(x) = 3x^2
let f_prime = f.nth_derivative(0, 1);
assert_eq!(f_prime.total_degree(), 2);

// f''(x) = 6x
let f_double_prime = f.nth_derivative(0, 2);
assert_eq!(f_double_prime.total_degree(), 1);

// f'''(x) = 6 (constant)
let f_triple_prime = f.nth_derivative(0, 3);
assert_eq!(f_triple_prime.total_degree(), 0);
assert!(!f_triple_prime.is_zero());

// f''''(x) = 0
let f_fourth = f.nth_derivative(0, 4);
assert!(f_fourth.is_zero());

pub fn gradient(&self) -> Vec<Polynomial>

Computes the gradient (vector of partial derivatives) with respect to all variables.

For a multivariate polynomial f(x₁, x₂, …, xₙ), returns the vector: ∇f = [∂f/∂x₁, ∂f/∂x₂, …, ∂f/∂xₙ]

Returns

A vector of polynomials, one for each variable, ordered by variable index.

Example

use oxiz_math::polynomial::Polynomial;
// f(x,y) = x²y + 2xy + y²
let f = Polynomial::from_coeffs_int(&[
    (1, &[(0, 2), (1, 1)]), // x²y
    (2, &[(0, 1), (1, 1)]), // 2xy
    (1, &[(1, 2)]),          // y²
]);

let grad = f.gradient();
// ∂f/∂x = 2xy + 2y
// ∂f/∂y = x² + 2x + 2y
assert_eq!(grad.len(), 2);

pub fn hessian(&self) -> Vec<Vec<Polynomial>>

Computes the Hessian matrix (matrix of second-order partial derivatives).

For a multivariate polynomial f(x₁, x₂, …, xₙ), returns the symmetric matrix: H[i,j] = ∂²f/(∂xᵢ∂xⱼ)

The Hessian is useful for:

• Optimization (finding local minima/maxima)
• Convexity analysis
• Second-order Taylor approximations

Returns

A vector of vectors representing the Hessian matrix. The matrix is symmetric: H[i][j] = H[j][i]

Example

use oxiz_math::polynomial::Polynomial;
// f(x,y) = x² + xy + y²
let f = Polynomial::from_coeffs_int(&[
    (1, &[(0, 2)]),          // x²
    (1, &[(0, 1), (1, 1)]), // xy
    (1, &[(1, 2)]),          // y²
]);

let hessian = f.hessian();
// H = [[2, 1],
//      [1, 2]]
assert_eq!(hessian.len(), 2);
assert_eq!(hessian[0].len(), 2);

pub fn jacobian(polys: &[Polynomial]) -> Vec<Vec<Polynomial>>

Computes the Jacobian matrix for a vector of polynomials.

For a vector of polynomials f = (f₁, f₂, …, fₘ) each depending on variables x = (x₁, x₂, …, xₙ), the Jacobian is the m×n matrix: J[i,j] = ∂fᵢ/∂xⱼ

Arguments

• polys - Vector of polynomials representing the function components

Returns

A matrix where each row i contains the gradient of polynomial i

Example

use oxiz_math::polynomial::Polynomial;
// f₁(x,y) = x² + y
// f₂(x,y) = x + y²
let f1 = Polynomial::from_coeffs_int(&[(1, &[(0, 2)]), (1, &[(1, 1)])]);
let f2 = Polynomial::from_coeffs_int(&[(1, &[(0, 1)]), (1, &[(1, 2)])]);

let jacobian = Polynomial::jacobian(&[f1, f2]);
// J = [[2x, 1 ],
//      [1,  2y]]
assert_eq!(jacobian.len(), 2); // 2 functions

pub fn integrate(&self, var: Var) -> Polynomial

Compute the indefinite integral (antiderivative) of the polynomial with respect to a variable.

For a polynomial p(x), returns ∫p(x)dx. The constant of integration is implicitly zero.

Arguments

• var - The variable to integrate with respect to

Examples

use oxiz_math::polynomial::Polynomial;
use num_bigint::BigInt;
use num_rational::BigRational;

// f(x) = 3x^2
let f = Polynomial::from_coeffs_int(&[(3, &[(0, 2)])]);

// ∫f(x)dx = x^3
let integral = f.integrate(0);

// Verify: derivative of integral should be original
let derivative = integral.derivative(0);
assert_eq!(derivative, f);

pub fn definite_integral( &self, var: Var, lower: &BigRational, upper: &BigRational, ) -> Option<BigRational>

Compute the definite integral of a univariate polynomial over an interval (a, b).

For a univariate polynomial p(x), returns ∫ₐᵇ p(x)dx = F(b) - F(a) where F is the antiderivative of p.

Arguments

• var - The variable to integrate with respect to
• lower - Lower bound of integration
• upper - Upper bound of integration

Returns

The definite integral value, or None if the polynomial is not univariate in the given variable

Examples

use oxiz_math::polynomial::Polynomial;
use num_bigint::BigInt;
use num_rational::BigRational;

// ∫[0,2] x^2 dx = [x^3/3] from 0 to 2 = 8/3
let f = Polynomial::from_coeffs_int(&[(1, &[(0, 2)])]);
let result = f.definite_integral(0, &BigRational::from_integer(BigInt::from(0)),
                                    &BigRational::from_integer(BigInt::from(2)));
assert_eq!(result, Some(BigRational::new(BigInt::from(8), BigInt::from(3))));

pub fn find_critical_points(&self, var: Var) -> Vec<(BigRational, BigRational)>

Find critical points of a univariate polynomial by solving f’(x) = 0.

Critical points are values where the derivative equals zero, which correspond to local maxima, minima, or saddle points.

Arguments

• var - The variable to find critical points for

Returns

A vector of isolating intervals containing the critical points. Each interval contains exactly one root of the derivative.

Examples

use oxiz_math::polynomial::Polynomial;

// f(x) = x^3 - 3x = x(x^2 - 3)
// f'(x) = 3x^2 - 3 = 3(x^2 - 1) = 3(x-1)(x+1)
// Critical points at x = -1 and x = 1
let f = Polynomial::from_coeffs_int(&[
    (1, &[(0, 3)]),  // x^3
    (-3, &[(0, 1)]), // -3x
]);

let critical_points = f.find_critical_points(0);
assert_eq!(critical_points.len(), 2); // Two critical points

pub fn trapezoidal_rule( &self, var: Var, lower: &BigRational, upper: &BigRational, n: u32, ) -> BigRational

Numerically integrate using the trapezoidal rule.

Approximates ∫ₐᵇ f(x)dx using the trapezoidal rule with n subintervals. The trapezoidal rule approximates the integral by summing the areas of trapezoids.

Arguments

• var - The variable to integrate with respect to
• lower - Lower bound of integration
• upper - Upper bound of integration
• n - Number of subintervals (more subintervals = higher accuracy)

Examples

use oxiz_math::polynomial::Polynomial;
use num_bigint::BigInt;
use num_rational::BigRational;

// ∫[0,1] x^2 dx = 1/3
let f = Polynomial::from_coeffs_int(&[(1, &[(0, 2)])]);
let approx = f.trapezoidal_rule(0,
    &BigRational::from_integer(BigInt::from(0)),
    &BigRational::from_integer(BigInt::from(1)),
    100);

let exact = BigRational::new(BigInt::from(1), BigInt::from(3));
// With 100 intervals, approximation should be very close

pub fn simpsons_rule( &self, var: Var, lower: &BigRational, upper: &BigRational, n: u32, ) -> BigRational

Numerically integrate using Simpson’s rule.

Approximates ∫ₐᵇ f(x)dx using Simpson’s rule with n subintervals (n must be even). Simpson’s rule uses parabolic approximation and is generally more accurate than the trapezoidal rule for smooth functions.

Arguments

• var - The variable to integrate with respect to
• lower - Lower bound of integration
• upper - Upper bound of integration
• n - Number of subintervals (must be even for Simpson’s rule)

Examples

use oxiz_math::polynomial::Polynomial;
use num_bigint::BigInt;
use num_rational::BigRational;

// ∫[0,1] x^2 dx = 1/3
let f = Polynomial::from_coeffs_int(&[(1, &[(0, 2)])]);
let approx = f.simpsons_rule(0,
    &BigRational::from_integer(BigInt::from(0)),
    &BigRational::from_integer(BigInt::from(1)),
    100);

let exact = BigRational::new(BigInt::from(1), BigInt::from(3));
// Simpson's rule should give very accurate results for polynomials

pub fn eval_at(&self, var: Var, value: &BigRational) -> Polynomial

Evaluate the polynomial at a point (substituting a value for a variable).

pub fn eval(&self, assignment: &FxHashMap<Var, BigRational>) -> BigRational

Evaluate the polynomial completely (all variables assigned).

(Numeric substitution of assignment into the polynomial’s terms — no code/expression execution of any kind is involved.)

Panics

Panics if assignment does not cover every variable occurring in self. Callers that cannot guarantee a complete assignment ahead of time (e.g. because it was built incrementally, or comes from untrusted/partial input) should use Self::try_eval instead, which reports the same condition as None rather than aborting.

pub fn try_eval( &self, assignment: &FxHashMap<Var, BigRational>, ) -> Option<BigRational>

Evaluate the polynomial completely (all variables assigned), returning None instead of panicking if assignment is missing a variable that occurs in self.

pub fn eval_horner(&self, var: Var, value: &BigRational) -> BigRational

Evaluate a univariate polynomial using Horner’s method.

Horner’s method evaluates a polynomial p(x) = a_n x^n + … + a_1 x + a_0 as (((…((a_n) x + a_{n-1}) x + …) x + a_1) x + a_0), which requires only n multiplications instead of n(n+1)/2 for the naive method.

This method is more efficient than eval() for univariate polynomials.

Arguments

• var - The variable to evaluate
• value - The value to substitute for the variable

Returns

The evaluated value

Panics

Panics if the polynomial is genuinely multivariate in var — i.e. it still contains a variable other than var after substituting value for var, so it cannot be reduced to a scalar. Callers that cannot guarantee univariate-in-var input (e.g. because it comes from untrusted or incrementally built terms) should use Self::try_eval_horner instead, which reports the same condition as None rather than aborting.

pub fn try_eval_horner( &self, var: Var, value: &BigRational, ) -> Option<BigRational>

Evaluate a univariate polynomial using Horner’s method, returning None instead of panicking when the polynomial is genuinely multivariate in var (i.e. it still contains another variable after substituting value for var, so no scalar result exists).

This is the non-panicking counterpart of Self::eval_horner; see its documentation for the numeric contract.

pub fn substitute(&self, var: Var, replacement: &Polynomial) -> Polynomial

Substitute a polynomial for a variable.

pub fn integer_content(&self) -> BigInt

Integer content: GCD of all coefficients (as integers). Assumes all coefficients are integers.

pub fn primitive(&self) -> Polynomial

Make the polynomial primitive (divide by integer content).

pub fn gcd_univariate(&self, other: &Polynomial) -> Polynomial

GCD of two polynomials (univariate, using Euclidean algorithm).

pub fn pseudo_remainder(&self, divisor: &Polynomial, var: Var) -> Polynomial

Pseudo-remainder for univariate polynomials. Returns r such that lc(b)^d * a = q * b + r for some q, where d = max(deg(a) - deg(b) + 1, 0).

Division by the zero polynomial is mathematically undefined; rather than panicking on a crafted degenerate divisor (a soundness/DoS concern for any caller reachable from untrusted input), this treats it as “no reduction possible” and returns self unchanged, mirroring Self::exact_remainder_univariate’s convention for the same case.

pub fn pseudo_div_univariate( &self, divisor: &Polynomial, ) -> (Polynomial, Polynomial)

Pseudo-division for univariate polynomials. Returns (quotient, remainder) such that lc(b)^d * a = q * b + r where d = deg(a) - deg(b) + 1.

Division by the zero polynomial is mathematically undefined; rather than panicking on a crafted degenerate divisor, this returns (0, self) – “no division performed” – matching the convention used for the deg(a) < deg(b) case below.

pub fn subresultant_prs(&self, other: &Polynomial, var: Var) -> Vec<Polynomial>

Compute the subresultant polynomial remainder sequence (PRS).

The subresultant PRS is a more efficient variant of the pseudo-remainder sequence that avoids coefficient explosion through careful normalization. It’s useful for GCD computation and other polynomial algorithms.

Returns a sequence of polynomials [p0, p1, …, pk] where:

• p0 = self
• p1 = other
• Each subsequent polynomial is derived using the subresultant algorithm

Reference: “Algorithms for Computer Algebra” by Geddes, Czapor, Labahn

pub fn resultant(&self, other: &Polynomial, var: Var) -> Polynomial

Resultant of two polynomials with respect to a variable.

The resultant Res_var(p, q) eliminates var from {p = 0, q = 0}, returning a polynomial in the remaining variables that vanishes exactly when p and q share a common root (over an algebraic closure) or their leading coefficients in var both vanish.

• Univariate inputs (both operands contain only var): computed via the exact subresultant PRS with rational scalar normalization.
• Multivariate inputs (either operand contains another variable): computed exactly as the determinant of the (deg_p + deg_q)-square Sylvester matrix over the polynomial coefficient ring, evaluated division-free via the Samuelson–Berkowitz algorithm (see sylvester_resultant). This replaces the former primitive()-based approximation, which could return a mathematically wrong (non-zero) value for genuinely multivariate operands.

The sign/value convention is the classical Res(p, q) = lc(p)^{deg q} · Π_{p(α)=0} q(α); in particular Res(x - a, x - b) = a - b.

pub fn discriminant(&self, var: Var) -> Polynomial

Discriminant of a polynomial with respect to a variable. discriminant(p) = resultant(p, dp/dx) / lc(p)

pub fn is_definitely_positive(&self) -> bool

Check if the polynomial has a positive sign for all variable assignments. This is an incomplete check that only handles obvious cases.

pub fn is_definitely_negative(&self) -> bool

Check if the polynomial has a negative sign for all variable assignments. This is an incomplete check.

pub fn make_monic(&self) -> Polynomial

Make the polynomial monic (leading coefficient = 1).

pub fn square_free(&self) -> Polynomial

Compute the square-free part of a polynomial (removes repeated factors).

pub fn sturm_sequence(&self, var: Var) -> Vec<Polynomial>

Compute the Sturm sequence for a univariate polynomial. The Sturm sequence is used for counting real roots in an interval.

The chain is built as p_{i+1} = -rem(p_{i-1}, p_i) using the exact rational remainder. The pseudo-remainder must not be used here: it scales by lc(p_i)^k, whose sign flips when the leading coefficient is negative, which would break the Sturm sign invariant and yield wrong root counts (e.g. -x^2 + 1 on (-2, 2) would report 0 roots instead of 2).

pub fn count_roots_in_interval( &self, var: Var, a: &BigRational, b: &BigRational, ) -> usize

Count the number of real roots in an interval using Sturm’s theorem. Returns the number of distinct real roots in (a, b).

pub fn cauchy_bound(&self, var: Var) -> BigRational

Compute Cauchy’s root bound for a univariate polynomial. Returns B such that all roots have absolute value <= B.

Cauchy bound: 1 + max(|a_i| / |a_n|) for i < n

This is a simple, conservative bound that’s fast to compute.

pub fn fujiwara_bound(&self, var: Var) -> BigRational

Compute Fujiwara’s root bound for a univariate polynomial. Returns B such that all roots have absolute value <= B.

Fujiwara bound: 2 * max(|a_i/a_n|^(1/(n-i))) for i < n

This bound is generally tighter than Cauchy’s bound but more expensive to compute.

Reference: Fujiwara, “Über die obere Schranke des absoluten Betrages der Wurzeln einer algebraischen Gleichung” (1916)

pub fn lagrange_positive_bound(&self, var: Var) -> BigRational

Compute Lagrange’s bound for positive roots of a univariate polynomial. Returns B such that all positive roots are <= B.

This can provide a tighter bound than general root bounds when analyzing only positive roots, useful for root isolation optimization.

pub fn isolate_roots(&self, var: Var) -> Vec<(BigRational, BigRational)>

Isolate real roots of a univariate polynomial. Returns a list of intervals, each containing exactly one root. Uses Descartes’ rule of signs to optimize the search.

pub fn newton_raphson( &self, var: Var, initial: BigRational, lower: BigRational, upper: BigRational, max_iterations: usize, tolerance: &BigRational, ) -> Option<BigRational>

Refine a root approximation using Newton-Raphson iteration.

Given an initial approximation and bounds, refines the root using the formula: x_{n+1} = x_n - f(x_n) / f’(x_n)

Arguments

• var - The variable to solve for
• initial - Initial approximation of the root
• lower - Lower bound for the root (for verification)
• upper - Upper bound for the root (for verification)
• max_iterations - Maximum number of iterations
• tolerance - Convergence tolerance (stop when |f(x)| < tolerance)

Returns

The refined root approximation, or None if the method fails to converge or if the derivative is zero.

Examples

use num_bigint::BigInt;
use num_rational::BigRational;
use oxiz_math::polynomial::Polynomial;

// Solve x^2 - 2 = 0 (find sqrt(2) ≈ 1.414...)
let p = Polynomial::from_coeffs_int(&[
    (1, &[(0, 2)]),  // x^2
    (-2, &[]),       // -2
]);

let initial = BigRational::new(BigInt::from(3), BigInt::from(2)); // 1.5
let lower = BigRational::from_integer(BigInt::from(1));
let upper = BigRational::from_integer(BigInt::from(2));

let root = p.newton_raphson(0, initial, lower, upper, 10, &BigRational::new(BigInt::from(1), BigInt::from(1000000)));
assert!(root.is_some());

pub fn as_dense_i64(&self, var: Var) -> Option<Vec<i64>>

Extract the dense integer coefficient vector for a univariate polynomial.

Returns Some(coeffs) where coeffs[k] is the coefficient of var^k, if and only if every coefficient is an integer that fits in i64. Returns None for multivariate, non-integer, or out-of-range coefficients.

pub fn add_fast_i64(&self, other: &Polynomial, var: Var) -> Polynomial

Add two univariate polynomials using the SIMD-accelerated i64 fast path.

Falls back to the generic Polynomial::add when coefficients do not fit in i64 or the polynomials are not univariate in the same variable.

pub fn scale_fast_i64(&self, scalar: i64, var: Var) -> Polynomial

Scale a univariate polynomial by an integer scalar using the SIMD fast path.

Falls back to generic scalar-multiply when coefficients do not fit in i64.

pub fn dot_product_i64(&self, other: &Polynomial, var: Var) -> Option<i64>

Compute the i64 dot product between two same-degree dense coefficient arrays.

Returns Some(value) when both polynomials are expressible as dense univariate i64 coefficient arrays of equal length, None otherwise.

Trait Implementations

impl Add for Polynomial

type Output = Polynomial

The resulting type after applying the + operator.

fn add(self, other: Polynomial) -> Polynomial

Performs the + operation.

impl Add<&Polynomial> for &Polynomial

type Output = Polynomial

The resulting type after applying the + operator.

fn add(self, other: &Polynomial) -> Polynomial

Performs the + operation.

impl Clone for Polynomial

fn clone(&self) -> Polynomial

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for Polynomial

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl Display for Polynomial

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl Eq for Polynomial

impl Mul for Polynomial

type Output = Polynomial

The resulting type after applying the * operator.

fn mul(self, other: Polynomial) -> Polynomial

Performs the * operation.

impl Mul<&Polynomial> for &Polynomial

type Output = Polynomial

The resulting type after applying the * operator.

fn mul(self, other: &Polynomial) -> Polynomial

Performs the * operation.

impl Neg for Polynomial

type Output = Polynomial

The resulting type after applying the - operator.

fn neg(self) -> Polynomial

Performs the unary - operation.

impl Neg for &Polynomial

type Output = Polynomial

The resulting type after applying the - operator.

fn neg(self) -> Polynomial

Performs the unary - operation.

impl PartialEq for Polynomial

fn eq(&self, other: &Self) -> bool

Equality operator ==.

fn ne(&self, other: &Rhs) -> bool

Inequality operator !=.

impl Sub for Polynomial

type Output = Polynomial

The resulting type after applying the - operator.

fn sub(self, other: Polynomial) -> Polynomial

Performs the - operation.

impl Sub<&Polynomial> for &Polynomial

type Output = Polynomial

The resulting type after applying the - operator.

fn sub(self, other: &Polynomial) -> Polynomial

Performs the - operation.