November 2014 

NASA/TM–2014-218548 

A Formally-Verified Decision Procedure for 

Univariate Polynomial Computation Based on 

Sturm’s Theorem 

Anthony J. Narkawicz and César A. Muñoz 

Langley Research Center, Hampton, Virginia 
NASA STI Program . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. 

It collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI 

program provides access to the NASA Aeronautics 

and Space Database and its public interface, the 

NASA Technical Report Server, thus providing one 

of the largest collections of aeronautical and space 

science STI in the world. Results are published in 

both non-NASA channels and by NASA in the 

NASA STI Report Series, which includes the 

following report types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

Programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but 

having less stringent limitations on manuscript 

length and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. 

Collected papers from scientific and 

technical conferences, symposia, seminars, 

or other meetings sponsored or cosponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from 

NASA programs, projects, and missions, 

often concerned with subjects having 

substantial public interest. 

• TECHNICAL TRANSLATION.

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission. 

Specialized services also include organizing 

and publishing research results, distributing 

specialized research announcements and feeds, 

providing information desk and personal search 

support, and enabling data exchange services. 

For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page 

at http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443-757-5803 

• Phone the NASA STI Information Desk at 

443-757-5802 

• Write to: 

 STI Information Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320
National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

November 2014 

NASA/TM–2014-218548 

A Formally-Verified Decision Procedure for 

Univariate Polynomial Computation Based on 

Sturm’s Theorem 

Anthony J. Narkawicz and César A. Muñoz 

Langley Research Center, Hampton, Virginia 
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 

Acknowledgments 

Authors would like to thank Leonardo de Moura, Manuel Eberl, and anonymous referees 

for their comments on previous versions of this document. 

The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 

official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics 

and Space Administration. 
Abstract

Sturm’s Theorem is a well-known result in real algebraic geometry that provides a

function that computes the number of roots of a univariate polynomial in a semiopen interval. This paper presents a formalization of this theorem in the PVS

theorem prover, as well as a decision procedure that checks whether a polynomial is

always positive, nonnegative, nonzero, negative, or nonpositive on any input interval. The soundness and completeness of the decision procedure is proven in PVS.

The procedure and its correctness properties enable the implementation of a PVS

strategy for automatically proving existential and universal univariate polynomial

inequalities. Since the decision procedure is formally verified in PVS, the soundness of the strategy depends solely on the internal logic of PVS rather than on an

external oracle. The procedure itself uses a combination of Sturm’s Theorem, an

interval bisection procedure, and the fact that a polynomial with exactly one root

in a bounded interval is always nonnegative on that interval if and only if it is

nonnegative at both endpoints.

1
1 Introduction

Problems involving polynomial inequalities appear in applications such as air traffic

conflict detection and resolution algorithms [16], floating point analysis [5], and

uncertainty and reliability analysis of dynamic and control systems [3, 10]. Solving

these problems in a rigorous way is fundamental to the logical correctness of these

safety-critical systems. However, formal reasoning about polynomials and other

nonlinear functions in an interactive theorem prover is challenging.

Fortunately, significant advances have been made in this area in recent years.

In addition to related work described in later sections, the authors developed formalizations in PVS [24] of multivariate Bernstein polynomials [22] and a generic

branch and bound algorithm [23], both of which yield well-known numerical approximation methods. These PVS developments include formally-verified semi-decision

procedures for checking validity and satisfiability of nonlinear properties involving

variables ranging over real intervals. The procedures are integrated into the PVS

theorem prover as automated proof-producing strategies such as bernstein, which

uses Bernstein polynomials, and interval, which uses interval arithmetic. To automatically discharge a formula such as

x120 − 2x60 + 1.001 > 0, (1)

whenever x ∈ [0, 3], the user simply needs to invoke one of those strategies in PVS.

In particular, the user does not need the insight that the polynomial in Formula (1)

is equal to (x60 − 1)2 + 0.001. While these strategies are powerful, they inherit the

downsides of other numerical approximation methods. For instance, neither of these

strategies succeed in discharging Formula (1) when the variable x is unbounded even

though the inequality still holds for any real number x. Moreover, in many cases,

only approximations of the given formulas are used for reasoning. In general, Bernstein polynomials and other numerical tools like interval arithmetic can compute

a tight bound for the range of a polynomial but not the exact range. Thus, the

strategies bernstein and interval cannot prove that

x120 − 2x60 + 1 ≥ 0, (2)

whenever x ∈ [0, 3], because 0 is the precise minimum of the polynomial on this

interval.

This paper addresses shortcomings of numerical approximation methods for the

special case of univariate polynomials. In particular, this paper presents a decision

procedure, which is formally verified in PVS, that can be used to check satisfiability

and validity of single-variable strict and non-strict polynomial inequalities where

the variable is restricted to any interval, even an infinite or an open interval. For

example, it can prove Formula (2) for values of x in [0, 3], (0, 3), (−∞, 0], [3, ∞),

(−∞, ∞), etc.

The decision procedure presented in this paper is based on Sturm’s theorem.

This is a theorem from real algebraic geometry that can be used to explicitly compute

the number of roots of a univariate polynomial in a bounded, semi-open interval

(a, b]. In the case where the polynomial does not have a multiple root at either a
or b, Sturm’s Theorem provides a computable function σ such that the number of

roots is equal to σ(a) − σ(b). Moreover, the function σ can be evaluated at −∞

and ∞, making it possible to explicitly count the number of roots in the intervals

(−∞, b] and (a, ∞). Thus, Sturm’s Theorem can be used to prove that a univariate

polynomial is always positive (respectively, always negative) on an interval, provided

that the polynomial does not have a multiple root at either endpoint. This proof can

be done by checking that there are no roots in the interval and that the polynomial

is nonnegative (respectively, nonpositive) at the endpoints.

Sturm’s Theorem allows many useful results to be proved, but it does not directly

apply to problems where an endpoint is a multiple root or where the polynomial

bound is exact, e.g., (x − 1

3 )2 ≥ 0 for all x ∈ (0, ∞). This paper expands the

functionality of Sturm’s Theorem into a decision procedure that addresses both of

these cases and applies to any connected interval, open or closed. To deal with

endpoints that are multiple roots, the decision procedure finds a small positive

number  that is less than the distance between any two roots of the polynomial. It

then perturbs by  either endpoint that is a multiple root. The new interval contains

the same number of roots as the original. This allows counting roots in any connected

interval. Finally, to deal with the fact that Sturm’s Theorem can not directly help

prove the nonnegativity (respectively, nonpositivity) of a polynomial that has roots

in a given interval, a branching algorithm is defined that progressively subdivides

the original interval into smaller intervals until each smaller interval contains at

most one root of the polynomial. Then nonnegativity (respectively, nonpositivity)

is checked by whether or not the function is nonnegative (respectively, nonpositive)

at each endpoint of these subintervals. As noted above, the final result is a decision

procedure that can be used to prove bounds on any univariate polynomial on any

connected interval. This decision procedure, which combines Sturm’s Theorem and

a branching algorithm, appears to be new to the literature. Sturm’s Theorem is

used most often to count roots rather than to prove nonnegativity or positivity.

The decision procedure presented in this paper has been formally proved in PVS.

This paper presents not only the decision procedure itself, but also a proof-producing

strategy in PVS, called sturm, for discharging existentially and universally quantified univariate polynomial inequalities. It can handle any such problem, it always

terminates, and it is always correct. The soundness of the strategy only relies on the

PVS deduction engine. In particular, the strategy does not depend on any trusted

external oracle.

The formal development presented in this paper is electronically available in the

contribution Sturm of the NASA PVS Library.1 Examples are provided in the theory

Sturm@examples.pvs. All theorems presented in this paper are formally verified in

PVS. For readability, standard mathematical notation is used throughout this paper.

The reader is referred to the formal development for implementation details.

1http://shemesh.larc.nasa.gov/fm/ftp/larc/PVS-library.

2
2 Sturm’s Theorem

Sturm’s Theorem [28] is easily stated. Let p be a univariate polynomial. A function

σ is defined on the extended real line R ̄ = R ∪ {−∞, ∞} by setting σ(x) to be equal

to the number of sign changes in the sequence

p0(x), p1(x), p2(x), ... ,pm(x), (3)

where

p0(x) = p(x),

p1(x) = p

(x),

∀ j > 1 : ∃ c ∈ R+ : pj (x) = −c · rem(pj−2, pj−1)(x),

pm = 0, and

∃ y ∈ R : pm−1(y) = 0.

(4)

Here, rem(pj−2, pj−1) is the remainder after polynomial division of pj−2 by pj−1.

Such a sequence is called a Sturm sequence of polynomial p. When counting the

number of sign changes in the Sturm sequence, any zeros are ignored. For example, if

m = 7 and p0(x) = 4, p1(x) = −3, p2(x) = −5, p3(x)=0, p4(x) = 18, p5(x) = −4,

p6(x) = −1 and p7(x) = 0, there are sign changes between p0(x) and p1(x), between

p2(x) and p4(x), and between p4(x) and p5(x). In this case, the number of sign

changes in the sequence is given by σ(x) = 3.

Note that the evaluation of a polynomial at −∞ or ∞ is equal to −∞ or ∞

depending on the degree of the polynomial and its leading coefficient. Thus, it is

possible to count the number of sign changes in the sequence of polynomials even if

x is equal to −∞ or ∞. Also, it is worth noting that the introduction of the positive

real number c in Formula (4) is not included in most statements of Sturm’s Theorem

where c is set to 1. However, a specific sequence will be constructed in the algorithm

presented later that improves the computational efficiency of the standard method

and uses different values for c, an idea credited to de Moura and Passmore [6].

Sturm’s Theorem is stated as follows.

Theorem 1 (Sturm). For a, b ∈ R ̄ with a<b, if either p(a) = 0 or p

(a) = 0, and

if either p(b) = 0 or p

(b) = 0, then the number of roots of p in the interval (a, b] is

given by σ(a) − σ(b).

Numerous proofs of Sturm’s Theorem can be found on the internet. The authors

would be content to avoid redundancy and leave the proof out of this paper, except

for the fact that there is often a gulf between how a theorem is proved on paper

and how it is proved in an interactive theorem prover. The intended audience is

both mathematicians and users of interactive theorem provers. The proof below

breaks the main idea down into basic principles, and the authors hope that it can

be replicated easily in another interactive theorem prover. It is not a verbatim copy

of the PVS proof of the theorem, but it does use the same approach. This proof

is the same in essence as several others that are found online, most notably the

one found in Sottile’s course notes [27]. The main idea behind the PVS proof is

given here, which is inspired by the proof given in [27], although it differs in its key

3
approach in that it does not use a previous result about Sylvester sequences. The

proof that is given here is both more elementary and more detailed.

Proof of Theorem 1. First, note that if either a = −∞ or b = ∞, then there is

another interval (a∗, b∗] with real numbers as endpoints that is contained in (a, b]

with the same number of roots of p as the interval (a, b], and such that for every

j, pj (a∗) has the same sign as pj (a) and pj (b∗) has the same sign as pj (b). Thus,

σm(a∗) = σm(a) and σm(b∗) = σm(b). Hence, it is sufficient to prove the result for

the interval (a∗, b∗]. This shows that the proof can be reduced to the case where a

and b are real numbers (not infinite). Thus, it is henceforth assumed that a and b

are simply real numbers with a<b.

It is relatively easy to prove in an interactive theorem prover that there is some

sequence

a = a0 < a1 < ··· < ak = b

of real numbers such that each closed interval [ai, ai+1] contains at most one point

that is the root of any pj , and that this root is in the interior interval (ai, ai+1) unless

either the root is a and i = 0 or the root is b and i = k. If Sturm’s Theorem can

be proved on each semi-open interval (ai, ai+1], then the main result will trivially

follow. Thus, instead of proving Sturm’s Theorem when there are multiple roots, it

is sufficient to assume that (1) there is at most one point that is a root of any of the

pj ’s in the closed interval [a, b], and (2) if that root is equal to either a or b then it

is not a root of both p and p

.

If there are no roots of any pj in this closed interval then the intermediate value

theorem implies that each pj (a) has the same sign as pj (b), so σ(a)−σ(b) = 0. This

leaves only the case that there exists exactly one point that is a root of some pj in

[a, b] and where this root can only be either a or b, if it is not a root of both p and

p

. Denote this root r. It can be verified from the definition of the sequence {pj}

that the following properties hold:

1. For all j<m, there exists a nonnegative integer kj and another polynomial

gj such that p(x)=(x − r)kj · gj (x) for all x, where gj (r) = 0.

2. If k0 > 0, then k1 = k0 − 1.

3. If k0 > 1, then kj = k0 − 1 and pj (r) = 0 whenever 1 ≤ j<m.

4. If two consecutive elements of the sequence {pj} both have a root at r, then

every pj has a root at r, r = a, r = b, k0 > 1, and kj > 0 for all j.

5. If pm−1(r) = 0, then kj > 0 and pj (r) = 0 for all j.

6. If k0 ≤ 1, then pm−1 never changes sign on [a, b].

7. If pj (r) = 0, then pj (a) and pj (b) have the same nonzero sign.

8. If k0 ≤ 1, 1 ≤ j<m and pj (r) = 0, then pj−1(a) and pj+1(b) are nonzero and

have opposite signs, and pj−1(a) and pj+1(b) are nonzero and have opposite

signs.

4
9. If k0 > 1 and k0 is even, then p0(a) and p0(b) are nonzero and have the same

sign, and for 0 <j<m, pj (a) and pj (b) have opposite signs.

10. If k0 > 1 and k0 is odd, then p0(a) and p0(b) are nonzero and have opposite

signs, and for 0 <j<m, pj (a) and pj (b) are nonzero and have the same sign.

11. If r = a and r = b, then p0(a) and p1(a) are nonzero have opposite signs, and

p0(b) and p1(b) are nonzero and have the same sign.

12. If r = a and p0(a) = 0, then both p1(a) and p0(b) have the same sign as p1(b),

which is nonzero.

13. If r = b and p1(b) = 0, then p0(a) and p1(a) are nonzero and have opposite

signs, while p1(a) has the same sign as p1(b).

For the purpose of the proof, define σj (x) to be the number of sign changes in the

first j terms of the sequence in Formula (3), so that σ = σm. The proof is carried

out in each case by induction on the subscript j in σj .

First consider the case where r = a. Then k0 ≤ 1 by Property 4. It is proved

by induction on j that if pj (a) = 0, then pt(a) and pt(b) have the same sign and

that σt(a) − σt(b) = 0 for either t = j or t = j + 1. The result will follow because

pm−1(a) is nonzero by properties 5 and 4, and σm = σm−1. The base case is either

j = 0 or j = 1. If p0(a) = 0, then the base case is j = 1, and Property 12 implies

that p1(a) and p1(b) have the same sign and that σ1(a) = σ1(b) = 0. If p0(a) = 0,

then the base case is j = 0, and Property 7 implies that p0(a) and p0(b) have the

same sign and that σ0(a) = σ0(b) = 0. Now by induction, suppose that pj (r) = 0.

Then, by Property 7, pj (a) and pj (b) have the same nonzero sign. If pj−1(r) = 0,

then pj−1(a) and pj−1(b) have the same sign as well, so it follows immediately that

σj (a) − σj−1(a) = σj (b) − σj−1(b), and the result follows by induction. Conversely,

if pj−1(r) = 0, then Property 8 implies that pj−1(a) and pj+1(b) are nonzero and

have opposite signs, and pj−1(b) and pj+1(b) are nonzero and have opposite signs.

It follows immediately that σj+1(a) = σj−1(a) + 1 and σj+1(b) = σj−1(b) + 1.

Furthermore, Property 4 implies that pj+1(a) = 0. By Property 7, pj+1(a) and

pj+1(b) have the same sign. This completes the proof in the case where r = a.

The case where r = b follows directly from the cases where r = a and where

a<r<b. To see this, suppose that r = b and note that there is some  > 0

such that no pj has a root in the semi-open interval (b, b + ]. By applying the

result to the interval (a, b + ], which contains r in its interior, it can deduced that

σm(a)−σm(b+) = 1. Similarly, by applying the result to the interval (b, b+], whose

greatest lower bound is r, it can deduced that σm(b) − σm(b + ) = 0. Combining

these results yields σm(a)−σm(b) = 1. This reduces the argument to the case where

a<r<b, which is assumed for the remainder of the proof.

First suppose that k0 > 0, which by Property 3 implies that pj (r) = 0 for all j.

It is proved by induction on j that σj (a) − σj (b) = 1 for all j ≥ 1. For the base

case, it follows immediately from Property 11 that σ1(a) = 1 and σ1(b) = 0, and

therefore σ1(a) − σ1(b) = 1. Define γ to be 1 if k0 is even and −1 if k0 is odd. For

the inductive case, note that properties 9 and 10 imply that for all j ≥ 1, the sign

of σj (a) is equal to −γ multiplied by the sign of σj (b). From this, it can be seen

5
that pj+1(a) and pj (a) have opposite sign if and only if pj+1(b) and pj (b) do as well,

implying that σj+1(a)−σj (a) = σj+1(b)−σj (b), and the result follows. This reduces

the proof to the case where a<r<b and k0 ≤ 1.

It is now proved by induction that for all j ≥ 1, σj (a) − σj (b) = 1 whenever

pj (r) = 0. The result follows because Property 6 implies that pm−1(r) = 0 and σm =

σm−1. For the base case, it follows immediately from Property 11 that σ1(a)=1

and σ1(b) = 0, and therefore σ1(a)−σ1(b) = 1. For the inductive case, suppose that

j > 1 and that pj (r) = 0.

If pj−1(r) = 0, then Property 7, when applied to both j − 1 and j, gives that

pj−1(a) and pj−1(b) have the same nonzero sign, and that pj (a) and pj (b) have the

same nonzero sign. From this, it can be seen that pj−1(a) and pj (a) have opposite

sign if and only if pj−1(b) and pj (b) do as well, implying that σj (a) − σj−1(a) =

σj (b) − σj−1(b), and the result follows.

Finally, suppose that pj−1(r) = 0. Property 4 implies that pj−2(r) and pj (r)

are both nonzero. Property 8 implies that pj−2(a) and pj (a) are nonzero and have

opposite signs and that pj−2(b) and pj (b) are nonzero and have opposite signs.

Thus, σj (a) = σj−2(a) + 1 and σj (b) = σj−2(b) + 1. If j − 2 ≥ 1, then the inductive

hypothesis can be applied to deduce that σj (a) − σj (b) = 1. This completes the

proof except in the case where j = 2. However, then j −1 = 1, and so by hypothesis

p1(r) = 0. Property 1 would then imply that k1 > 0, and so by Property 2, k0 > 1,

a contradiction, since it is assumed in this part of the proof that k0 ≤ 1.

Sturm’s Theorem is less efficient than other methods at isolating roots. However, as this paper shows, Sturm’s Theorem works well for determining whether the

polynomial is always positive, negative, non-zero, etc., on a given interval.

3 Sturm Sequence of Integer Polynomials

An important element in the definition of a decision procedure for polynomial reasoning is a function that explicitly computes the remainder between two polynomials.

For the development presented in this paper, pseudo division is used rather than

standard polynomial division. This removes the need for dividing coefficients in the

algorithm, which means that if the coefficients of the original polynomials are integers, then the coefficients of their remainder after division will also be integers. In

fact, the pseudo division algorithm in PVS is defined only for integer polynomials.

The remainder after pseudo division of a polynomial g by a polynomial h is

not equal to the remainder after standard division, but is a power of the leading

coefficient of h multiplied by the standard remainder. This power can be either positive or negative. In the PVS development, the pseudo remainder is multiplied by a

constant so that this multiple is always negative. This ensures that the conditions

in Formula (4), defining the polynomial remainder sequence, all hold. Moreover,

pseudo division can be costly in terms of memory, because avoiding divisions can

cause integer coefficients to become quite large in the pseudo remainder. Thus, in

the PVS development, the pseudo remainder of two integer polynomials is also multiplied by the reciprocal of the greatest common division of all its coefficients, which

still ensures that the output is an integer polynomial. This multiplication helps to

6
mitigate coefficient growth, and since the greatest common divisor is positive, it does

not affect the conditions in Formula (4). In the PVS development, the function gcd

on polynomials is defined such that gcd(p) is equal to the greatest common divisor

of the coefficients of p whenever p is nonzero.

Formally, a univariate polynomial function p is a function from real numbers to

real numbers such that

p(x) ≡ n

i=0

cixi, (5)

where cn = 0. The natural number n is the degree of p and real numbers ci are the

coefficients of p. If ci is an integer for every i ≤ n then p will be called an integer

polynomial. If each such ci is a rational number, then p will be called a rational

polynomial. In order to specify functions and properties involving polynomials in

an interactive theorem prover, it is necessary to define a data structure to represent

them. The development presented in this paper uses a degree-dense representation

where a polynomial is a list of numerical coefficients of type T, i.e., the i-th element

of the list represents the coefficient of the xi monomial. The type T can be real,

for real polynomials, int, for integer polynomials, or rat, for rational polynomials.

Nothing in this paper fundamentally depends on this particular representation. For

readability, when there is no ambiguity, lower case letters f, g, . . . , p, q, . . . will be

used for both the mathematical concept of polynomial, e.g., Formula (5), and the

polynomial data structure used in the definition of PVS functions, e.g., Formula (6)

below. In the latter case, the expressions deg(p) and coeff(p, i) refer to the degree

of p and the i-th coefficient of the polynomial p, respectively. The same remark

applies to p

, etc. It denotes the mathematical definition of the derivative of p and

the data structure that represents this polynomial.

The adjusted pseudo remainder can be calculated for two integer polynomials g

and h as follows.

adjusted rem(g, h) ≡

let r = pseudo rem(g, h),

d = gcd(p) in

if coeff(h, deg(h)) > 0 ∨ mod(deg(g) − deg(h)+1, 2) = 0

then − r/d

else r/d endif.

(6)

Here, pseudo rem(g, h) is the remainder after pseudo division of g by h. Its definition

is common enough that it is left out here, although it can be found in the PVS

development. Also, note that coeff(h, deg(h)) is the leading coefficient of h, i.e.,

the nonzero coefficient with the highest degree. It can be verified that since g and

h both have integer coefficients, the polynomial adjusted rem(g, h) does as well.

The function comp rem seq, defined below, computes the remainder sequence

for any two integer polynomials g and h such that the degree of h is less than the

degree of g. It has a list  of integer polynomials as an input, which is also used as

7
an accumulator to store the sequence that is recursively computed by the function.

comp rem seq(g, h, ) ≡

if length()=0 then comp rem seq(g, h, cons(g, ))

elsif deg(head()) = 0 then 

elsif length()=1 then comp rem seq(g, h, cons(h, ))

else let p = adjusted rem(head(tail()), head()) in

comp rem seq(g, h, cons(p, ))

endif.

(7)

The function sturm seq, defined below, computes a Sturm sequence of an integer

polynomial p.

sturm seq(p) ≡ comp rem seq(p, p

, φ), (8)

where φ refers to the empty list. A curious reader may have noticed that sturm seq

computes a Sturm sequence in reverse order, i.e., if  = sturm seq(p), pm in Formula (3) corresponds to head(), pm−1 corresponds to head(tail()), and so on. It

can be easily checked that the sign changes of  and reverse() are equal.

4 Decision Procedure for Integer Polynomials

This section presents a decision procedure for computing the sign of an integer

polynomial p, i.e., its positivity, nonpositivity, negativity, nonnegativity, or nonzero

property, on a nonempty interval, which may or may not have infinite endpoints.

This decision procedures depends on a function that explicitly computes the number

of roots of p in that interval. One immediately obvious problem with using the

function σ from Sturm’s Theorem to define such function is that it will not work

when either the lower bound or the upper bound of the interval is a multiple root

of the polynomial, i.e., when the polynomial and its derivative both have roots at

that point.

The problem of multiple roots at the endpoints of an interval can be addressed

by perturbing such bounds outward by a small amount so that the new interval

contains exactly the same number of roots as the original but has endpoints that

are not multiple roots. To see how such a perturbation can be computed, let r be

a root of p, so p(r) = 0. A function can be defined with p and r as inputs that

computes the width of a small interval around r that contains no other roots of p.

The first step is to compute the degree of the first successive derivative of p that

does not vanish at r, which is accomplished with the following recursive function.

md(p, r) ≡ if deg(p)=0 ∨ p(r) = 0 then deg(p)

else md(p

, r) endif. (9)

The function md is well defined since the degree of the polynomial that is passed as

parameter strictly decreases at each recursive call.

By using Taylor’s Theorem, p is approximated in a small neighborhood of

r by p(n−md(p,r))(r) · (x − r)n−md(p,r)

. The error is bounded by a constant times

8
(x − r)n−md(p,r)+1, where p(n−md(p,r))(r) is the (n − md(p, r))-th derivative of p at

r. The interval around r containing no other roots is determined by computing a

neighborhood of r on which this derivative is always nonzero. Since it is always

nonzero, it can be proved by induction and the mean value theorem that all lesser

derivatives vanish only at r on this neighborhood. This neighborhood is computed

as follows. First, the following function is defined that takes as inputs polynomial

p, x ∈ R,, and  ∈ R > 0, and returns a δ ∈ R > 0 so that any input within δ of x

will have value less than  from x when p is applied to both points.

pcc(p, x, ) ≡

let n = deg(p),

c = n max

i=0 |coeff(p, i)| + 1 in

if n = 0 then 1/2

else min

⎛

⎝



c

⎛

⎝1 +n

i=1



i

j=1

 i

j − 1



|x|

j−1

⎞

⎠ ,

1

2

⎞

⎠

endif.

(10)

Using the function pcc, it is now possible to define a radius around the point r

in which the polynomial p has no other roots. This radius is computed with the

function root rad.

root rad(p, r) ≡ let n = deg(p) in

pcc(p(n−md(p,r)), r, |p(n−md(p,r))(r)|). (11)

Note that the function root rad can be used to compute the number of roots in

any interval, not simply an interval without a multiple root at either endpoint. To

see this by example, note that if (a, b] is an interval of real numbers such that b is

a multiple root of p but a is not, then the interval (a, b + root rad(p, b)/2] contains

the same number of roots as (a, b] but neither endpoint is a multiple root.

The next step in the development is to define a function called roots cl int

with two extended real numbers a and b as inputs, with a<b, along with an integer

polynomial p. The function returns the number of roots in the closed interval [a, b].

The other input to this function is a list  of polynomials, which in practice is set

to sturm seq(p).

roots cl int(p, a, b, ) ≡

let a∗ = if a = −∞ ∨ p(a) = 0 ∨ p

(a) = 0 then a

else a − root rad(p, a)/2 endif

b∗ = if b = ∞ ∨ p(b) = 0 ∨ p

(b) = 0 then b

else b + root rad(p, b)/2 endif

c = if a = −∞ ∧ p(a)=0 ∧ p

(a) = 0 then 1

else 0 endif

in sigma(, a∗) − sigma(, b∗) + c.

(12)

9
The definition of roots cl int uses the function sigma, which implements the

function σ from §2. This function takes as inputs a finite list  of polynomials and

an extended real number r, and returns the number of sign changes in that list when

each element is evaluated at r. The introduction of the number c in the definition of

roots cl int addresses the limitation in Sturm’s Theorem, which only gives a way

to count the number of roots in a half open interval that does not include its lower

bound. The term c adjusts this number based on whether the lower bound is equal

to the newly computed a∗ and is also a root of p but not of its derivative. Finally, it

should also be noted that the numbers a and b in the definition of roots cl int are

extended real numbers and therefore may be equal to −∞ or ∞. There is no built-in

support for extended reals in PVS. Hence, a simple data structure for representing

these numbers has been defined. Using this data structure, the evaluation of a

polynomial p at an extended real r ∈ {−∞, ∞} is defined as ±∞, where the sign

of the infinite number is given by sign(coeff(p, deg(p))), if deg(p) is even, and

sign(r) · sign(coeff(p, deg(p))), otherwise. Furthermore, it is assumed that −∞

is smaller than any other extended real number and ∞ is greater than any other

extended real number. No other operations are defined on infinite extended real

numbers.

The following theorem states the correctness of the function roots cl int. It

has been formally proved in PVS, where it is called roots closed int def.

Theorem 2. Let a, b ∈ R ̄, with a<b, p be an integer polynomial, and S = {r ∈

R | a ≤ r ≤ b and p(r)=0}. If N = roots cl int(p, a, b, sturm seq(p)) then

N ≥ 0 and there exists a bijection β : {0, 1,...,N − 1} → S.

The next step in the PVS development is to use the function roots cl int to

define a function roots interval that computes the number of roots of an integer

polynomial in any interval, whether it is closed, open, half open and half closed,

unbounded, and whether or not the polynomial has a multiple root at an endpoint of the interval. The precise definition of this function is omitted from this

description, because it is straightforward to define it directly in terms of the function roots cl int. The reader is referred to the PVS development Sturm in the

NASA PVS Library for the precise definition of this function. Basically, the number

of roots in the closure of the given interval is computed using roots cl int, and

this number is then adjusted upward or downward depending on whether the lower

bound or the upper bound of the interval is a root of the polynomial, and whether

the interval itself actually contains this bound. The function roots interval takes

as inputs an integer polynomial p, two extended real numbers a and b as inputs,

with a<b, as well as two Boolean values cont lb and cont ub, which refer to

whether the interval contains its lower or upper bound, respectively. That is, if

a = −10, b = 7, cont lb = true, and cont ub = false, then the corresponding

half open interval is [−10, 7). The following theorem has been proved in PVS and

can be found in the PVS development with the name number roots interval def.

10
Theorem 3. Let a, b ∈ R ̄, with a<b, p be an integer polynomial, and

S = {r ∈ R | a ≤ r ≤ b and p(r)=0 and

(r = a =⇒ cont lb = true) and

(r = b =⇒ cont ub = true)}.

If N = roots interval(p, a, b, cont lb, cont ub, sturm seq(p)) then N ≥ 0 and

there exists a bijection β : {0, 1,...,N − 1} → S.

It can now be noted that if a polynomial is always positive on a given interval, it

is trivial to prove that this is true using the function roots interval. All that must

be checked is that the function roots interval returns 0, so that the polynomial

has no roots on the interval, and that it is positive at any fixed point in the interval,

such as (a + b)/2 in the case where a, b ∈ R. Thus, the difficulty when determining

whether the polynomial p satisfies p(r) R 0 for all r in a given interval, when R is a

relation in {>, <, =, ≥, ≤} lies in the case when R is nonstrict, i.e., when R ∈ {≥, ≤}.

Moreover, if a decision procedure can be defined for when R is ≥, then it can be

defined for when R is ≤ as well by just replacing p with −p. Thus, the next step is to

define a specific decision procedure that determines whether an integer polynomial

is always nonnegative on a bounded closed interval. The PVS function nonneg int

takes as inputs an integer polynomial p, real numbers x and y (not extended real

numbers), and a list  of polynomials, which in practice is set to sturm seq(p). It

returns a Boolean, which is equal to true if and only if p(r) ≥ 0 for r in the closed,

bounded interval [x, y].

The function nonneg int works by recursively subdividing the interval [x, y] into

left and right halves, until each subinterval is small enough that it contains at most

one root of the polynomial p. Then, the polynomial p is evaluated at each endpoint

of each subinterval. One important fact that is used in the verification of this

function is that a continuous function with at most one root in a closed, bounded

interval is always nonnegative on that interval if and only if it is nonnegative at its

endpoints. This result follows directly from the intermediate value theorem. The

recursive function nonneg int is defined below.

nonneg int(p, x, y, ) ≡

if x>y then true

elsif x = y then p(x) ≥ 0

elsif roots cl int(p, x, y, ) ≤ 1 then p(x) ≥ 0 ∧ p(y) ≥ 0

else nonneg int(p, x,(x + y)/2, ) ∧

nonneg int(p,(x + y)/2, y, )

endif.

(13)

The function nonneg int is a decision procedure for nonnegativity on closed, bounded

intervals. However, it should be noted that a polynomial is always nonnegative on

any given interval if and only if it is nonnegative on the closure of that interval.

For instance, p is always nonnegative on the open interval (0, 1) if and only if it is

nonnegative on the closed interval [0, 1]. Thus, the function nonneg int can be used

11
as a decision procedure on any bounded interval, even if it is open. To extend this

function to unbounded intervals, a number is computed that can be guaranteed to

bound all of the roots of the integer polynomial p. This bound is used to reduce any

unbounded interval to a bounded interval on which the polynomial is nonnegative

if and only if it is nonnegative on the original unbounded interval.

root bound(p) ≡ max

⎛

⎝2,

deg(



p)−1

i=0

|coeff(p, i)|

|coeff(p, deg(p))|

⎞

⎠ . (14)

It can be proved that any root of p must lie in the bounded interval

(−root bound(p), root bound(p)).

The function nonneg int can now be used to define a function always nonneg that

takes as inputs an integer polynomial p and two extended (so possibly unbounded)

real numbers a and b with a<b. It returns true precisely when p is always

nonnegative on (a, b), which, as noted above, is equivalent to p being nonnegative on

the closure of this interval. Thus, this function does not need inputs that determine

whether the interval contains its endpoints, since that information is irrelevant to

nonnegativity.

always nonneg(p, a, b) ≡

let  = sturm seq(p),

M = root bound(p) in

if a = −∞ ∧ b = ∞ then nonneg int(p, a, b, )

elsif b = −∞ then nonneg int(p, min(−M, b − 1), b, )

elsif a = ∞ then nonneg int(p, a, max(M, a + 1), )

else nonneg int(p, −M,M, )

endif.

(15)

The following theorem has been proved in PVS and can be found in the PVS development under the name always nonnegative def.

Theorem 4. If a, b ∈ R ̄, with a<b, and p is an integer polynomial, then

always nonneg(p, a, b) = true

if and only if every real number x, with a ≤ x ≤ b, satisfies p(x) ≥ 0.

Note again that the continuity of p implies that this theorem is true if either (or

both) of the ≤ signs is replaced with a < sign.

The decision procedure that determines whether p(x) R 0 for all x in a given

interval, and where R is a relation in {>, <, =, ≥, ≤}, can now be defined as follows.

If R is either ≥ or ≤, then the function always nonneg is used for either p or −p

(respectively). Otherwise, it is simply checked that the function roots interval

returns 0 on the interval and that a given point r in the interval satisfies p(r) R 0.

12
The correctness of the procedure then follows from the intermediate value theorem.

As above, the interval here is represented by two extended real numbers a and b, as

well as two Boolean values cont lb and cont ub, which refer to whether the interval

contains a and b, respectively. The point r can be chosen in a number of ways. For

this development, it is simply assumed that there is a function choose with inputs a,

b, cont lb, and cont ub as inputs that picks a point in this interval. In the actual

PVS development, r is defined as the midpoint of the non-empty proper interval

[c, d], where c and d are real numbers and c<d. The numbers c, d are defined as

follows. If a and b are both real numbers then c = a and d = b; if both a and b

are infinite numbers, c = −1 and d = 1; if a is a real number and b is +∞, then

c = a and d = a + 1; otherwise, c = b − 1 and d = b. The input p to the function

compute poly sat, defined below, must be an integer polynomial.

compute poly sat(p, a, b, cont lb, cont ub, R) ≡

if R = (≥) then always nonneg(p, a, b)

elsif R = (≤) then always nonneg(−p, a, b)

else roots interval(p, a, b, cont lb, cont ub, sturm seq(p))

∧ R(p(choose(p, a, b, cont lb, cont ub)), 0)

endif.

(16)

The following theorem has been formally proved in PVS and can be found in the

PVS development under the name poly sat correct.

Theorem 5. If p is a nonzero integer polynomial, a and b are extended real numbers,

and R is a relation in {>, <, =, ≥, ≤}, then

compute poly sat(p, a, b, cont lb, cont ub, R) = true

if and only if p(x) R 0 for all real numbers x in I, where I is the interval corresponding to a, b and the Boolean values cont lb and cont ub.

Example 1. If a = 0, b = 3, cont lb = false, and cont ub = false, then the

corresponding interval is I is (0, 3). If p is the integer polynomial in §1 given by

x120 − 2x60 + 1, then p(x) ≥ 0 for all x in I. In fact, p factors as (x60 − 1)2. It can

be checked that

compute poly sat(p, 0, 3, false, false, ≥) = true,

compute poly sat(p, 0, 3, false, false, >) = false. (17)

5 Decision Procedure for Rational Polynomials

If the polynomial p has integer coefficients and R ∈ {>, <, =, ≥, ≤}, the function

compute poly sat, defined in §4, decides whether or not p(x) R 0 holds on any given

interval. Thus, if p has rational coefficients, the relation p(x) R 0 can be decided

on any given interval by multiplying this relational formula by the product of the

denominators of coefficients in p and then using the procedure compute poly sat

on the resulting polynomial.

13
In order to multiply the polynomial p by the product of all of its coefficients’

denominators, a function must be defined that computes the denominator of a given

coefficient. In a theorem prover like PVS, where rational numbers are built-in,

defining such as function is not straightforward. For example, in PVS, the terms 1

2

and 2

4 are indistinguishable. However, there are several ways in which a function that

computes the numerator and denominator of a rational number can be defined in

PVS. For example, the PVS function that computes the rational decomposition uses

a recursive function called compute pos rat that continually simplifies a positive

rational number, making its numerator and denominator smaller and removing the

integer part at each steps. In the final recursive call, the function is called on an

integer, making the answer trivial. The function takes a positive rational number

and return as a pair of positive integers, the first being the numerator and the second

being the denominator.

compute pos rat(r) ≡

if floor(r) = r then (floor(r), 1)

elsif r > 1 then let (a, b) = compute pos rat(r − floor(r)) in

(b · floor(r) + a, b)

else let (a, b) = compute pos rat(1/r) in (b, a) endif.

(18)

The function compute pos rat can then be used to define two other functions

numerator and denominator, which are simply given by the first and second componet of the output of compute pos rat. It has been proved in PVS that if r is

any positive rational number, then r = numerator(r)/denominator(r). A challenging part of the proof is showing that the function compute pos rat terminates. In

PVS, showing termination requires that the function be given a measure, which is

a function on the inputs of the function that returns a natural number and strictly

decreases in value every time the function is called recursively. In PVS, the measure

function on r is defined by

pos rat meas(r) ≡ if r < 1 then 10g else 10g − 1 endif, (19)

where g is the least integer such that there exists positive rational numbers a and b

with r = a/b and g = a + b.

The final step before defining a decision procedure for rational polynomials is

to define a function that takes as an input a rational polynomial p and returns an

integer polynomial that is a positive multiple of the first. This integer polynomial

is obtained by multiplying p by the product of the denominators of its coefficients.

Each coefficient supplies one integer to this product. If c is one of the coefficients,

then this number is 1 if c is zero and it is |denominator(c)| otherwise. A function

rat2poly is defined in PVS that takes a rational polynomial and coverts it to an

integer polynomial. This process works by recursively considering each coefficient

of the polynomial. At each step, it multiplies the polynomial by the denominator

of the coefficient in question, and it also stores the current greatest common divisor

of all resulting integer coefficients that it has simplified so far in the recursion. At

the end of the recursion, all of the coefficients are divided by this greatest common

14
divisor to simplify the answer. The fact that rat2poly(p) is a positive multiple of

p is specified in the output type of the function rat2poly. Using this function, the

decision procedure for rational polynomials is defined below.

compute poly sat rat(p, a, b, cont lb, cont ub, R) ≡

compute poly sat(rat2poly(p), a, b, cont lb, cont ub, R). (20)

The following theorem has been formally proved in PVS and is found in the PVS

development under the name poly sat rat correct.

Theorem 6. Let p be a nonzero rational polynomial, a, b ∈ R ̄, and R be a relation

in {>, <, =, ≥, ≤}.

compute poly sat rat(p, a, b, cont lb, cont ub, R) = true

if and only if p(x) R 0 for all real numbers x in I, where I is the interval corresponding to a, b and the Boolean values cont lb and cont ub.

Example 2. If a = −∞, b = 3, cont lb = false, and cont ub = false, then the

corresponding interval I is (−∞, 3). If p is the rational polynomial in §1 given by

x120 − 2

3x60 + 19 , then p(x) ≥ 0 for all x in I. In fact, p factors as (x60 − 13 )2. It can

be seen that

compute poly sat rat(p, −∞, 3, false, false, ≥) = true,

compute poly sat rat(p, −∞, 3, false, false, >) = false. (21)

6 Automated Strategy

The decision procedure implemented by the function compute poly sat rat and its

correctness property (Theorem 6) can be used to discharge, in a methodical way,

PVS sequents involving univariate polynomial inequalities.

In a nutshell, a PVS sequent is a logical judgement of the form Γ  Δ, where

Γ, called the antecedent, and Δ, called the consequent, are lists of logical formulas.

The intuitive meaning of a sequent is that from the conjunction of formulas in Γ, the

disjunction of formulas in Δ can be deduced. PVS proof commands are logical rules

that transforms a sequent into a set of sequents, with the objective of producing

sequents where the formula false appears in the antecedent or the formula true

appears in the consequent. When all sequents generated by a proof command are

of one of these forms, the initial sequent is discharged. Hence, a proof in PVS can

be represented by a tree of proof commands that discharge an initial sequent.

The following example describes a PVS proof of the polynomial inequality given

in Example 2 using the development presented in §5.

Example 3. Consider the sequent below displayed in a vertical layout. Formulas in

a PVS sequent are numbered using the notation {n}, where n < 0 in the antecedent

and n > 0 in the consequent. In this sequent, x is a free real variable (actually, a

Skolem real constant in PVS terminology).

{−1} x < 3



{ 1 } x120 − 2

3x60 + 19 ≥ 0

15
1. The first step in the PVS proof is to execute a proof command that instantiates

Theorem poly sat rat correct with a = −∞, b = 3, cont lb = false,

cont ub = false, and p = p, where p is a list of rational numbers that is

0 everywhere except in the positions 0, 60, and 120, where it has the values

1

9 , −23 , and 1, respectively. It should be noted that p is a PVS data structure

that represents the polynomial p(x) = x120 − 2

3x60 + 19 . The notation p(x) will

be used to represent the evaluation of the polynomial represented by the data

structure p on x. That proof command yields the following sequent.

{−1} compute poly sat rat(p, −∞, 3, false, false, ≥) ⇐⇒

(∀(x: R): x ∈ (−∞, 3) =⇒ p(x) ≥ 0)

{−2} x < 3



{ 1 } x120 − 2

3x60 + 19 ≥ 0

2. Next, a proof rule that evaluates

compute poly sat rat(p, −∞, 3, false, false, ≥)

is executed. This evaluation, which only involves computable functions and

ground terms, is efficiently performed by the PVS ground evaluator. The following sequent is obtained.

{−1} true ⇐⇒ (∀(x: R): x ∈ (−∞, 3) =⇒ p(x) ≥ 0)

{−2} x < 3



{ 1 } x120 − 2

3x60 + 19 ≥ 0

3. Sequent formula {−1} can be easily reduced using propositional simplification

to

{−1} ∀(x: R): x ∈ (−∞, 3) =⇒ p(x) ≥ 0

{−2} x < 3



{ 1 } x120 − 2

3x60 + 19 ≥ 0

4. Next, the quantified variable x in sequent formula {−1} is instantiated with

the Skolem constant x appearing in sequent formulas {−2} and {1}. This

instantiation yields the sequent.

{−1} x ∈ (−∞, 3) =⇒ p(x) ≥ 0

{−2} x < 3



{ 1 } x120 − 2

3x60 + 19 ≥ 0

5. The elimination of the implication in sequent formula {−1} yields two sequents.

{−1} x < 3



{ 1 } x ∈ (−∞, 3)

{ 2 } x120 − 2

3x60 + 19 ≥ 0

{−1} p(x) ≥ 0

{−2} x < 3



{ 1 } x120 − 2

3x60 + 19 ≥ 0

16
6. The two sequents generated in the previous step can be easily discharged by

unfolding the definitions of the operations “∈” and “p(x)”, respectively.

The PVS proof in Example 3 illustrates a well-known technique in theorem

proving known as computational reflection [12]. Given an object theory, such as

the theory of univariate rational polynomials, a computable target theory is defined

such that elements of the object theory are embedded using data structures in

the target theory, e.g., formulas x ∈ (−∞, 3) and p(x) ≥ 0 are embeddings of

x < 3 and x120 − 2

3x60 + 19 ≥ 0, respectively. Functions and theorems, such as

compute poly sat rat and poly sat rat correct, that relate both theories enable

the proof of formulas in the object theory by evaluating terms in the target theory,

e.g., Step 2 in Example 3.

Computational reflection is particularly well-adapted for the development of automated proof strategies in interactive theorem provers. First, it produces proofs

whose size is independent of the size of the initial sequent. For instance, the method

in Example 3 can be used for discharging inequalities on any univariate polynomial

of any degree and any number of monomials. Second, proofs by computational reflection are small since the most involved logical steps are done once for all in the

proof of theorems such as poly sat rat correct. Finally, proofs by computational

reflection are efficient since they depend on computation rather than on deduction.

A PVS strategy is a procedure that, when executed by the theorem prover’s

Lisp engine, produces a PVS proof, i.e., a tree of proof commands. The PVS

strategy language includes combinators for sequencing, branching, and backtracking

of proof commands. The language also provides mechanisms to inspect the internal

representation of PVS syntactic elements within a proof context. For instance,

in the PVS specification language, the expression x120 − 2

3x60 + 19 ≥ 0 is just a

Boolean expression. In the strategy language, it is possible to parse the internal

Lisp representation of this expression and determine that it is a real order relation

whose left hand side is a univariate rational polynomial of degree 120 over variable

x with monomials x120, −2

3x60, and 19 . PVS strategies, as tactics in the LCF-style

of interactive theorem provers, preserve the logical consistency of the proof system.

That is, any strategy within a proof can be unfolded in a tree of proof commands

that only includes basic PVS proof rules.

The development presented in this paper includes a PVS strategy called sturm

that implements an enhanced version of the method described in Example 3. Using this strategy, the initial sequent in this example is automatically discharged in

less than one second. More generally, the strategy sturm automatically discharges

sequents having one of the following forms

1. X1,...,Xm, Γ  p1(x) ✸ p2(x), Δ

2. X1,...,Xm, p1(x) ✸p2(x), Γ  Δ

3. Γ  ∀(x: T): X1 ∧ ... ∧ Xm =⇒ p1(x) ✸ p2(x), Δ

4. ∀(x: T): X1 ∧ ... ∧ Xm =⇒ p1(x) ✸ p2(x), Γ  Δ

5. Γ  ∃(x: T): X1 ∧ ... ∧ Xm ∧ p1(x) ✸p2(x), Δ

17
6. ∃(x: T): X1 ∧ ... ∧ Xm ∧ p1(x) ✸p2(x), Γ  Δ

where

• T is a subtype of R,

• Γ and Δ are arbitrary lists of formulas, which are ignored by the strategy,

• m ≥ 0 and for 1 ≤ j ≤ m, Xm denotes a Boolean expression of one of the

forms a ≺ x, x ≺ a, |x| ≺ a, or x ∈ [a, b], where ≺∈ {<, ≤}, and a, b are

numerical rational constants,

• p1(x) and p2(x) denote univariate rational polynomial such that the polynomial p(x) = p1(x) − p2(x) is not a constant.

• ✸ is a real order in {<, ≤, >, ≥, =} and ✸is a real order in {<, ≤, >, ≥, =}.

The strategy sturm works on a formula of interest, which is the underlined formula in the forms above. By default, the strategy assumes that the formula of

interest is the first formula in the consequent, but the user can specify a different

formula though an optional parameter in the strategy. As specified above, the formula of interest can appear in the antecedent or in the consequent, be nonquantified

or quantified, and the quantifier can be existential or universal.

The strategy proceeds as follows. First, it examines the sequent and determines

whether or not it has one of the supported forms. If this is not the case, the strategy prints an error message and does nothing else. If the sequent is supported,

the strategy parses the formula of interest and constructs a PVS data structure

p that represents the univariate polynomial p(x) = p1(x) − p2(x). This part of

the strategy uses a function developed by B. Di Vito (NASA) for parsing PVS

polynomial expressions [22]. The parse does not assume any particular polynomial

normal form. In particular, real expressions such as (x-1)2, (x-1)*(x-1), and

x*x-2*x+1 are all parsed into the same data structure p. Next, the strategy computes PVS data structures for instantiating the variables a, b, cont lb, cont ub in

Theorem poly sat rat correct. These data structures are constructed by examining the relational formulas Xj , 1 ≤ j ≤ m, and, in the case of quantified forms, the

type T. It is noted that the correctness of the strategy is not compromised by the

constructions of these objects. Indeed, as shown in Example 3, the strategy checks

that the semantics of the original expressions and their data structure embeddings

coincide.

The instantiation of relation R in Theorem poly sat rat correct depends on

the form of the sequent. In the case of forms 1, 3, and 4, R is instantiated with

✸. Otherwise, R is instantiated with ¬✸

. From this point on, the strategy follows the method described in Example 3. Some minor modifications are needed

for the cases where the formula of interest is quantified. Sequents of the forms 3

and 6 can be transformed into sequents of the forms 1 and 2 by introducing the

quantified variable (this process is called Skolemization, in PVS terminology) and

propositional simplification. Sequents of the forms 4 and 5 use the fact that Theorem poly sat rat correct is a double implication. In this case, the evaluation

18
of the function compute poly sat rat returns false and the proof proceeds accordingly. Finally, if the sequent holds, the strategy succeeds and the sequent is

discharged. If the sequent does not hold, the strategy prints a message stating that

the sequent is not provable.

Example 4. Turan’s inequality for Legendre’s polynomials is a theorem that states

that pn(x)2 > pn−1(x) · pn+1(x) for all x such that −1 <x< 1, where pj is the j-th

Legendre polynomial (of degree j). For instance, the 10-th Legendre polynomial is

given by 1

256 · (46189x10 − 109395x8 + 90090x6 − 30030x4 + 3465x2 − 63). Turan’s

inequality for n = 10 easily reduces to showing that

3969

65536 +

63063

4096 x6 +

1792791

4096 x10 +

3002285

4096 x18 +

6600165

4096 x14 − 7276565536

x4

− 3558555

32768 x8 − 1020776965536

x20 − 35043645

32768 x12 − 9585189965536

x16 > 0,

for all x ∈ (−1, 1). The strategy sturm in PVS automatically discharges this theorem

in less than a second. In fact, the strategy is able to prove that the polynomial on

the left hand side of this inequality is positive over this interval when it is cubed and

then has degree 60, a proof that takes the decision procedure less than 10 seconds to

run.

The Appendix includes examples of PVS formulas, including Example 3 and

Example 4, that are automatically discharged by the strategy sturm.

7 Related Work

The authors are aware of three bodies of work that are closely related to this paper

in their use of Sturm’s Theorem for nonlinear reasoning in a theorem prover.

The first related work was accomplished by John Harrison in the 90’s [13]. He

proved Sturm’s Theorem in HOL and used it there for root isolation. The intended

application of that work was to guarantee error bounds of polynomial approximations to transcendental functions. Root isolation was used on the derivative of the

error to find the places where an error may be maximized. In that work, Sturm’s

Theorem was proved in the case where the polynomial is square free, which simplifies the proof. In that paper, Harrison writes “[. . . ] we would prefer the polynomial

to have no multiple real roots [. . . ] Sturm’s Theorem is easier to prove for polynomials without multiple real roots - this is actually the only form we have proved in

HOL.” Harrison was not interested in proving statements like p(x) ≥ 0, so a decision

procedure for such problems was out of the scope of his paper.

The recent work by Eberl [8] is the most similar to this paper. It was apparently

being completed concurrently with this work. Eberl completed a formal proof in

Isabelle/HOL of Sturm’s Theorem and used that to define a proof method sturm in

Isabelle/HOL for solving simply quantified polynomial inequalities similar to those

solved by PVS’s strategy sturm. In the case of non-strict universally-quantified

inequalities, Eberl’s proof method relies on unverified ML code to generate the

interval splitting as a witness. In practice, the use of ML code improves efficiency

19
and does not compromises soundness. However, it may compromise completeness.

In contrast, the PVS’s strategy sturm relies on a PVS algorithm that is proven

to be correct and complete. Another distinction between the current paper and

the work of Eberl is in the use of pseudo division instead of regular division. PVS

handles A final point worth noting regarding the differences between the work in

this paper and that by Eberl is that the proofs of Sturm’s Theorem are different.

Eberl proves Sturm’s Theorem in the square free case, similar to Harrison’s proof in

HOL [13]. In the non-square free case, the proof proceeds by dividing each term in

the remainder sequence by the greatest common divisor of the original polynomial

and its derivative. The resulting sequence is not a Sturm sequence in the standard

sense, but it maintains similar properties regarding root counting. The proof in

the PVS development presented in this paper is a direct approach that considers

the highest power of a linear divisor that divides each polynomial in the sequence

and analyzes whether the polynomial swaps signs at the corresponding root. The

distinctions between these proof methods primarily amount to user preference.

Sturm sequences and several speed enhancements such as the use of pseudo

division are implemented in the SMT solver Z3 [6]. That implementation was the

inspiration for the work presented in this paper.2 Z3 is a highly efficient tool. In

some cases, for example when a formula is satisfiable, Z3 can produce models, which

can be understood as proof certificates for existential formulas. However, in general,

Z3 statements are not supported by formal proofs. Hence, in formal verifications

efforts, Z3 is used as a trusted oracle.

Z3 is used as an external algebraic decision method (EADM) in Metitarski, a

theorem prover for real numbers [1]. In a recent work, Denman and Mu ̃noz [7]

developed the PVS proof rule metit that integrates Metitarski as an external oracle

into PVS theorem prover. In contrast to the work presented in this paper, metit

is not implemented as a proof producing strategy. Nevertheless, the integration of

Metitarski/Z3 in PVS, while unproven, is quite useful and has helped the authors

to check results that were impossible with previous strategies in PVS.

There is a much larger collection of works on the general problem of reasoning about nonlinear arithmetic. Sophisticated implementations of the Cylindrical

Algebraic Decomposition (CAD) [2] procedure are available in the Redlog system3

and in the QEPCAD library.4 The systems RealPaver [11] and dReal [9] integrate

powerful methods based on interval constraint propagation [11]. MetiTarski [1] and

RAHD (Real Algebra in High Dimensions) [25] are specialized theorem provers for

the theory of real closed fields. MetiTarski is designed to prove universally quantified inequalities involving real-valued functions such as trascendental functions.

RAHD combines several decision methods for the existential theory of real closed

fields. Both systems use a CAD procedure for quantifier elimination among many

other proof strategies.

Proof tactics that implement H ̈ormander’s quantifier elimination method are

available in Coq and HOL Light [18, 19]. These tactics are theoretically interesting.

2The authors would like to give credit to Leonardo de Moura for directing their work in this

direction and advising them in this effort. 3http://redlog.dolzmann.de. 4http://www.usna.edu/cs/~qepcad/B/QEPCAD.html.

20
H ̈ormander’s method is known to be more inefficient than CAD. An implementation

of CAD that will eventually yield a proof producing tactic is available in Coq [17].

Another approach to solve multivariate polynomial inequalities in theorem provers

is based on polynomial sum of square (SOS) decompositions through semidefinite

programming. Such an approach has been implemented in HOL Light [14] and

seems to be more promising than quantifier elimination for polynomials with many

variables. Semidefinite programming is a somewhat complicated numerical procedure that is usually implemented with floating point numbers. Because of numerical

approximation errors, it is difficult to integrate this method into theorem provers.

Recent developments in SOS address this issue by producing rational polynomial

decompositions [15, 21]. Proof producing strategies for proving real-number properties based on interval arithmetic and branch and bound methods are available in

PVS [4], Coq [20], HOL Light [26]. As stated in the introduction, the authors have

developed semi-decision procedures for multivariate polynomials, based on Bernstein

polynomials [22], and Boolean expressions involving real-valued functions, based on

interval arithmetic [23]. Those algorithms are quite powerful and can prove tight

bounds on complex polynomials with up to 16 variables and degree 4. However,

since Bernstein polynomials and interval arithmetic yield numerical approximation

methods, they cannot always prove exact bounds on a polynomial. That is, in general, they cannot prove that p(x) ≥ 0 for x in a fixed, bounded interval unless it is

also true that p(x) > c for some positive c. This limitation is key to the authors’

development of Sturm’s Theorem in PVS.

8 Conclusion

This paper presented a formalization Sturm’s Theorem in PVS along with a decision procedure for deciding the sign of univariate rational polynomials where the

polynomial variable ranges over an interval, which can be any connected set of real

numbers. The decision procedure, which is is shown to be complete and correct in

PVS, is used to implement a proof strategy for automatically discharging sequents

involving univariate polynomial inequalities. The correctness of this strategy only

depends on the correctness of PVS deduction engine. The strategy is based on

computational reflection, which is a theorem proving technique for building efficient

strategies. Although the strategy uses concrete data structures for representing

polynomials and infinite numbers, these data structures are invisible to the user.

The strategy can be used to discharge sequents involving native PVS real number

expressions.

In its standard form, Sturm’s Theorem is only applied when the end points of the

interval are not multiple roots and the polynomial inequality is strict. The authors

address these limitations by using a branching algorithm that progressively subdivides the interval, noting that a polynomial with at most one root in a given interval

is nonnegative on that interval if and only if it is nonnegative at both endpoints.

Furthermore, the PVS development presented in this paper uses the remainder of

polynomials after pseudo division rather than standard division. This technique

ensures that integer polynomials beget more integer polynomials. Moreover, be21
fore working with a polynomial, the PVS decision procedure first divides it by the

greatest common divisor of its coefficients. These features substantially improve the

efficiency of the decision procedure.

The authors’s main motivation for developing Sturm’s Theorem in PVS comes

from NASA’s verification effort on aircraft separation assurance systems.5 The

PVS proofs of correctness of these system are nontrivial and require considerable

algebraic manipulations. In one particular case of an actual verification effort, a

176-line sequent involving a 16-variable polynomial was generated. That sequent

was automatically checked using a PVS proof rule that integrates Metitarski and Z3

as trusted external oracles into the PVS theorem prover [7]. This example shows the

usefulness of techniques based on Sturm sequences, which are implemented in Z3.

The ultimate goal of the work presented here is to build a similar algorithm to that

used in Z3, directly in PVS with a formal proof of its correctness. The work with

univariate polynomials presented in this paper is a first step in that direction. The

next step is to define an algorithm that not only handles multivariate polynomials,

but also handles arbitrary Boolean expressions involving those polynomials. This

work will be guided by the algorithm developed by de Moura in Z3, but the authors

expect subtleties to arise, since it will have to be designed with its formal verification

in mind. That is, having to prove its correctness may change the way that the

algorithm is defined in PVS.

References

1. Behzad Akbarpour and Lawrance C. Paulson. MetiTarski: An automatic theorem prover for real-valued special functions. Journal of Automated Reasoning,

44(3):175–205, 2010.

2. George Collins. Quantifier elimination for real closed fields by cylindrical algebraic decomposition. In Second GI Conference on Automata Theory and Formal

Languages, volume 33 of Lecture Notes in Computer Science, pages 134–183,

Kaiserslautern, 1975. Springer-Verlag.

3. Luis G. Crespo, C ́esar A. Mu ̃noz, Anthony J. Narkawicz, Sean P. Kenny, and

Daniel P. Giesy. Uncertainty analysis via failure domain characterization: Polynomial requirement functions. In Proceedings of European Safety and Reliability

Conference, Troyes, France, September 2011.

4. Marc Daumas, David Lester, and C ́esar Mu ̃noz. Verified real number calculations: A library for interval arithmetic. IEEE Transactions on Computers,

58(2):1–12, February 2009.

5. Florent de Dinechin, Christoph Lauter, and Guillaume Melquiond. Certifying

the floating-point implementation of an elementary function using Gappa. IEEE

Transactions on Computers, 60(2):242–253, February 2011.

5http://shemesh.larc.nasa.gov/fm/fm-atm-cdr.html

22
6. Leonardo de Moura and Grant Passmore. Computation in real closed infinitesimal and transcendental extensions of the rationals. In Automated Deduction - CADE-24, 24th International Conference on Automated Deduction, Lake

Placid, New York, June 9-14, 2013, Proceedings, 2013.

7. William Denman and C ́esar Mu ̃noz. Automated real proving in PVS via metitarski, 2014. Accepted for publication at 19th International Symposium on

Formal Methods (FM 2014).

8. Manuel Eberl. A formal proof of Sturm’s theorem in Isabelle/HOL. Manuscript.

Available from http://afp.sf.net/entries/Sturm_Sequences.shtml, 2014.

9. Sicun Gao, Soonho Kong, and Edmund M. Clarke. dReal: An SMT solver for

nonlinear theories over the reals. In Maria Paola Bonacina, editor, Automated

Deduction - CADE-24 - 24th International Conference on Automated Deduction,

Lake Placid, NY, USA, June 9-14, 2013. Proceedings, volume 7898 of Lecture

Notes in Computer Science, pages 208–214. Springer, 2013.

10. J. Garloff. Application of Bernstein expansion to the solution of control problems. Reliable Computing, 6:303–320, 2000.

11. Laurent Granvilliers and Fr ́ed ́eric Benhamou. RealPaver: An interval solver

using constraint satisfaction techniques. ACM Transactions on Mathematical

Software, 32(1):138–156, March 2006.

12. John Harrison. Metatheory and reflection in theorem proving: A survey and

critique. Technical Report CRC-053, SRI Cambridge, Millers Yard, Cambridge,

UK, 1995.

13. John Harrison. Verifying the accuracy of polynomial approximations in HOL.

In Elsa L. Gunter and Amy Felty, editors, Theorem Proving in Higher Order

Logics: 10th International Conference, TPHOLs’97, volume 1275 of Lecture

Notes in Computer Science, pages 137–152, Murray Hill, NJ, 1997. SpringerVerlag.

14. John Harrison. Verifying nonlinear real formulas via sums of squares. In Theorem Proving in Higher Order Logics, volume 4732 of Lecture Notes in Computer

Science, pages 102–118. Springer, 2007.

15. Erich L. Kaltofen, Bin Li, Zhengfeng Yang, and Lihong Zhi. Exact certification

in global polynomial optimization via sums-of-squares of rational functions with

rational coefficients. In Lorenzo Robbiano and John Abbott, editors, Approximate Commutative Algebra, Texts and Monographs in Symbolic Computation.

Springer Vienna, 2010.

16. James Kuchar and Lee Yang. A review of conflict detection and resolution

modeling methods. IEEE Transactions on Intelligent Transportation Systems,

1(4):179–189, December 2000.

23
17. Assia Mahboubi. Implementing the cylindrical algebraic decomposition within

the Coq system. Mathematical Structures in Computer Science, 17(1):99–127,

February 2007.

18. Assia Mahboubi and Lo ̈ıc Pottier. Elimination des quantificateurs sur les r ́eels

en Coq. In Journ ́ees Francophone des Langages Applicatifs (JFLA), 2002.

19. Sean McLaughlin and John Harrison. A proof-producing decision procedure

for real arithmetic. In Robert Nieuwenhuis, editor, Proceedings of the 20th

International Conference on Automated Deduction, proceedings, volume 3632 of

Lecture Notes in Computer Science, pages 295–314, 2005.

20. Guillaume Melquiond. Proving bounds on real-valued functions with computations. In Alessandro Armando, Peter Baumgartner, and Gilles Dowek, editors,

Automated Reasoning, 4th International Joint Conference, IJCAR 2008, Sydney, Australia, August 12-15, 2008, Proceedings, volume 5195 of Lecture Notes

in Computer Science, pages 2–17. Springer, 2008.

21. David Monniaux and Pierre Corbineau. On the generation of Positivstellensatz

witnesses in degenerate cases. In Proceedings of Interactive Theorem Proving

(ITP). Lecture Notes in Computer Science, 2011.

22. C ́esar Mu ̃noz and Anthony Narkawicz. Formalization of a representation of

Bernstein polynomials and applications to global optimization. Journal of Automated Reasoning, 51(2):151–196, August 2013.

23. Anthony Narkawicz and C ́esar Mu ̃noz. A formally verified generic branching

algorithm for global optimization. In Ernie Cohen and Andrey Rybalchenko,

editors, Fifth Working Conference on Verified Software: Theories, Tools and

Experiments (VSTTE 2013), volume 8164 of Lecture Notes in Computer Science, pages 326–343. Springer, 2014.

24. Sam Owre, John Rushby, and Natarajan Shankar. PVS: A prototype verification system. In Deepak Kapur, editor, Proceeding of the 11th International

Conference on Automated Deduction (CADE), volume 607 of Lecture Notes in

Artificial Intelligence, pages 748–752. Springer, June 1992.

25. Grant Olney Passmore and Paul B. Jackson. Combined decision techniques

for the existential theory of the reals. In L. Dixon, editor, Proceedings of Calculemus/Mathematical Knowledge Managment, number 5625 in LNAI, pages

122–137. Springer-Verlag, 2009.

26. Alexey Solovyev and Thomas C. Hales. Formal verification of nonlinear inequalities with Taylor interval approximations. In Guillaume Brat, Neha Rungta, and

Arnaud Venet, editors, Proceedings of the 5th International Symposium NASA

Formal Methods, volume 7871 of Lecture Notes in Computer Science, pages

383–397, 2013.

27. Frank Sottile. Chapter 2: Real solutions to univariate polynomials. Course

Notes.

24
28. C. Sturm. M ́emoire sur la r ́esolution des  ́equations num ́eriques. In JeanClaude Pont, editor, Collected Works of Charles Fran ̧cois Sturm, pages 345–390.

Birkhuser Basel, 2009.

25
Appendix A

PVS Examples of Strategy sturm

examples : THEORY

BEGIN

IMPORTING Sturm@strategies

x : VAR real

Example_0 : LEMMA x*(1-x) <= 1/4

%|- Example_0 : PROOF (sturm) QED

Example_1 : LEMMA

0 < x AND x < 3 IMPLIES x^120 - 2*x^60 + 1 >= 0

%|- Example_1 : PROOF (sturm) QED

Example_2_3 : LEMMA

x < 3 IMPLIES x^120 - (2/3)*x^60 + 1/9 >= 0

%|- Example_2_3 : PROOF (sturm) QED

Example_4: LEMMA

abs(x) < 1 IMPLIES

3969/65536 + 63063/4096 * x^6 + 1792791/4096 * x^10 +

3002285/4096 * x^18 + 6600165/4096 * x^14

- 72765/65536 * x^4 - 3558555/32768 * x^8

- 10207769/65536 * x^20 - 35043645/32768 * x^12

- 95851899/65536 * x^16 > 0

%|- Example_4 : PROOF (sturm) QED

Example_5: LEMMA

abs(x) < 1 IMPLIES

(3969/65536 + 63063/4096 * x^6 + 1792791/4096 * x^10 +

3002285/4096 * x^18 + 6600165/4096 * x^14

- 72765/65536 * x^4 - 3558555/32768 * x^8

- 10207769/65536 * x^20 - 35043645/32768 * x^12

- 95851899/65536 * x^16)^3 > 0

%|- Example_5 : PROOF (sturm) QED

Example_6: LEMMA

EXISTS (x) : -1 <= x AND x <= 0 AND

32*x^5 - 160*x^3 + 120*x > 16*x^4 - 48*x^2 +12

%|- Example_6 : PROOF (sturm) QED

END examples

26
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

A Formally-Verified Decision Procedure for Univariate Polynomial 

Computation Based on Sturm's Theorem 

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Narkawicz, Anthony J.; Munoz, Cesar A.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20489

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 64

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

Sturm's Theorem is a well-known result in real algebraic geometry that provides a function that computes the number of roots 

of a univariate polynomial in a semiopen interval. This paper presents a formalization of this theorem in the PVS theorem 

prover, as well as a decision procedure that checks whether a polynomial is always positive, nonnegative, nonzero, negative, 

or nonpositive on any input interval. The soundness and completeness of the decision procedure is proven in PVS. The 

procedure and its correctness properties enable the implementation of a PVS strategy for automatically proving existential and 

universal univariate polynomial inequalities. Since the decision procedure is formally verified in PVS, the soundness of the 

strategy depends solely on the internal logic of PVS rather than on an external oracle. The procedure itself uses a combination

of Sturm's Theorem, an interval bisection procedure, and the fact that a polynomial with exactly one root in a bounded interval

is always nonnegative on that interval if and only if it is nonnegative at both endpoints. 

15. SUBJECT TERMS

Polynomial inequalities; Prototype verification system; Sturm's theorem; non-linear arithmetic decision procedure

18. NUMBER

 OF 

 PAGES

32

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

 534723.02.15.07 

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TM-2014-218548

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 11 - 2014