Sturm Sequences

Overview

The Sturm sequence module implements sophisticated algorithms for:

Counting real roots in intervals
Isolating real roots
Finding complex roots in rectangular regions
Root refinement using Newton’s method
Computing Sturm sequences for polynomials

These methods are particularly effective for polynomials with symbolic or exact coefficients.

Complex Sturm Sequences

csturm_seq

modpoly&

First polynomial (typically P)

modpoly&

Second polynomial (typically P’)

listquo

vecteur&

Output: list of quotients Q_k

coeffP

vecteur&

Output: coefficients coeffP_k

coeffR

vecteur&

Output: coefficients coeffR satisfying the recurrence relation

contextptr

GIAC_CONTEXT

Context pointer

return

gen

GCD of r0 and r1 (without content)

gen csturm_seq(modpoly & r0, modpoly & r1, vecteur & listquo, 
               vecteur & coeffP, vecteur & coeffR, GIAC_CONTEXT);

Computes the Sturm sequence of two polynomials r0 and r1. Returns the GCD and computes quotients and coefficients such that:

coeffR * r_(k+2) = Q_k * r_(k+1) - coeffP_k * r_k

This forms the basis for counting sign changes and locating roots.

csturm_square

gen

Polynomial (as symbolic expression)

gen

Lower-left corner of rectangle (complex number)

gen

Upper-right corner of rectangle (complex number)

pgcd

gen&

Output: GCD if a factor is found

contextptr

GIAC_CONTEXT

Context pointer

return

int

Number of complex roots in rectangle, or -1 if factor found

int csturm_square(const gen & p, const gen & a, const gen & b, 
                  gen& pgcd, GIAC_CONTEXT);

Counts the number of complex roots of polynomial p inside the rectangle defined by corners a and b. Returns -1 if a factor is discovered during computation.

// Count roots in rectangle [0,2] × [0,2i]
gen a = 0;
gen b = gen(2, 2); // 2+2i
gen pgcd;
int nroots = csturm_square(x*x*x - 1, a, b, pgcd, contextptr);

Complex Root Finding

complex_roots

modpoly

Polynomial with numeric coefficients in Q[i]

gen

Lower bound for real part

gen

Lower bound for imaginary part

gen

Upper bound for real part

gen

Upper bound for imaginary part

roots

vecteur&

Output: vector of complex roots

pgcd

gen&

Output: factor of P if found

eps

double

Target precision

return

bool

False if a factor was found, true otherwise

bool complex_roots(const modpoly & P, const gen & a0, const gen & b0, 
                   const gen & a1, const gen & b1, vecteur & roots, 
                   gen & pgcd, double eps);

Finds all complex roots of P in the rectangle [a0,a1] × [b0,b1]i using complex Sturm sequences. Returns false if a factor is discovered.

complex_roots (high-level)

modpoly

Polynomial with numeric coefficients

gen

Lower bound for real part

gen

Lower bound for imaginary part

gen

Upper bound for real part

gen

Upper bound for imaginary part

complexe

bool

Whether to find complex roots

eps

double

Target precision

use_proot

bool

Whether to use proot algorithm as fallback

return

vecteur

Vector of roots

vecteur complex_roots(const modpoly & P, const gen & a0, const gen & b0, 
                      const gen & a1, const gen & b1, bool complexe, 
                      double eps, bool use_proot);

High-level interface for finding roots at specified precision. Automatically handles rectangular search regions.

complexroot / realroot

gen

Polynomial or expression

complexe

bool

Whether to find complex roots (internal version)

contextptr

GIAC_CONTEXT

Context pointer

return

gen

Roots of the polynomial

gen complexroot(const gen & g, bool complexe, GIAC_CONTEXT);
gen _complexroot(const gen & g, GIAC_CONTEXT);
gen _realroot(const gen & g, GIAC_CONTEXT);

Finds complex or real roots of polynomials using advanced root isolation methods.

Real Root Isolation

vas

modpoly

Polynomial

gen

Left endpoint of interval

gen

Right endpoint of interval

eps

double

Target precision

vasres

vecteur&

Output: list of intervals or rational roots

with_mult

bool

Whether to include multiplicities

contextptr

GIAC_CONTEXT

Context pointer

return

bool

True if successful

bool vas(const modpoly & P, const gen & a, const gen & b, double eps, 
         vecteur & vasres, bool with_mult, GIAC_CONTEXT);

Isolates and finds real roots of P in interval [a,b] at precision eps. Returns intervals containing roots or exact rational roots.

The VAS (Vincent-Akritas-Strzeboński) algorithm is particularly efficient for polynomials with large degree or exact coefficients.

newton_improve

vecteur

Polynomial coefficients

vecteur

Derivative polynomial coefficients

Preal

bool

Whether P has real coefficients

vecteur&

Approximate roots to refine

vradius

vecteur&

Output: radius of convergence for each root

int

Index of root to refine

kmax

int

Maximum Newton iterations

int

Current precision in bits

epsn

int

Target precision in bits

epsg2surdeg2

gen

eps^2 / degree(P)^2 as a gen

epsg

gen

Target precision as gen

return

bool

True if refinement successful

bool newton_improve(const vecteur & P, const vecteur & P1, bool Preal, 
                    vecteur & v, vecteur & vradius, int i, int kmax, 
                    int n, int epsn, const gen & epsg2surdeg2, 
                    const gen & epsg);

Refines approximate root v[i] using Newton’s method to achieve target precision. Updates vradius[i] with certified error bound.

newton_complex_1root

modpoly

Polynomial

gen

Lower bound for real part

gen

Lower bound for imaginary part

gen

Upper bound for real part

gen

Upper bound for imaginary part

complexroots

vecteur&

Output: found root (appended to vector)

eps

double

Target precision

return

bool

True if a root was found

bool newton_complex_1root(const modpoly & P, const gen & a0, const gen & b0, 
                          const gen & a1, const gen & b1, 
                          vecteur & complexroots, double eps);

Finds one complex root inside the rectangle using Newton’s method starting from the center.

Polynomial Root Algorithms

aberth

vecteur

Polynomial coefficients

vecteur&

Input/output: approximate roots

rayon

vecteur&

Output: error radii for each root

int

Degree of polynomial

eps

double

Target precision

isolate

int

Isolation mode flag

do_exact

bool

Whether to compute exact roots when possible

contextptr

GIAC_CONTEXT

Context pointer

return

bool

True if successful

bool aberth(const vecteur & P0, vecteur & R, vecteur & rayon, int N, 
            double eps, int isolate, bool do_exact, GIAC_CONTEXT);

Aberth’s method for simultaneous approximation of all polynomial roots. Highly efficient for well-conditioned polynomials.

The Aberth method simultaneously refines all roots, taking into account their mutual interactions. This often converges faster than individual Newton iterations.

mps_solve

vecteur

Polynomial coefficients

vecteur&

Output: approximate roots

rayon

vecteur&

Output: error radii

eps

double

Target precision

isolate

int

Isolation mode

secular

bool

Whether to use secular equation algorithm

contextptr

GIAC_CONTEXT

Context pointer

return

int

Status code

int mps_solve(const vecteur & P, vecteur & R, vecteur & rayon, double eps, 
              int isolate, bool secular, GIAC_CONTEXT);

Multiprecision polynomial solver using MPS (Multiprecision Polynomial Solver) algorithms.

Rational Root Finding

crationalroot / rationalroot

polynome&

Polynomial

complexe

bool

Whether to find complex rational roots (Gaussian rationals)

return

vecteur

Vector of rational (or Gaussian rational) roots

vecteur crationalroot(polynome & p, bool complexe);
gen _crationalroot(const gen & g, GIAC_CONTEXT);
gen _rationalroot(const gen & g, GIAC_CONTEXT);

Finds all rational roots (or Gaussian rational roots for complex version) using the rational root theorem.

Coordinate Transformations

change_scale

modpoly&

Polynomial to transform (modified in place)

gen

Scale factor

longlong

Power of 2 scale (if not 1 shifted left 31, use l=2^lb)

void change_scale(modpoly & p, const gen & l, longlong lb=1<<31);
void back_change_scale(modpoly & p, const gen & l, longlong lb=1<<31);

Transforms P(x) → P(l·x) for change_scale, and the inverse for back_change_scale. Useful for normalizing root search regions.

linear_changevar

modpoly

Polynomial

gen

Linear coefficient

gen

Constant term

return

modpoly

Transformed polynomial

modpoly linear_changevar(const modpoly & p, const gen & a, const gen & b);
modpoly inv_linear_changevar(const modpoly & p, const gen & a, const gen & b);

Transforms P(x) → P(a·x + b) and the inverse transformation.

ab2a0b0a1b1

gen

First corner of rectangle (complex)

gen

Opposite corner of rectangle (complex)

gen&

Output: minimum real part

gen&

Output: minimum imaginary part

gen&

Output: maximum real part

gen&

Output: maximum imaginary part

contextptr

GIAC_CONTEXT

Context pointer

void ab2a0b0a1b1(const gen & a, const gen & b, gen & a0, gen & b0, 
                 gen & a1, gen & b1, GIAC_CONTEXT);

Extracts rectangle boundaries from two complex corner points.

Utility Functions

horner_minmax

vecteur

Polynomial coefficients

Preal

bool

Whether P has real coefficients

gen

Evaluation point

eps

gen

Perturbation radius

Prmin

gen&

Output: minimum value of P in [r-eps, r+eps]

Prmax

gen&

Output: maximum value of P in [r-eps, r+eps]

void horner_minmax(const vecteur & P, bool Preal, const gen & r, 
                   const gen & eps, gen & Prmin, gen & Prmax);

Evaluates P at r±eps using Horner’s method and returns interval bounds.

round2

gen&

Value to round (modified in place)

int

Number of bits precision

void round2(gen & x, int n);
void in_round2(gen & x, const gen & deuxn, int n);

Rounds x to n bits of precision. Useful for controlling numerical stability.

keep_in_rectangle

croots

vecteur

Complex roots to filter

gen

Minimum real part

gen

Minimum imaginary part

gen

Maximum real part

gen

Maximum imaginary part

embed

bool

Whether to embed in surrounding space

contextptr

GIAC_CONTEXT

Context pointer

return

vecteur

Filtered roots within rectangle

vecteur keep_in_rectangle(const vecteur & croots, const gen A0, const gen & B0,
                          const gen & A1, const gen & B1, bool embed, 
                          GIAC_CONTEXT);

Filters a list of roots to keep only those within the specified rectangle.

symb2poly_num

gen

Symbolic expression

contextptr

GIAC_CONTEXT

Context pointer

return

vecteur

Numeric polynomial coefficients

vecteur symb2poly_num(const gen & g, GIAC_CONTEXT);

Converts symbolic polynomial to numeric coefficient vector.

Algorithm Overview

Sturm Sequences

Classical Sturm sequences count real roots in an interval by tracking sign changes. Complex Sturm sequences extend this to rectangular regions in the complex plane.The algorithm computes a sequence P₀, P₁, P₂, … where P₀ is the polynomial, P₁ is its derivative, and subsequent terms follow the Euclidean algorithm. Sign changes at interval endpoints indicate root presence.

VAS Algorithm

The Vincent-Akritas-Strzeboński algorithm efficiently isolates real roots using continued fractions and Descartes’ rule of signs. It’s particularly effective for high-degree polynomials.The method recursively subdivides intervals and applies transformations to isolate roots to arbitrary precision.

Aberth's Method

Aberth’s method simultaneously approximates all roots of a polynomial by solving:

z_i^(k+1) = z_i^(k) - P(z_i^(k))/P'(z_i^(k)) / (1 - Σ_{j≠i} 1/(z_i^(k) - z_j^(k)))

This accounts for interactions between nearby roots and often converges cubically.

Newton Refinement

Once roots are isolated, Newton’s method provides rapid local convergence:

z^(k+1) = z^(k) - P(z^(k))/P'(z^(k))

With multiprecision arithmetic, this achieves arbitrary accuracy.

Equation Solving - General equation solving functions
Polynomial arithmetic - See polynomial module
Root finding - See numeric module

Core API

Algebra

Calculus

Linear Algebra

Solvers

Utilities

Sturm Sequences

Overview

Complex Sturm Sequences

csturm_seq

csturm_square

Complex Root Finding

complex_roots

complex_roots (high-level)

complexroot / realroot

Real Root Isolation

vas

Newton Refinement

newton_improve

newton_complex_1root

Polynomial Root Algorithms

aberth

mps_solve

Rational Root Finding

crationalroot / rationalroot

Coordinate Transformations

change_scale

linear_changevar

ab2a0b0a1b1

Utility Functions

horner_minmax

round2

keep_in_rectangle

symb2poly_num

Algorithm Overview

Build docs developers (and LLMs) love

Core API

Algebra

Calculus

Linear Algebra

Solvers

Utilities

​Overview

​Complex Sturm Sequences

​csturm_seq

​csturm_square

​Complex Root Finding

​complex_roots

​complex_roots (high-level)

​complexroot / realroot

​Real Root Isolation

​vas

​Newton Refinement

​newton_improve

​newton_complex_1root

​Polynomial Root Algorithms

​aberth

​mps_solve

​Rational Root Finding

​crationalroot / rationalroot

​Coordinate Transformations

​change_scale

​linear_changevar

​ab2a0b0a1b1

​Utility Functions

​horner_minmax

​round2

​keep_in_rectangle

​symb2poly_num

​Algorithm Overview

​Related Functions

Build docs developers (and LLMs) love

Overview

Complex Sturm Sequences

csturm_seq

csturm_square

Complex Root Finding

complex_roots

complex_roots (high-level)

complexroot / realroot

Real Root Isolation

vas

Newton Refinement

newton_improve

newton_complex_1root

Polynomial Root Algorithms

aberth

mps_solve

Rational Root Finding

crationalroot / rationalroot

Coordinate Transformations

change_scale

linear_changevar

ab2a0b0a1b1

Utility Functions

horner_minmax

round2

keep_in_rectangle

symb2poly_num

Algorithm Overview

Related Functions