Polynomial<T> — Polynomial Template

Overview

Polynomial<T> is a polynomial template class in calx. The coefficient type T can be double, Rational, Complex<double>, or any type satisfying the PolynomialCoeff concept.

Coefficients are stored in ascending order: coeffs_[i] is the coefficient of $x^i$.

• Any coefficient type — from double to Rational (exact) or Complex<T>
• Horner evaluation — $O(n)$ evaluation of $p(x)$
• Polynomial division — quotient and remainder (divmod)
• Euclidean GCD — greatest common divisor of polynomials
• Differentiation / Integration — derivative(), integral()
• Orthogonal polynomials — Chebyshev, Legendre, Hermite, Laguerre, Bessel, Jacobi, Gegenbauer
• C++23 concepts — type-constrained via PolynomialCoeff

Build

Required Header

#include <math/core/Polynomial.hpp>

Header-only. No library linkage required.

Constructors

Polynomial<double> p = {1.0, 2.0, 3.0};  // 1 + 2x + 3x^2
Polynomial<Rational> q = {Rational(1,2), Rational(3,4)};  // 1/2 + 3/4 x

Accessors

Arithmetic Operators

Scalar operations: p + c, c * p, p / c. Compound assignment: +=, -=, *=, /=.

Division / Modulo

auto [q, r] = p.divmod(d);  // p = d * q + r, deg(r) < deg(d)

Division Mode — Switching operator/ Behavior

The return type of Polynomial / Polynomial can be switched at compile time via a macro.

Method 1: Define in source code

// Define before including RationalFunction.hpp
#define CALX_POLY_DIV_RATIONAL
#include <math/core/RationalFunction.hpp>

Polynomial<double> p = {1, 0, 1};  // x^2 + 1
Polynomial<double> q = {1, 1};     // x + 1

auto rf = p / q;       // RationalFunction<double>: (x^2+1)/(x+1)
auto quot = p.quo(q);  // Polynomial<double>: x - 1 (when quotient is needed)

Method 2: Configure via CMake

# CMakeLists.txt
option(CALX_POLY_DIV_RATIONAL
       "Polynomial / Polynomial returns RationalFunction" OFF)

if(CALX_POLY_DIV_RATIONAL)
    add_compile_definitions(CALX_POLY_DIV_RATIONAL)
endif()

This setting switches operator/ behavior project-wide. quo() and divmod() always perform polynomial division regardless of mode.

Evaluation

The argument type may differ from T (e.g., evaluate a Polynomial<double> at Complex<double>).

Passing a polynomial as the argument performs polynomial composition $p(q(x))$, equivalent to compose(p, q).

Polynomial<double> p({1, 0, 1});  // 1 + x^2
Polynomial<double> q({0, 1, 1});  // x + x^2
auto r = p(q);  // p(q(x)) = 1 + (x + x^2)^2 = 1 + x^2 + 2x^3 + x^4

The composed polynomial has degree $\deg(p(q)) = \deg(p) \cdot \deg(q)$.

Calculus

Algebraic Operations

Factorization & Resultant

Exact factorization of integer polynomials, plus resultant and discriminant. Include math/core/Polynomial_factorization.hpp.

#include <math/core/Polynomial.hpp>
#include <math/core/Polynomial_factorization.hpp>
#include <math/core/FactoredInteger.hpp>

factorize() — Irreducible factorization

Factor an integer polynomial $f(x) \in \mathbb{Z}[x]$ into irreducible factors.

Coefficient type T must be int or int64_t.

Pipeline:

1. Content extraction — GCD of all coefficients is pulled out as a constant factor.
2. Rational root theorem — search $f(b/a) = 0$ over $b \mid a_0$, $a \mid a_n$, remove all linear factors $(ax - b)$.
3. Yun's square-free factorization — use $\gcd(f, f')$ to split $f = \prod f_i^i$ where each $f_i$ is square-free and pairwise coprime.
4. Kronecker's method — for each square-free factor, interpolate trial factors from divisors of values at $\deg f / 2 + 1$ points.

// x^4 - 1 = (x - 1)(x + 1)(x^2 + 1)
Polynomial<int> f = {-1, 0, 0, 0, 1};
auto factors = factorize(f);

for (size_t i = 0; i < factors.size(); ++i) {
    std::cout << factors[i] << "^" << factors.exponent(i) << std::endl;
}
// (-1 + x)^1
// (1 + x)^1
// (1 + x^2)^1

The above is the core version (header-only, no linking required). Kronecker's method cuts off at 100,000 combinations, practical up to $\deg f \lesssim 10$.

Fast version (calx_factorize lib)

Link calx_factorize to automatically switch to a fast pipeline:

1. $\mathrm{GF}(p)$ factorization — Cantor-Zassenhaus (DDF + EDF)
2. Hensel lifting — linear + multi-factor binary-tree Hensel lift to $\mathbb{Z}[x]$
3. Factor recombination — $r \leq 12$: Zassenhaus enumeration, $r > 12$: van Hoeij lattice (LLL)

Int delegates to int64_t when all coefficients fit in 64 bits. Rational clears denominators (LCM) and delegates to the Int path.

// CMake: target_link_libraries(your_target PRIVATE calx_factorize)

squareFreeFactorization() — Yun's algorithm

Yun's algorithm (1976), implemented over subresultant PRS for $\gcd(f, f')$. Separates repeated factors with exponents; does not guarantee irreducibility.

Integer polynomial GCD and pseudo-remainder

Pseudo-remainder vs. ordinary remainder

The ordinary remainder assumes the coefficient ring is a field, so it does not work over the integers:

Example: dividing $f = x^2$ by $g = 2x + 1$ over the integers gives $q = \dfrac{x}{2} - \dfrac{1}{4}$, introducing fractions in the quotient (so neither quotient nor remainder exists in $\mathbb{Z}[x]$).

The pseudo-remainder avoids fractions by pre-multiplying by a power of the leading coefficient of $g$:

Multiplying by $\mathrm{lc}(g)$ in advance keeps all coefficients in the ring when dividing by $g$.

For the same example: with $\mathrm{lc}(g)^{2-1+1} = 2^2 = 4$,

so the pseudo-remainder is $r = 1$ — purely integer.

In calx, pseudoRemainder() exists to support intPolyGcd() (integer polynomial GCD) — it is the building block of the subresultant PRS that computes GCDs while keeping all coefficients integer.

Resultant & discriminant

// Discriminant of f(x) = x^2 - 2 is 4·1·2 = 8
Polynomial<int> f = {-2, 0, 1};
std::cout << discriminant(f) << std::endl;  // 8

A non-zero discriminant certifies that $f$ is square-free (no repeated roots).

Orthogonal Polynomials

Classical orthogonal polynomials generated by recurrence relations.

auto T5 = chebyshevT<double>(5);
std::cout << T5 << std::endl;
// 5x - 20x^3 + 16x^5

Type Traits

Example

#include <math/core/Polynomial.hpp>
#include <iostream>
using namespace calx;

int main() {
    // Construct polynomials
    Polynomial<double> p = {1.0, -3.0, 2.0};  // 1 - 3x + 2x^2
    Polynomial<double> q = {1.0, 1.0};         // 1 + x

    // Arithmetic
    std::cout << "p = " << p << std::endl;
    std::cout << "q = " << q << std::endl;
    std::cout << "p * q = " << p * q << std::endl;

    // Division
    auto [quot, rem] = p.divmod(q);
    std::cout << "p / q = " << quot << ", remainder = " << rem << std::endl;

    // Evaluation (Horner)
    std::cout << "p(2.0) = " << p(2.0) << std::endl;  // 3.0

    // Calculus
    std::cout << "p' = " << p.derivative() << std::endl;  // -3 + 4x
    std::cout << "integral(p) = " << p.integral() << std::endl;

    // GCD
    Polynomial<double> a = {-1.0, 0.0, 1.0};  // x^2 - 1 = (x-1)(x+1)
    Polynomial<double> b = {-1.0, 1.0};        // x - 1
    std::cout << "gcd(x^2-1, x-1) = " << gcd(a, b) << std::endl;

    // Orthogonal polynomials
    auto L3 = legendre<double>(3);
    std::cout << "P_3(x) = " << L3 << std::endl;
}

Related Mathematical Background

The following articles explain the mathematical concepts underlying the Polynomial class.

• Polynomial Representation and Basic Operations — Dense/sparse representation, arithmetic
• Polynomial GCD — Euclidean and subresultant algorithms
• Factorization over Finite Fields — Berlekamp, Cantor-Zassenhaus
• Factorization over Integers — Hensel lifting, LLL lattice basis reduction