Evaluating Winding Numbers and Counting Complex Roots Through Cauchy Indices in Isabelle/HOL

In complex analysis, the winding number measures the number of times a path (counter-clockwise) winds around a point, while the Cauchy index can approximate how the path winds. We formalise this approximation in the Isabelle theorem prover, and provide a tactic to evaluate winding numbers through Cauchy indices. By further combining this approximation with the argument principle, we are able to make use of remainder sequences to effectively count the number of complex roots of a polynomial within some domains, such as a rectangular box and a half-plane.


Introduction
The winding number, given by measures how the path γ winds around the complex point z. It is an important object in complex analysis, and its evaluation is ubiquitous among analytic proofs. However, when formally evaluating the winding number in proof assistants such as Isabelle/HOL and HOL Light, unexpected difficulties arise, as pointed out by Harrison [8] and Li et al. [14]. To address this problem, we formalise a theory of the Cauchy index on the complex plane, thereby approximating how the path winds. When the path is a cycle and comprises line segments and parts of circles, we can now evaluate the winding number by calculating Cauchy indices along those sub-paths.
In addition, by further combining our previous formalisation of the argument principle [14] (which associates the winding number with the number of complex roots), we build effective procedures to count the complex roots of a polynomial within some domains, such as a rectangle box and a half-plane.
In short, the main contributions of this paper are -a novel tactic to enable users to evaluate the winding number through Cauchy indices, -and novel verified procedures to count complex roots of a polynomial.
The Isabelle sources of this paper are available from the Archive of Formal Proofs [11,12]. Formulations in this paper, such as the definition of the Cauchy index and statements of some key lemmas, mainly follow Rahman and Schmeisser's book [19,Chapter 11] and Eisermann's paper [6]. Nevertheless, we were still obliged to devise some proofs on our own, as discussed later.
This paper continues as follows: we start with a motivating example (Sect. 2) to explain the difficulty of formal evaluation of the winding number in Isabelle/HOL. We then present an intuitive description of the link between the winding number and the Cauchy indices (Sect. 3), which is then developed formally (Sect. 4). Next, we present verified procedures that count the number of complex roots in a domain (Sect. 5), along with some limitations (Sect. 6) and make some general remarks on the formalisation (Sect. 7). Finally, we discuss related work (Sect. 8) and present conclusions (Sect. 9).

A Motivating Example
In the formalisation of Cauchy's residue theorem [14], we demonstrated an application of this theorem to formally evaluate an improper integral in Isabelle/HOL: The idea is to embed this integral into the complex plane, and, as illustrated in Fig. 1, to construct a linear path L r from −r to r and a semi-circular path C r centred at 0 with radius r > 1: Next, by letting and r → ∞, we can derive (1) through the following steps: Here L r + C r is formed by appending C r to the end of L r , and Res( f , i) is the residue of f at i. Equation (3) is because C r f = 0 as r → ∞. The application of the residue theorem is within (4); we exploit the fact that i and −i are the only two singularities of f over the complex plane, since While carrying out the formal proofs of (5), surprisingly, the most troublesome part of the proof is to evaluate the winding numbers: n(L r + C r , i) = 1 ( 6 ) n(L r + C r , −i) = 0.
Equations (6) and (7) are straightforward to humans, as it can be seen from Fig. 1 that L r + C r passes counterclockwise around the point i exactly one time, and around −i zero times. However, formally deriving these facts was non-trivial.
Example 1 (Proof of n(L r +C r , i) = 1) We defined an auxiliary semi-circular path C r where C r (t) = re iπ(t+1) for t ∈ [0, 1] as can be seen in Fig. 1a. As C r + C r forms a (full) circular path with i lying inside the circle, we had n(C r + C r , i) = 1.
In addition, we further proved that C r + C r and L r + C r are homotopic on the space of the complex plane except for the point i (i.e., on C − {i}), and hence n(L r + C r , i) = n(C r + C r , i) (9) by using the following Isabelle lemma: Lemma 1 (winding_number_homotopic_paths) fixes z::complex and γ 1 γ 2 ::"real ⇒ complex" assumes "homotopic_paths (-{z}) γ 1 γ 2 " shows "winding_number γ 1 z = winding_number γ 2 z" where winding_number γ 1 z encodes the winding number of γ 1 around z: n(γ 1 , z), and homotopic_paths encodes the homotopic proposition between two paths. Putting (8) and (9) together yields n(L r + C r , i) = 1, which concludes the whole proof.
Example 2 (Proof of n(L r + C r , −i) = 0) We started by defining a ray L r starting from −i and pointing towards the negative infinity of the imaginary axis: as illustrated in Fig. 1b. Subsequently, we showed that and then applied the following lemma in Isabelle Lemma 2 (winding_number_less_1) fixes z w::complex and γ ::"real ⇒ complex" assumes "valid_path γ " and "z / ∈ path_image γ " and "w = z" and not_intersection:" a::real. 0 < a ⇒ z + a*(w -z) / ∈ path_image γ " shows "|Re(winding_number γ z)| < 1" where -valid_path γ assumes that γ is piecewise continuously differentiable on [0, 1], z / ∈ path_image γ asserts that z is not on the path γ , -the assumption not_intersection asserts that the ray starting at z ∈ C and through w ∈ C ({z + a(w − z) | a > 0}) does not intersect with γ -for all a > 0, z + a(w − z) does not lie on γ .
Note that the real part of a winding number Re(n(γ , z)) measures the degree of the winding: in case of γ winding around z counterclockwise for exactly one turn, we have n(γ , z) = Re(n(γ , z)) = 1. Essentially, Lemma 2 claims that a path γ can only wind around z for less than one turn, |Re(n(γ , z))| < 1, if there is a ray starting at z and not intersecting with γ . Joining Lemma 2 with (10) leads to Moreover, as L r + C r is a closed path, By combining (11) and (12), we managed to derive n(L r + C r , −i) = 0. As can be observed in Examples 1 and 2, our proofs of n(L r + C r , i) = 1 and n(L r + C r , −i) = 0 were ad hoc, and involved the manual construction of auxiliary paths or rays (e.g., C R and L R ). Similar difficulties have also been mentioned by John Harrison when formalising the prime number theorem [8]. In the next section, we will introduce an idea to systematically evaluate winding numbers.
. Right: the image of f as a point travels through γ

The Intuition
The fundamental idea of evaluating a winding number n(γ , z 0 ) in this paper is to reduce the evaluation to classifications of how paths cross the line {z | Re(z) = Re(z 0 )}: continuously or not and in which direction. In a simple case, suppose a path γ crosses the line {z | Re(z) = Re(z 0 )} exactly once at the point γ (t 0 ) such that Im(γ (t 0 )) > Im(z 0 ) (see Fig. 2 (left)), and let θ be the change in the argument of a complex point travelling through γ . It should not be hard to observe that 0 < θ < 2π, and by considering Re(n(γ , z 0 )) = θ/(2π) we can have 0 < Re(n(γ , z 0 )) < 1, which is an approximation of Re(n(γ , z 0 )). That is, we have approximated Re(n(γ , z 0 )) by the way that γ crosses the line {z | Re(z) = Re(z 0 )}.
To make this idea more precise, let The image of f as a point travels through γ is as illustrated in Fig. 2 (right), where f jumps from +∞ to −∞ across t 0 . We can then formally characterise those jumps.
Specifically, we can conjecture that jump captures the way that γ crosses the line {z | Re(z) = Re(z 0 )} in Fig. 2, hence Re(n(γ , z 0 )) can be approximated using jump + and jump − : In more general cases, we can define Cauchy indices by summing up these jumps over an interval and along a path.
Definition 3 (Cauchy index along a path) Given a path γ : [0, 1] → C and a point z 0 ∈ C, the Cauchy index along γ about z 0 is defined as In particular, it can be checked that the Cauchy index Indp(γ , z 0 ) captures the way that γ crosses the line {z | Re(z) = Re(z 0 )}, hence leads to an approximation of Re(n(γ , z 0 )): More interestingly, by further knowing that γ is a loop we can derive Re(n(γ , z 0 )) = n(γ , z 0 ) ∈ Z and Indp(γ , z 0 )/2 ∈ Z, following which we come to the core proposition of this paper: Proposition 1 Given a valid path γ : [0, 1] → C and a point z 0 ∈ C, such that γ is a loop and z 0 is not on the image of γ , we have That is, under some assumptions, we can evaluate a winding number through Cauchy indices! A formal proof of Proposition 1 will be introduced in Sect. 4.1. Here, given the statement of the proposition, we can have alternative proofs for n(L r +C r , i) = 1 and n(L r +C r , −i) = 0. Example 3 (Alternative proof of n(L r + C r , i) = 1) As L r + C r is a loop, applying Proposition 1 yields which reduces n(L r + C r , i) to the evaluations of Indp(L r , i) and Indp(C r , i). In this case, by definition we can easily decide Indp(L r , i) = −1 and Indp(C r , i) = −1 as illustrated in Fig. 3a. Hence, we have and conclude the proof. Evaluating n(L r + C r , i) and n(L r + C r , −i) through the way that the path L r + C r crosses the imaginary axis Example 4 (Alternative proof of n(L r + C r , −i) = 0) As shown in Fig. 3b, we can similarly have by which the proof is completed.
Compared to the previous proofs presented in Examples 1 and 2, the alternative proofs in Examples 3 and 4 are systematic and less demanding to devise once we have a formalisation of Proposition 1, which is what we will introduce in the next section.

Evaluating Winding Numbers
The previous section presented an informal intuition to systematically evaluate winding numbers; in this section, we will report the formal development of this intuition. We will first present a mechanised proof of Proposition 1 (Sect. 4.1), which includes mechanised definitions of jumps and Cauchy indices (i.e., Definition 1, 2 and 3) and several related properties of these objects. After that, we build a tactic in Isabelle/HOL that is used to mechanise proofs presented in Example 3 and 4 (Sect. 4.2). Finally, we discuss some subtleties we encountered during the formalisation (Sect. 4.3).

A Formal Proof of Proposition 1
For jump − and jump + (see Definition 1), we have used the filter mechanism [9] to define a function jumpF: respectively. Here, at_left x, at_right x, at_top, and at_bot are all filters, where a filter is a predicate on predicates that satisfies certain properties. Filters are extensively used in the analysis library of Isabelle to encode varieties of logical quantification: for example, at_left x encodes the statement "for a variable that is sufficiently close to x from the left", and at_top represents "for a sufficiently large variable". Furthermore, LIM x (at_left x). f x :> at_top encoded the proposition and this encoding can be justified by the following equality in Isabelle: where ∀ z. ∃ b<x. ∀ y>b. y < x −→ z ≤ f y matches the usual definition of (13) in textbooks. We can then encode Ind b a ( f ) and Indp(γ , z 0 ) (see Definitions 2 and 3) as cindexE and cindex_pathE respectively: which actually hides an assumption: that only a finite number of points within the interval [a, b) contribute to the sum. This assumption is made explicit when cindexE is defined by summing jumps over the following set: {x. jumpF f (at_right x) = 0 ∧ a ≤ x ∧ x < b}.
If the set above is infinite (i.e., the sum x∈ [a,b) jump + ( f , x) is not mathematically welldefined) we have In other words, Isabelle/HOL deems the sum over an infinite set to denote zero. Due to the issue of well-defined sums, many of our lemmas related to cindexE should assume a finite number of jumps: which guarantees the well-definedness of cindexE. Now, suppose that we know that Indp is well-defined: there are only a finite number of jumps over the path. What strategy can we employ to formally prove Proposition 1? Naturally, we may want to divide the path into a finite number of segments (subpaths) separated by those jumps, and then perform inductions on these segments. To formalise the finiteness of such segments, we defined an inductive predicate: finite_Psegments P a s]] ⇒ finite_Psegments P a b" definition finite_ReZ_segments:: The idea behind finite_ReZ_segments is that a jump of takes place only if λt. Re(γ (t) − z 0 ) changes from 0 to = 0 (or vice versa). Hence, each of the segments of the path γ separated by those jumps has either λt.
As can be expected, the finiteness of jumps over a path can be derived by the finiteness of segments:
By assuming such a finite number of segments we have well-defined cindex_pathE, and can then derive some useful related properties:

Lemma 4 (cindex_pathE_subpath_combine)
fixes γ ::"real ⇒ complex" and z 0 ::complex assumes "finite_ReZ_segments γ z 0 "and "path γ " and "0≤a" and "a≤b" and "b≤c" and "c≤1" where subpath a b γ gives a sub-path of γ based on parameters a and b: Essentially, Lemma 4 indicates that we can combine Cauchy indices along consecutive parts of a path: given a path γ and three parameters a, b, More importantly, we now have an induction rule for a path with a finite number of segments: Lemma 5 (finite_ReZ_segments_induct) fixes γ ::"real ⇒ complex" and z 0 ::complex and P::"(real ⇒ complex) ⇒ complex ⇒ bool" assumes "finite_ReZ_segments γ z 0 " and sub0:" g z. (P (subpath 0 0 g) z)" and subEq: where P is a predicate that takes a path γ and a complex point z 0 , and -sub0 is the base case that P holds for a constant path; -subEq is the inductive case when the last segment is right on the line {x | Re(x) = Re(z)}: ∀t ∈ (s, 1). Re(g(t)) = Re(z); -subNEq is the inductive case when the last segment does not cross the line {x | Re(x) = Re(z)}: ∀t ∈ (s, 1). Re(g(t)) = Re(z).
Given a path γ with a finite number of segments, a complex point z 0 and a predicate P that takes a path and a complex number and returns a boolean, Lemma 5 provides us with an inductive rule to derive P(γ , z 0 ) by recursively examining the last segment. Before attacking Proposition 1, we can show an auxiliary lemma about Re(n(γ , z 0 )) and Indp(γ , z 0 ) when the end points of γ are on the line {z | Re(z) = Re(z 0 )}: Lemma 6 (winding_number_cindex_pathE_aux) fixes γ ::"real ⇒ complex" and z 0 :: complex assumes "finite_ReZ_segments γ z 0 " and "valid_path γ " and "z 0 / ∈ path_image γ " and "Re (γ 1) = Re z 0 " and "Re (γ 0) = Re z 0 " shows "2 * Re(winding_number γ z 0 ) = -cindex_pathE γ z 0 " Here, Lemma 6 is almost equivalent to Proposition 1 except for that more restrictions haven been placed on the end points of γ .

Proof of Lemma 6
As there are a finite number of segments along γ (i.e., finite_ReZ_ segments γ z 0 ), by inducting on these segments with Lemma 5 we end up with three cases.
For the inductive case when the last segment is right on the line We have and, by the induction hypothesis, 2 Re(n(g 1 , z)) = − Indp(g 1 , z).
For the other inductive case when the last segment does not cross the line {x | Re(x) = Re(z)}, without loss of generality, we assume and the shape of g is as illustrated in Fig. 4b. Similar to the previous case, by letting g 1 (t) = g(st) and g 2 (t) = g((1 − s)t), we have n(g, z) = n(g 1 , z) + n(g 2 , z) and, by the induction hypothesis, 2 Re(n(g 1 , z)) = − Indp(g 1 , z). Moreover, by observing the shape of g 2 we have (20) with (21) leads to 2 Re(n(g 2 , z)) = −Indp(g 2 , z), following which we finish the case by deriving 2 Re(n(g, z)) = − Indp(g, z) in a way analogous to (18).
In the case of failure, without loss of generality, we can assume Re(γ (t)) > Re(z 0 ) for all 0 ≤ t ≤ 1, and the shape of γ is as illustrated in Fig. 5a. As the path γ does not cross the line {z | Re(z) = Re(z 0 )}, we can evaluate where Ln is the principle value of a complex logarithm function with its branch being the negative real axis and −π < Im(Ln(z)) ≤ π for all z. Hence, n(γ , z 0 ) = − Indp(γ , z 0 )/2 which concludes the case. In the case of success, as illustrated in Fig. 5b, we have Re(γ (s)) = Re(z 0 ). We then define a shifted path γ s : such that Re(γ s (0)) = Re(γ s (1)) = Re(z 0 ). By applying Lemma 6, we obtain a relationship between Re(n(γ s , z 0 )) and Indp(γ s , z 0 ):

A Tactic for Evaluating Winding Numbers
With Proposition 1 formalised, we are now able to build a tactic to evaluate winding numbers using Cauchy indices. The idea has already been sketched in Examples 3 and 4. We have built a tactic eval_winding, for goals of the form where k is an integer and γ j (1 ≤ j ≤ n) is either a linear path: or a part of a circular path: The tactic eval_winding will transform (22) into where (23) ensures that the path γ 1 + γ 2 + · · · + γ n is a loop; (24) certifies that z 0 is not on the image of γ 1 + γ 2 + · · · + γ n . To achieve this transformation, eval_winding will first perform a substitution step on the left-hand side of Eq. (22) using Theorem 1. As the substitution is conditional, we will need to resolve four extra subgoals (i.e., (26), (27), (28) and (29) as follows) and Eq. (22) is transformed into (30): -cindex_pathE (γ 1 +++ γ 2 +++ ... +++γ n ) z 0 / 2 = k.

Lemma 7 (finite_ReZ_segments_joinpaths)
fixes γ 1 γ 2 :: "real ⇒ complex" and z 0 :: complex assumes "finite_ReZ_segments γ 1 z 0 " and "finite_ReZ_segments γ 2 z 0 " and "path γ 1 " and "path γ 2 " and "γ 1 1 = γ 2 0" shows "finite_ReZ_segments (γ 1 +++γ 2 ) z 0 " to eliminate the path join operations (+++) until the predicate finite_ReZ_segments is only applied to a linear path or a part of a circular path, and either of these two cases can be directly discharged because these two kinds of paths are proved to be divisible into a finite number of segments by the imaginary axis:

Lemma 9 (finite_ReZ_segments_part_circlepath)
"finite_ReZ_segments (part_circlepath z0 r st tt) z" In terms of other subgoals introduced when applying Lemma 7, such as path γ 1 , path γ 2 and γ 1 1 = γ 2 0, we can discharge them by the following introduction and simplification rules (all of which have been formally proved): As a result, eval_winding will eventually simplify the subgoal (26) to (23).
Finally, with respect to (30), we can similarly rewrite with a rule between the Cauchy index (cindex_pathE) and the path join operation (+++):
After building the tactic eval_winding, we are now able to convert a goal like Eq. (22) to (23), (24) and (25). In most cases, discharging (23) and (24) is straightforward. To derive (25), we will need to formally evaluate each Indp(γ j , z 0 ) (1 ≤ j ≤ n) when γ j is either a linear path or a part of a circular path.
When γ j is a part of a circular path, a similar lemma has been provided to facilitate the evaluation of Indp(γ j , z 0 ).

Subtleties
The first subtlety we have encountered during the formalisation of Proposition 1 is about the definitions of jumps and Cauchy indices, for which our first attempt followed the standard definitions in textbooks [2,16,19].
The impact of the difference between the current definition of the Cauchy index (i.e., Definition 2) and the classic one (i.e., Definition 5) is small when formalising the Sturm-Tarski theorem [10,13], where f is a rational function. In this case, the path γ intersects with the line {z | Re(z) = Re(z 0 )} a finite number of times, and for each intersection point (see Fig. 6a, b), by letting provided jump + ( f , a) = 0 and jump − ( f , b) = 0. That is, the classic Cauchy index and the current one are equal when f is a rational function and does not jump at both ends of the target interval.
Naturally, the disadvantages of Definition 5 are twofold: -The function λt. Re(γ (t)−z 0 ) cannot vanish at either end of the interval. That is, we need to additionally assume Re(γ (0) − z 0 ) = 0 as in Rahman and Schmeisser's formulation [19, Lemma 11.1.1 and Theorem 11.1.3], and Proposition 1 will be inapplicable in the case of Fig. 6c where Re(γ (0)) = Re(γ (1)) = Re(z 0 ). -The function λt. Im(γ (t) − z 0 )/Re(γ (t) − z 0 ) has to be rational, which makes Proposition 1 inapplicable for cases like in Fig. 6d (if we follow Definition 5). To elaborate, it can be observed in Fig. 6d that n(γ , z 0 ) = −1, while we will only get a wrong answer by following Definition 5 and evaluating via Proposition 1: Fortunately, Eisermann [6] recently proposed a new formulation of the Cauchy index that overcomes those two disadvantages, and this new formulation is what we have followed (in Definitions 1 and 2).
Another subtlety we ran into was the well-definedness of the Cauchy index. Such welldefinedness is usually not an issue and left implicit in the literature, because, in most cases, the Cauchy index is only defined on rational functions, where only finitely many points can contribute to the sum. When attempting to formally derive Proposition 1, we realised that this assumption needed to be made explicit, since the path γ can be flexible enough to allow the function f (t) = Im(γ (t) − z 0 )/Re(γ (t) − z 0 ) to be non-rational (e.g., Fig. 6d). In our first attempt of following Definition 5, the Cauchy index was formally defined as follows: and its well-definedness was ensured by the finite number of times that γ crosses the line {z | Re(z) = Re(z 0 )}: where the part Re (γ t -z 0 ) = 0 ensures that jump f t is non-zero only at finitely many points over the interval [0, 1]. When constrained by finite_axes_cross, the function behaves like a rational function. More importantly, the path γ , in this case, can be divided into a finite number of ordered segments delimited by those points over [0, 1], which makes an inductive proof of Proposition 1 possible. However, after abandoning our first attempt and switching to Definition 2, the well-definedness of the Cauchy index is assured by the finite number of jump + and jump − of f (i.e., Definition finite_jumpFs in Sect. 4.1), with which we did not know how to divide the path γ into segments and carry out an inductive proof. It took us some time to properly define the assumption of a finite number of segments (i.e., Definition finite_ReZ_segments) that implied the well-definedness using Lemma 3 and provided a lemma for inductive proofs (i.e., Lemma 5).

Counting the Number of Complex Roots
The previous section described a way to evaluate winding numbers via Cauchy indices. In this section, we will further explore this idea and propose verified procedures to count the number of complex roots of a polynomial in some domain such as a rectangle and a half-plane. Does a winding number have anything to do with the number of roots of a polynomial? The answer is yes. Thanks to the argument principle, we can calculate the number of roots by evaluating a contour integral: where p ∈ C[x], p (x) is the first derivative of p and N is the number of complex roots of p (counted with multiplicity) inside the loop γ . Also, by the definition of winding numbers, we have Combining Eqs. (31) and (32) gives us the relationship between a winding number and the number of roots of a polynomial: And the question becomes: can we evaluate n( p • γ, 0) via Cauchy indices?

Roots in a Rectangle
Let N be the number of complex roots of a polynomial p inside the rectangle defined by its lower left corner a 1 and upper right corner a 3 . As illustrated in Fig. 7, we can define four linear paths along the edge of the rectangle: where a 2 = Re(a 3 ) + i Im(a 1 ) and a 4 = Re(a 1 ) + i Im(a 3 ). Combining Proposition 1 with Eq. (33) yields Here, the path p • L j : [0, 1] → C (1 ≤ j ≤ 4) is (mostly) neither a linear path nor a part of a circular path, which indicates that the evaluation strategies of Sect. 4.2, such as Lemma 11, will no longer apply. Thankfully, the Sturm-Tarski theorem [10,13] came to our rescue. In general, the Sturm-Tarski theorem is about calculating Tarski queries through sign variations and signed remainder sequences: let p, q ∈ R[x], a and b be two extended real numbers such that a < b and are not roots of p, we have TaQ(q, p, a, b) = Var (SRemS( p, p q) where p is the first derivative of p, -the Tarski query TaQ(q, p, a, b) defined as follows: sgn(q(x)), -SRemS( p, q) is the signed remainder sequence started with p and q.
Note that when q = 1, (35) becomes the famous Sturm's theorem, which counts the number of distinct real roots over an interval. For example, by calculating we know that the polynomial x 2 − 3x + 2 has two distinct real roots within the interval (0, 3). In our previous formal proof of the Sturm-Tarski theorem [10,13], we used the Cauchy index to relate the Tarski query and the right-hand side of (35). Therefore, as a byproduct, we can also evaluate the Cauchy index through sign variations and signed remainder sequences: where p, q ∈ R[x], a, b are two extended real numbers such that a < b and are not roots of p.
Back to the case of Indp( p • L j , 0), we have and both Im( p(L j (t))) and Re( p(L j (t))) happen to be polynomials with real coefficients. Therefore, combining Eqs. (34) and (37) yields an approach to count the number of roots inside a rectangle. While proceeding to the formal development, the first problem we encountered was that the Cauchy index in Eq. (37) actually follows the classic definition (i.e., Definition 5), and is different from the one in Eq. (34) (i.e., Definitions 2 and 3). Subtle differences between these two formulations have already been discussed in Sect. 4.3. Luckily, Eisermann [6] has also described an alternative sign variation operator so that our current definition of the Cauchy index (i.e., Definition 2) can be computationally evaluated: Lemma 12 (cindex_polyE_changes_alt_itv_mods) fixes a b::real and p q::"real poly" assumes "a < b" and "coprime p q" shows "cindex_polyE a b q p = changes_alt_itv_smods a b p q / 2" Here, cindex_polyE a b q p encodes our current definition of the Cauchy index Ind b a (λt. q(t)/ p(t)), and changes_alt_itv_smods a b p q stands for Var (SRemS( p, q); a, b) where the alternative sign variation operator Var is defined as follows: Var ([ p 1 , p 2  Before implementing Eq. (34), we need to realise that there is a restriction in our strategy: roots are not allowed on the border (i.e., the image of the path L 1 + L 2 + L 3 + L 4 ). To computationally check this restriction, the following function is defined definition no_proots_line::"complex poly ⇒ complex ⇒ complex ⇒ bool" where "no_proots_line p a b = (proots_within p (closed_segment a b) = {})" which will return "true" if there is no root on the closed segment between a and b, and "false" otherwise. Here, closed_segment a b is defined as the set {(1−u)a+ub | 0 ≤ u ≤ 1} ⊆ C, and the function proots_within p s gives the set of roots of the polynomial p within the set s: definition proots_within::"'a::comm_semiring_0 poly ⇒ 'a set ⇒ 'a set" where "proots_within p s = {x∈s. poly p x=0}" The next step is to make the definition no_proots_line executable. This is achieved by proving a code equation, where the left-hand side of the equation is the target definition and the right-hand side is an executable expression. In the case of no_proots_line, the code equation is the following lemma:

Lemma 13 (no_proots_line_code[code])
"no_proots_line p a b = (if poly p a = 0 ∧ poly p b = 0 then (let p c = p • p [:a, b -a:]; p R = map_poly Re p c ; p I = map_poly Im p c ; g = gcd p R p I in if changes_itv_smods 0 1 g (pderiv g) = 0 then True else False) else False)" where • p is the polynomial composition operation and map_poly Re and map_poly Im, respectively, extract the real and imaginary parts of the complex polynomial p c .

Proof of Lemma 13
Supposing L : [0, 1] → C is a linear path from a to b: L(t) = (1 − t)a + tb, we know that p • L is still a polynomial with complex coefficients. Subsequently, we extract the real and imaginary parts ( p R and p I , respectively) of p • L such that

p(L(t)) = p R (t) + i p I (t).
If there is a root of p lying right on L, we will be able to obtain some t 0 ∈ [0, 1] such that hence, by letting g = gcd( p R , p I ) we have g(t 0 ) = 0. Therefore, the polynomial p has no (complex) root on L if and only if g has no (real) root within the interval [0, 1], and the latter can be computationally checked using Sturm's theorem.
Finally, we define the function proots_rectangle that returns the number of complex roots of a polynomial (counted with multiplicity) within a rectangle defined by its lower left and upper right corner: where proots_count p s denotes the number of roots of the polynomial p within the set s: definition proots_count::"'a::idom poly ⇒ 'a set ⇒ nat" where "proots_count p s = ( r∈proots_within p s. order r p)" The executability of the function proots_rectangle can be established with the following code equation:
Example 5 Given a rectangle defined by (−1, 2+2i) (as illustrated in Fig. 8) and a polynomial p with complex coefficients: we can now type the following command to count the number of roots within the rectangle: value "proots_rectangle [:-1, -2*i, 1:] (-i) (2+2*i)" which will return 2 as p has exactly two complex roots (i.e., i with multiplicity 2) in the area.

Roots in a Half-plane
For roots in a half-plane, we can start with a simplified case, where we count the number of roots of a polynomial in the upper half-plane of C: definition proots_upper::"complex poly ⇒ int" where "proots_upper p = proots_count p {z. Im z>0}" As usual, our next step is to set up the executability of proots_upper. To achieve that, we first define a linear path L r (t) = (1 − t)(−r ) + tr and a semi-circular path C r (t) = re iπ t , as illustrated in Fig. 9. Subsequently, let and by following Eq. (33) we have where N r is the number of roots of p inside the path L r + C r . Note that as r approaches positive infinity, N r will be the roots on the upper half-plane (i.e., proots_upper p), which is what we are aiming for. For this reason, it is natural for us to examine two cases: provided that the polynomial p is monic and does not have any root on the real axis. Next, for lim r →+∞ Re(n(C p (r ), 0)), we first derive a lemma about C r :
Finally, following Lemma 18, the executability of the function proots_upper is established:

Lemma 19 (proots_upper_code1[code])
"proots_upper p = (if p = 0 then (let p m = smult (inverse (lead_coeff p)) p; p I = map_poly Im p m ; p R = map_poly Re p m ; g = gcd p I p R in if changes_R_smods g (pderiv g) = 0 then (degree p -changes_R_smods p R p I ) div 2 else Code.abort (STR "proots_upper fails when there is a root on the border.") (λ_. proots_upper p) ) else Code.abort (STR "proots_upper fails when p=0.") (λ_. proots_upper p))" where smult (inverse (lead_coeff p)) p divides the polynomial p by its leading coefficient so that the resulting polynomial p m is monic. This corresponds to the assumption lead_coeff p=1 in Lemma 18. -changes_R_smods g (pderiv g) = 0 checks if p has no root lying on the real axis, which is due to the second assumption in Lemma 18. -changes_R_smods p R p I evaluates Ind +∞ −∞ λt.
As for the general case of a half-plane, we can have a definition as follows: has exactly two roots within the left half-plane of the vector (0, i), as shown in Fig. 10.
Despite our naive implementation, both proofs_half and proots_rectangle are applicable for small or medium examples. For most polynomials with coefficient bitsize up to 10 and degree up to 30, our complex root counting procedures terminate within minutes.

Limitations and Future Work
There are, of course, several improvements that can be made on both the evaluation tactic of Sect. 4.2 and the root counting procedures of Sect. 5. As the tactic is intended to be applied to winding numbers with variables, full automation with this tactic is unlikely in most cases, but we can always aim for better automation and an enhanced interactive experience for users (e.g., presenting unsolved goals in a more user-friendly way).
Regarding the two root-counting procedures in Sect. 5, a key limitation is that they do not allow cases where any root is on the border. There are two possible solutions to this problem: -To generalise the definition of winding numbers. The current formulation of winding numbers in Isabelle/HOL follows the one in complex analysis: which becomes undefined when the point z is on the image of the path γ . With more general formulations of winding numbers, such as the algebraic version by Eisermann [6], we may be able to derive a stronger version of the argument principle that allows zeros on the border.
-To deploy a more sophisticated strategy to count the number of times that the path winds.
Recall that the underlying idea in this paper is to reduce the evaluation of winding numbers to classifications of how paths cross some line. The Cauchy index merely provides one classification strategy, which we considered simple and elegant enough for formalisation. In contrast, Collins and Krandick [4] propose a much more sophisticated strategy for such classifications. Their strategy has, in fact, been widely implemented in modern systems, such as Mathematica and SymPy, to count the number of complex roots.
Neither of these two solutions are straightforward to incorporate, hence we leave them for future investigation. Besides rectangles and half-planes, it is also possible to similarly count the number of roots in an open disk and even a sector: where arg(−) returns the argument of a complex number. Informal proofs of root counting in these domains can be found in Rahman and Schmeisser [19,Chapter 11].

Potential Applications
Rahman and Schmeisser's book [19,Chapter 11] and Eisermann's paper [6] are the two main sources that our development is built upon. Nevertheless, there are still some differences in formulations: -Rahman and Schmeisser formulated the Cauchy index as in Definitions 4 and 5, and we used their formulation in our first attempt. However, after we realised the subtleties discussed in Sect. 4.3, we abandoned this formulation and switched to Eisermann's (i.e., Definition 2). As a result, the root counting procedures presented in this paper are more general than the ones in their book, having fewer preconditions. -Eisermann formulated a winding number n(γ , z 0 ) in a real-algebraical sense where γ is required to be a piecewise polynomial path (i.e., each piece from the path needs to be a polynomial). In comparison, n(γ , z 0 ) in Isabelle/HOL follows the classic definition in complex analysis, and places fewer restrictions on the shape of γ (i.e., piecewise continuously differentiable is less restrictive than being a piecewise polynomial) but does not permit z 0 to be on the image of γ (while Eisermann's formulation does). Consequently, Eisermann's root counting procedure works in more restrictive domains (i.e., he only described the rectangle case in his paper) but does not prevent roots on the border.
Another point worth mentioning is the difference between informal and formal proofs. In this development, we generally treated their lemma statements as bald facts: we had to discover our own proofs. For instance, when proving Proposition 1, we defined an inductive data type for segments and derived an induction rule for it, which was nothing like the informal proof. Such situations also happened when we justified the root counting procedure in a half-plane. Overall, the formal proofs are about 12,000 lines. Interestingly, the root-counting procedure in a half-plane is also related to the stability problems in the theory of dynamical systems. For instance, let A ∈ R n×n be a square matrix with real coefficients and y : [0, +∞) → R n be a function that models the system state over time. A linear dynamical system can be described as an ordinary differential equation: with an initial condition y(0) = y 0 . The system of (42) is considered stable if all roots of the characteristic polynomial of A lie within the open left half-plane (i.e., {z | Re(z) < 0}), and this stability test is usually referred as the Routh-Hurwitz stability criterion [1, Section 23], [16,Chapter 9]. As has been demonstrated in Example 6, counting the number of roots in the left half-plane is within the scope of the procedure proots_half. For this reason, we believe that the development in this paper will be beneficial for reasoning about dynamical systems in Isabelle/HOL. It is worth mentioning that root counting in a rectangle is usually coupled with a classic problem in computer algebra, namely, complex root isolation. The basic idea is to keep bisecting a rectangle (vertically or horizontally) into smaller ones until a sub-rectangle contains exactly one root or none (provided the target polynomial is square-free). Following this idea, it is possible to build a simple and verified procedure for complex root isolation similar to Wilf's algorithm [20]: we start with a large rectangle and then repeatedly apply the verified procedure to count roots during the rectangle bisection phase. However, compared to modern complex procedures [4,21], this simplistic approach suffers from several drawbacks: -Our root counting procedure is based on remainder sequences, which are generally considered much slower than those built upon Descartes' rule of signs. -Modern isolation procedures are routinely required to deliver isolation boxes whose sizes meet some user-specified limit, hence they usually keep refining the isolation boxes even after the roots have been successfully isolated. The bisection strategy still works in the root refinement stage, but dedicated numerical approaches such as Newton's iteration are commonly implemented for efficiency reasons. -Modern isolation procedures sometimes prefer a bit-stream model in which the coefficients of the polynomial are approximated as a bit stream. This approach is particularly beneficial when the coefficients have extremely large bit-width or consist of algebraic numbers. -Modern implementations usually incorporate numerous low-level optimisations, such as hash tables, which are hard to implement as verified procedures in a theorem prover.
Therefore, it is unlikely that our verified root counting procedures will ever deliver high performance. Nevertheless, they can be used to certify results from untrusted external root isolation programs, as in the certificate-based approach to solving univariate polynomial problems [13].

Related Work
Formalisations of the winding number (from an analytical perspective) are available in Coq [3], HOL Light [7] and Isabelle/HOL. To the best of our knowledge, our tactic of evaluating winding numbers through Cauchy indices is novel. As both HOL Light and Isabelle/HOL have a relatively comprehensive library of complex analysis (i.e., at least including Cauchy's integral theorem), our evaluation tactic could be useful when deriving analytical proofs in these two proof assistants. The ability to count the real roots of a polynomial only requires Sturm's theorem, so this capability is widely available among major proof assistants including PVS [18], Coq [15], HOL Light [17] and Isabelle [5,10,13]. However, as far as we know, our procedures to count complex roots are novel, as they require a formalisation of the argument principle [14], which is only available in Isabelle at the time of writing.
This can be then solved by individually evaluating Indp(γ 1 , z 0 ), . . . , Indp(γ n , z 0 ). As open variables may occur in those Cauchy indices, the evaluation of them is unlikely to be fully automatic, but we provide lemmas (e.g., Lemma 11) to mitigate the laborious process. The tactic winding_eval has greatly helped us with the motivating proofs shown in Sect. 2, and we believe that it should be also beneficial in similar situations when dealing with winding numbers in a formal framework.
We have further related Cauchy indices to the argument principle and developed novel verified procedures to count the complex roots of a polynomial within the areas of rectangles and half-planes. Despite the limitations of not allowing roots on the border (which we will solve in future work), the ability to formally count complex roots is believed to lay the foundations for conducting stability analysis (e.g., the Routh-Hurwitz stability criterion) in the framework of the Isabelle theorem prover. duction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.