The Sturm step

The Sturm step of ratpoints, _ratpoints_compute_sturm of sturm.c, which restricts the search to the intervals where the polynomial f may be ≥ 0, computed in integers and on the doubles of the C code as sturm.c computes it. First the remainder chain of f and its derivative f', row by row and with its loops (Ratpoints.Sturm.chain): a remainder is zero exactly when f has a repeated factor; the chain is none then, and ratpoints refuses the form as not squarefree, as Ratpoints.findPoints does with the Sturm step on. Then the numbers of sign changes of the chain at -∞ and +∞ from the leading coefficients (Ratpoints.Sturm.countsAtInfinity), with the exit when f has no real root and is negative; the bisection of the real line from [-∞, +∞], which drops an interval without a root on which f is negative, keeps one without a root on which it is positive and one at the depth of the step, and splits the others at their midpoint (Ratpoints.Sturm.midpoint), where it evaluates the chain (Ratpoints.Sturm.evalSign, Ratpoints.Sturm.signsAt), the left half first (Ratpoints.Sturm.iterate); the kept intervals, with dyadic ends n / 2 ^ e converted to doubles and an interval that begins where the last one ends joined to it (Ratpoints.Sturm.addInterval); and their intersection with the search intervals of the caller (Ratpoints.Sturm.intersect). The whole step is Ratpoints.Sturm.computeSturm, which the search calls through Ratpoints.sturmRegion.

Main definitions

• Ratpoints.Sturm.reduce: one reduction of a remainder by the row before it, which cancels its leading coefficient.
• Ratpoints.Sturm.chain: the remainder chain of f and f', or none when a remainder is zero.
• Ratpoints.Sturm.evalSign, Ratpoints.Sturm.signsAt: the sign of a row at a dyadic number, and the signs of f and f' and the number of sign changes of the chain there.
• Ratpoints.Sturm.countsAtInfinity: the numbers of sign changes of the chain at ±∞.
• Ratpoints.Sturm.midpoint, Ratpoints.Sturm.addInterval, Ratpoints.Sturm.iterate: the bisection.
• Ratpoints.Sturm.meet, Ratpoints.Sturm.intersect: the intersection with the caller's intervals.
• Ratpoints.Sturm.computeSturm: the Sturm step.

def Ratpoints.Sturm.reduce (s t : Array Int) : Array Int

One reduction of t by s, for deg t ≥ deg s, each given by the array of its coefficients up to its degree, constant term first: u t - v X^(deg t - deg s) s, where u and v are the leading coefficients of s and t divided by their positive gcd and signed so that u > 0; the leading coefficient of the result is 0.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sturm.chain (c : List Int) (d : Nat) : Id (Option (Array (Array Int)))

The remainder chain of the polynomial f of degree d ≥ 1 with the coefficients c₀, …, c_d, the entries of c, and c_d ≠ 0, in rows of integers: s₀ = f, s₁ = f', and for k = 2, …, d the row s_k from the remainder t of s_(k-2) by s_(k-1). Starting from t = s_(k-2), while the degree of t is at least that of s_(k-1), t is reduced by s_(k-1) (Ratpoints.Sturm.reduce) and its zero leading coefficients are dropped; the chain is none when t becomes zero, which happens exactly when f has a repeated factor (Ratpoints.Sturm.run_chain_eq_none_iff). Otherwise s_k is -t divided by the gcd of its coefficients, whose loop stops when the gcd is 1, and the chain ends when s_k is constant. Each row is the array of its coefficients up to its degree, constant term first; for d = 1 the loop over k does not run, and the chain is f, f'.

C: sturm.c:99-163

Deviation: the rows as arrays of their coefficients up to their degree, for the table of (d + 1)^2 GMP integers and the array of degrees, with no mpz_init and mpz_clear (a chain has at most d + 1 rows, Ratpoints.Sturm.IsChain.size_le, and its row k + 1 at most d - k entries, Ratpoints.Sturm.IsChain.size_getElem!_add_one_add_le): an array has one entry more than the degree of its row, and a remainder that becomes zero is the empty array, where the degree in the C code reaches -1; a reduction step drops the leading coefficient, which it makes 0, where the C code lowers the degree; the rows are handed back, and none for the return of -1, where _ratpoints_compute_sturm keeps them in its table for the rest of the Sturm step (Ratpoints.Sturm.computeSturm), and the index of the last row, which the C code leaves in k and sets to d after the loop for d = 1, is the size of the array less one; a list of integers for the array of GMP integers of the coefficients, local values for the scratch integers work[0], work[1] and work[2], and Lean's / for mpz_fdiv_q, which agree on the divisions of the chain, all exact.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sturm.evalSign (row : Array Int) (num denexp : Int64) : Id Int

The sign -1, 0 or 1 of the value at num / 2 ^ denexp of the polynomial whose coefficients, constant term first, are the entries of row, as eval_sign computes it: by Horner's rule in integers, from the leading coefficient w = c_deg, for n = deg - 1, …, 0 in turn w = w num + c_n 2 ^ e with e = denexp, 2 denexp, …, so that w is 2 ^ (deg denexp) times the value, deg the size of row less one.

C: sturm.c:67-84

Deviation: the row as an array of integers, for the GMP integers cofs and the long degree, which is the size of the array less one; integers for the GMP integers work[0] and work[1], the multiplication by a power of 2 for mpz_mul_2exp, and Int.sign for mpz_sgn; Int64 for the longs num, denexp and e; the steps counted up from 0 to deg - 1, for the index n that counts down; no scratch integers of the arguments.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sturm.countsAtInfinity (rows : Array (Array Int)) (k : Nat) : Id (Nat × Nat)

The numbers of sign changes of the rows 0, …, k of the chain rows at +∞ and at -∞, the pair (count1, count2), from the signs of their leading coefficients: two neighbouring rows count at +∞ when their leading coefficients have different signs, and at -∞ when they do after the sign of each row of odd degree is flipped.

C: sturm.c:165-175

Deviation: natural numbers for the long counts, which stay at most k; the rows as arrays, their degrees the sizes less one, for the arrays sturm and sturm_degs, and Int.sign for mpz_sgn; the counts returned, for the locals of _ratpoints_compute_sturm.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sturm.longMax : Int64

LONG_MAX = 2 ^ 63 - 1, as _ratpoints_compute_sturm computes it, the unsigned long -1 shifted right by one bit: the numerator of the end at +∞ of the bisection, and its negation that of the end at -∞; converted to doubles, they become ±2 ^ 63, beyond every x = a / b with |a| ≤ LONG_MAX.

C: sturm.c:191-192

Deviation: one constant for the two locals max and min = -max.

Equations

• Ratpoints.Sturm.longMax = (~~~0 >>> 1).toInt64

def Ratpoints.Sturm.midpoint (nl nr del der : Int64) : Int64 × Int64

The point nm / 2 ^ dem at which the bisection splits the interval from nl / 2 ^ del to nr / 2 ^ der: 0 for [-∞, +∞]; for [-∞, nr] the point -1 if nr = 0 and else 2 nr, and for [nl, +∞] the point 1 if nl = 0 and else 2 nl, both with dem = 0 (the ends at infinity are ∓Ratpoints.Sturm.longMax, with the exponent 0); for finite ends with del = der, which are then both 0, the half of nl + nr, with dem = 0 when the sum is even and dem = 1 otherwise; and else the sum of the ends over the larger denominator, divided by 2, so with the larger exponent plus 1.

C: sturm.c:234-251

Deviation: the point returned as a pair, for the locals nm and dem of iterate; Int64 for the longs, >>> on Int64 the arithmetic shift of the C code; the arithmetic of Int64 wrapping where that of the longs overflows, which it does not up to the depth 55 at which the search clamps the step, where every finite end is n / 2 ^ e with |n| ≤ 2 ^ 53.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sturm.signsAt (rows : Array (Array Int)) (k : Nat) (nm dem : Int64) : Id (Int × Int × Nat)

The signs s0 and s1 of f and f' at nm / 2 ^ dem (Ratpoints.Sturm.evalSign) and the number cmid of sign changes of the rows 0, …, k of the chain rows there, zeros skipped: the first two rows count when their signs are opposite; then the sign s of the last row that is not zero, s0 when s1 = 0, and each further row counts when its sign is -s, and replaces s when it is not zero.

C: sturm.c:252-260

Deviation: a function of its own for the evaluation of the chain at the midpoint in iterate, the signs and the count returned, for its locals s0, s1 and cmid; integers for the long signs and a natural number for the long count; the rows as arrays, for the arrays sturm and sturm_degs.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sturm.addInterval (iv : Array (Double × Double)) (nl nr del der : Int64) : Array (Double × Double)

The interval from nl / 2 ^ del to nr / 2 ^ der added to the intervals iv that the bisection keeps: its ends converted to doubles, ldexp((double)n, -e) (Ratpoints.Double.ofInt64, Ratpoints.Double.ldexp), which is exact up to the depth 55 and gives ±2 ^ 63 for the ends at infinity; the interval appended, or, when the last interval kept ends at its lower end (== of the doubles, Ratpoints.Double.beq), that interval extended to its upper end.

C: sturm.c:209-227

Deviation: a function of its own for the branch of iterate that keeps an interval; the intervals kept as an array, passed and returned, for the array ivlocal on the stack and the pointer iptr of _ratpoints_compute_sturm, which iterate reaches as a nested function of gcc, with no bound 1 + ⌊degree / 2⌋ on their number, the size of ivlocal, which the C code relies on and which holds assuming Sturm's theorem (Ratpoints.Sturm.size_run_iterate_le_of_chain_eq_some); the doubles as their bits (Ratpoints.Double); Int64 for the longs, converted to integers for the int casts of the exponents.

Equations

• One or more equations did not get rendered due to their size.

@[irreducible]

def Ratpoints.Sturm.iterate (rows : Array (Array Int)) (k iter : Nat) (nl nr del der : Int64) (cleft cright : Nat) (sl sr : Int) (depth : Nat) (iv : Array (Double × Double)) : Id (Array (Double × Double))

The bisection of _ratpoints_compute_sturm, the nested function iterate, on the interval from nl / 2 ^ del to nr / 2 ^ der, an end at infinity ∓Ratpoints.Sturm.longMax with the exponent 0, with the numbers cleft and cright of sign changes of the chain rows (index of the last row k) at its ends, as the step hands them on, and the signs sl and sr of f just inside its ends, at the depth depth, adding the intervals it keeps to iv: an interval with cleft = cright and sl < 0, without a root of f and with f < 0 on it, is dropped; one with cleft = cright and sl > 0, or at the depth iter, is kept (Ratpoints.Sturm.addInterval); and any other is split at its midpoint m (Ratpoints.Sturm.midpoint), where the chain is evaluated (Ratpoints.Sturm.signsAt), and its left half bisected before its right half, the count at m and the sign of f at m handed to both halves, and when m is a root of f, the count plus 1 and -f'(m) to the left half and the sign of f'(m) to the right half, the count and the sign just left and just right of the simple root m.

C: sturm.c:197-208, sturm.c:228-233, sturm.c:261-267

Deviation: a recursive function of its own, for the nested function of gcc, with the intervals kept passed and returned (Ratpoints.Sturm.addInterval), for the array ivlocal and the pointer iptr that it reaches, and the chain, the index of its last row and the depth of the step as arguments, for the arrays sturm and sturm_degs, k and iter that it reaches; natural numbers for the long counts, depths and the long iter, which is from 0 to 55 (the search clamps it, find_points.c:4682), integers for the long signs, and Int64 for the other longs; the recursion ends since the depth grows towards iter.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sturm.meet (a b : Double × Double) : Option (Double × Double) × Bool

One step of the intersection of the interval a that the bisection keeps with the interval b of the caller, both closed, as the loop of _ratpoints_compute_sturm makes it: the interval it notes, if any, and whether it goes on with the next interval kept (true) or the next interval of the caller (false), the one that ends first. If a begins first (a.low ≤ b.low): a is forgotten if it ends before b begins, and otherwise [b.low, a.up] is noted with the next a if a.up ≤ b.up, and b with the next b otherwise; symmetrically if b begins first. The ends are compared as doubles (Ratpoints.Double.le, Ratpoints.Double.lt) and copied, so that a -0.0 of the caller can stay.

C: sturm.c:289-322

Deviation: a function of its own for the body of the loop, the interval noted and the cursor to move returned, for the write into ivlist[n] and the increments of n, n1 and n2.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sturm.intersect (loc ivs : Array (Double × Double)) : Id (Array (Double × Double))

The intersection of the intervals loc that the bisection keeps with the intervals ivs of the caller, both increasing and disjoint, by the loop of _ratpoints_compute_sturm: from the first of each, while both lists have one left, a step (Ratpoints.Sturm.meet) notes their intersection, if they meet, and goes on with the next interval of the list whose interval ends first. The noted intervals may be degenerate, [a, a].

C: sturm.c:276, sturm.c:285-288, sturm.c:323, sturm.c:325

Deviation: the intervals of the caller read from the array it gives, where the C code reads them from a copy on the heap (lines 274-275, 286) and writes the intersection over them; the intersection returned as a new array, with no bound on its length, for the writes into ivlist, the caller's array, which must have room for 1 + ⌊degree / 2⌋ intervals more than it holds (lines 52-53), as the C code relies on and as suffices assuming Sturm's theorem (Ratpoints.Sturm.size_run_computeSturm_le), and for the number num_inter of the arguments that it sets; natural numbers for the long indices.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Sturm.computeSturm (c : List Int) (d iter : Nat) (ivlist : Array (Double × Double)) : Id (Option (Array (Double × Double)))

The Sturm step _ratpoints_compute_sturm on the polynomial f of degree d ≥ 1 with the coefficients c₀, …, c_d, the entries of c, and c_d ≠ 0, to the depth iter, for the search intervals ivlist of the caller, increasing and disjoint: none when the remainder chain is none (Ratpoints.Sturm.chain), for a polynomial with a repeated factor; the empty array when the chain has as many sign changes at -∞ as at +∞ (Ratpoints.Sturm.countsAtInfinity), so that f has no real root, and c₀ < 0, so that f < 0 everywhere; and otherwise the intervals that the bisection keeps from [-∞, +∞] (Ratpoints.Sturm.iterate, with the counts at ∓∞ at its ends and the sign of c_d, negated at -∞ for odd d, as the signs just inside them) intersected with ivlist (Ratpoints.Sturm.intersect), which is empty when they do not meet.

C: sturm.c:86-97, sturm.c:177-190, sturm.c:193-194, sturm.c:269-273, sturm.c:326-332

Deviation: the coefficients, the degree, the depth and the search intervals as arguments, for the fields cof, degree, sturm, domain and num_inter of the arguments; natural numbers for the longs degree and iter, which is from 0 to 55 (the search clamps it, find_points.c:4682); the result as an option of an array, for the long that the C code returns and the intervals it writes into domain with their number num_inter: none for -1, the empty array for 0, with num_inter = 0, and else the intervals, their number the return value; the table of the chain made by Ratpoints.Sturm.chain, whose rows the step reads, for the table on the stack with the degrees, which _ratpoints_compute_sturm fills itself, the index of the last row the size of the chain less one, for the k of line 163; the array of the intervals kept made by Ratpoints.Sturm.iterate, for the array ivlocal on the stack of 1 + ⌊degree / 2⌋ intervals, which suffice assuming Sturm's theorem (Ratpoints.Sturm.size_run_iterate_le_of_chain_eq_some); integers for the GMP integers and Int.sign for mpz_sgn; no mpz_init and mpz_clear, no copy of the caller's intervals on the heap and no failure for lack of memory (RATPOINTS_NO_MEMORY), which the model cannot meet.

Equations

• One or more equations did not get rendered due to their size.

Examples