David Monniaux

Pierre Corbineau

On the Generation of Positivstellensatz Witnesses in Degenerate Cases *

One can reduce the problem of proving that a polynomial is nonnegative, or more generally of proving that a system of polynomial inequalities has no solutions, to finding polynomials that are sums of squares of polynomials and satisfy some linear equality (Positivstellensatz ). This produces a witness for the desired property, from which it is reasonably easy to obtain a formal proof of the property suitable for a proof assistant such as Coq.

The problem of finding a witness reduces to a feasibility problem in semidefinite programming, for which there exist numerical solvers. Unfortunately, this problem is in general not strictly feasible, meaning the solution can be a convex set with empty interior, in which case the numerical optimization method fails. Previously published methods thus assumed strict feasibility; we propose a workaround for this difficulty.

We implemented our method and illustrate its use with examples, including extractions of proofs to Coq.

Introduction

Consider the following problem: given a conjunction of polynomial equalities, and (wide and strict) polynomial inequalities, with integer or rational coefficients, decide whether this conjunction is satisfiable over R; that is, whether one can assign real values to the variables so that the conjunction holds. A particular case is showing that a given polynomial is nonnegative.

The decision problem for real polynomial inequalities can be reduced to quantifier elimination: given a formula F , whose atomic formulas are polynomial (in)equalities, containing quantifiers, provide another, equivalent, formula F ′ , whose atomic formulas are still polynomial (in)equalities, containing no quantifier. An algorithm for quantifier elimination over the theory of real closed fields (roughly speaking, (R, 0, 1, +, ×, ≤) was first proposed by Tarski [START_REF] Tarski | A Decision Method for Elementary Algebra and Geometry[END_REF][START_REF] Seidenberg | A New Decision Method for Elementary Algebra[END_REF], but this algorithm had non-elementary complexity and thus was impractical. Later, the cylindrical algebraic decomposition (CAD) algorithm was proposed [START_REF] Collins | Quantifier elimination for real closed fields by cylindrical algebraic decomposition[END_REF], with a doubly exponential complexity, but despite improvements [START_REF] Collins | Quantifier Elimination by Cylindrical Algebraic Decomposition -Twenty Years of Progress[END_REF] CAD is still slow in practice and there are few implementations available.

Quantifier elimination is not the only decision method. Basu et al. [2, Theorem 3] proposed a satisfiability testing algorithm with complexity s k+1 d O(k) , where s is the number of distinct polynomials appearing in the formula, d is their maximal degree, and k is the number of variables. We know of no implementation of that algorithm. Tiwari [START_REF] Tiwari | An Algebraic Approach for the Unsatisfiability of Nonlinear Constraints[END_REF] proposed an algorithm based on rewriting systems that is supposed to answer in reasonable time when a conjunction of polynomial inequalities has no solution.

Many of the algebraic algorithms are complex, which leads to complex implementations. This poses a methodology problem: can one trust their results? The use of computer programs for proving lemmas used in mathematical theorems was criticized in the case of Thomas Hales' proof of the Kepler conjecture. Similarly, the use of complex decision procedures (as in the proof assistant PVS 1 ) or program analyzers (as, for instance, Astrée2 ) in order to prove the correctness of critical computer programs is criticized on grounds that these verification systems could themselves contain bugs.

One could formally prove correct the implementation of the decision procedure using a proof assistant such as Coq [START_REF][END_REF][START_REF] Mahboubi | Implementing the CAD algorithm inside the Coq system[END_REF]; but this is likely to be long and difficult. An alternative is to arrange for the procedure to provide a witness of its result. The answer of the procedure is correct if the witness is correct, and correctness of the witness can be checked by a much simpler procedure, which can be proved correct much more easily.

Unsatisfiability witnesses for systems of complex equalities or linear rational inequalities are already used within DPLL(T ) satisfiability modulo theory decision procedures [16, ch. 11] [START_REF] Dutertre | Integrating Simplex with DPLL(T)[END_REF]. It is therefore tempting to seek unsatisfiability witnesses for systems of polynomial inequalities.

In recent years, it was suggested [START_REF] Parrilo | Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization[END_REF] to use numerical semidefinite programming to look for proof witnesses whose existence is guaranteed by a Positivstellensatz [START_REF] Krivine | Anneaux préordonnés[END_REF][START_REF] Stengle | A Nullstellensatz and a Positivstellensatz in Semialgebraic Geometry[END_REF][START_REF] Schweighofer | On the complexity of Schmüdgen's Positivstellensatz[END_REF]. The original problem of proving that a system of polynomial inequalities has no solution is reduced to: given polynomials P i and R, derived from those in the original inequalities, find polynomials Q i that are sums of squares such that i P i Q i = R. Assuming some bounds on the degrees of Q i , this problem is in turn reduced to a semidefinite programming pure feasibility problem [START_REF] Boyd | Convex Optimization[END_REF][START_REF] Vandenberghe | Semidefinite Programming[END_REF], a form of convex optimization. The polynomials Q i then form a witness, from which a machine-checkable formal proof, suitable for tools such as Coq [START_REF][END_REF] or Isabelle [START_REF] Harrison | Verifying nonlinear real formulas via sums of squares[END_REF], may be constructed.

Unfortunately, this method suffers from a caveat: it applies only under a strict feasibility condition [START_REF] Peyrl | Computing sum of squares decompositions with rational coefficients[END_REF]: a certain convex geometrical object should not be degenerate, that is, it should have nonempty interior. Unfortunately it is very easy to obtain problems where this condition is not true. Equivalently, the method of rationalization of certificates [START_REF] Kaltofen | Exact certification of global optimality of approximate factorizations via rationalizing sums-of-squares with floating point scalars[END_REF] has a limiting requirement that the rationalized moment matrix remains positive semidefinite.

In this article, we explain how to work around the degeneracy problem: we propose a method to look for rational solutions to a general SDP feasibility problem. We have implemented our method and applied it to some examples from the literature on positive polynomials, and to examples that previously published techniques failed to process.

Witnesses

For many interesting theories, it is trivial to check that a given valuation of the variables satisfies a quantifier-free formula. A satisfiability decision procedure will in this case tend to seek a satisfiability witness and provide it to the user when giving a positive answer.

In contrast, if the answer is that the problem is not satisfiable, the user has to trust the output of the satisfiability testing algorithm, the informal meaning of which is "I looked carefully everywhere and did not find a solution." In some cases, it is possible to provide unsatisfiability witnesses: solutions to some form of dual or auxiliary problem that show that the original problem had no solution.

Nonnegativity Witnesses

To prove that a polynomial P is nonnegative, one simple method is to express it as a sum of squares of polynomials. One good point is that the degree of the polynomials involved in this sum of squares can be bounded, and even that the choice of possible monomials is constrained by the Newton polytope of P , as seen in §3.

Yet, there exist nonnegative polynomials that cannot be expressed as sums of squares, for instance this example due to Motzkin [START_REF] Reznick | Some Concrete Aspects Of Hilbert's 17th Problem[END_REF]:

M = x 6 1 + x 4 2 x 2 3 + x 2 2 x 4 3 -3x 2 1 x 2 2 x 2 3 (1)
However, Artin's answer to Hilbert's seventeenth problem is that any nonnegative polynomial can be expressed as a sum of squares of rational functions. 3It follows that such a polynomial can always be expressed as the quotient Q 2 /Q 1 of two sums of squares of polynomials, which forms the nonnegativity witness, and can be obtained by solving P.Q 1 -Q 2 = 0 for Q 1 = 0 (this result is also a corollary of Th. 1).

Unsatisfiability Witnesses for Polynomial Inequalities

For the sake of simplicity, we shall restrict ourselves to wide inequalities (the extension to mixed wide/strict inequalities is possible). Let us first remark that the problem of testing whether a set of wide inequalities with coefficients in a subfield K of the real numbers is satisfiable over the real numbers is equivalent to the problem of testing whether a set of equalities with coefficients K is satisfiable over the real numbers: for each inequality P (x 1 , . . . , x m ) ≥ 0, replace it by P (x 1 , . . . , x m )µ 2 = 0, where µ is a new variable. Strict inequalities can also be simulated as follows: P i (x 1 , . . . , x m ) = 0 is replaced by P i (x 1 , . . . , x m ).µ = 1 where µ is a new variable. One therefore does not gain theoretical simplicity by restricting oneself to equalities.

Stengle [START_REF] Stengle | A Nullstellensatz and a Positivstellensatz in Semialgebraic Geometry[END_REF] proved two theorems regarding the solution sets of systems of polynomial equalities and inequalities over the reals (or, more generally, over real closed fields): a Nullstellensatz and a Positivstellensatz ; a similar result was proved by Krivine [START_REF] Krivine | Anneaux préordonnés[END_REF]. Without going into overly complex notations, let us state consequences of these theorems.

Let K be an ordered field (such as Q) and K ′ be a real closed field containing K (such as the real field R), and let X be a list of variables X 1 , . . . , X n . A * 2 denotes the squares of elements of A. The multiplicative monoid generated by A is the set of products of zero of more elements from A. The ideal generated by A is the set of sums of products of the form P Q where Q ∈ K[X] and P ∈ A. The positive cone generated by A is the set of sums of products of the form p.P.Q 2 where p ∈ K, p > 0, P is in the multiplicative monoid generated by A, and Q ∈ K[X]. Remark that we can restrict P to be in the set of products of elements of A where no element is taken twice, with no loss of generality.

The result [START_REF] Lombardi | Une borne sur les degrés pour le Théorème des zéros réel effectif[END_REF][START_REF] Lombardi | Théorème des zéros réels effectif et variantes[END_REF][START_REF] Coste | Dynamical method in algebra: effective Nullstellensätze[END_REF] of interest to us is:

Theorem 1. Let F > , F ≥ , F = , F = be sets of polynomials in K[X]
, to which we impose respective sign conditions > 0, ≥ 0, = 0, = 0. The resulting system is unsatisfiable over K ′ n if and only if there exist an equality in K[X] of the type S + P + Z = 0, with S in the multiplicative monoid generated by F > ∪ F * 2 = , P belongs to the positive cone generated by F ≥ ∪ F > , and Z belongs to the ideal generated by F = .

(S, P, Z) then constitute a witness of the unsatisfiability of the system. 4For a simple example, consider the following system, which obviously has no solution:

-2

+ y 2 ≥ 0 1 -y 4 ≥ 0 (2) A Positivstellensatz witness is y 2 (-2 + y 2 ) + 1(1 -y 4 ) + 2y 2 + 1 = 0. Another is 2 3 + y 2 3 (-2 + y 2 ) + 1 3 (1 -y 4 ) + 1 = 0. Consider the conjunction C: P 1 ≥ 0∧• • •∧P n ≥ 0 where P i ∈ Q[X 1 , . . . , X m ].
Consider the set Π of products of the form i P wi i for w ∈ {0, 1} n -that is, the set of all products of the P i where each P i appears at most once. Obviously, if one can exhibit nonnegative functions Q j such that Tj ∈Π Q j T j + 1 = 0, then C does not have solutions. Theorem 1 guarantees that if C has no solutions, then such functions Q j exist as sum of squares of polynomials (we simply apply the theorem with F > = F = = ∅ and thus S = {1}). We have again reduced our problem to the following problem: given polynomials T j and R, find sums-ofsquares polynomials Q j such that j Q j T j = R. Because of the high cost of enumerating all products of the form i P wi i , we have first looked for witnesses of the form Tj ∈S Q j P j + 1 = 0.

Solving the Sums-of-Squares Problem

In §2.1 and §2.2, we have reduced our problems to: given polynomials (P j ) 1≤j≤n and R in Q[X 1 , . . . , X m ], find polynomials that are sums of squares Q j such that j

P j Q j = R (3)
We wish to output the

Q j as Q j = nj i=1 α ji L 2
ji where α ji ∈ Q + and L ji are polynomials over Q. We now show how to solve this equation.

Reduction to Semidefinite Programming

Lemma 1. Let P ∈ K[X, Y, . . . ] be a sum of squares of polynomials i P 2

i . Let M = {m 1 , . . . , m |M| } be a set such that each P i can be written as a linear combination of elements of M (M can be for instance the set of monomials in the P i ). Then there exists a |M | × |M | symmetric positive semidefinite matrix Q with coefficients in K such that P (X, Y, . . .

) = [m 1 , . . . , m |M| ]Q[m 1 , . . . , m |M| ] T , not- ing v T the transpose of v.
Assume that we know the M j , but we do not know the matrices Qj . The equality j P j (M j Qj (M j ) T ) = R directly translates into a system (S) of affine linear equalities over the coefficients of the Qj : j (M j Qj (M j ) T )P j -R is the zero polynomial, so its coefficients, which are affine linear combinations of the coefficients of the Qj matrices, should be zero; each of these combinations thus yields an affine linear equation. The additional requirement is that the Qj are positive semidefinite.

One can equivalently express the problem by grouping these matrices into a block diagonal matrix Q and express the system (S) of affine linear equations over the coefficients of Q. By exact rational linear arithmetic, we can obtain a system of generators for the solution set of (S): Q ∈ -F 0 + vect(F 1 , . . . , F m ). The problem is then to find a positive semidefinite matrix within this search space; that is, find α 1 , . . . , α m such that -F 0 + i α i F i 0. This is the problem of semidefinite programming: finding a positive semidefinite matrix within an affine linear variety of symmetric matrices, optionally optimizing a linear form [START_REF] Vandenberghe | Semidefinite Programming[END_REF][START_REF] Boyd | Convex Optimization[END_REF].

For instance, the second unsatisfiability witness we gave for constraint system 2 is defined, using monomials {1, y}, 1 and {1, y}, by:

      2 3 0 0 1 3 1 3 0 0 0 0      
It looks like finding a solution to Equ. 3 just amounts to a SDP problem. There are, however, several problems to this approach:

1. For the general Positivstellensatz witness problem, the set of polynomials to consider is exponential in the number of inequalities.

2. Except for the simple problem of proving that a given polynomial is a sum of squares, we do not know the degree of the Q j in advance, so we cannot5 choose finite sets of monomials M j . The dimension of the vector space for Q j grows quadratically in |M j |.

3. Some SDP algorithms can fail to converge if the problem is not strictly feasible -that is, the solution set has empty interior, or, equivalently, is not full dimensional (that is, it is included within a strict subspace of the search space).

4. SDP algorithms are implemented in floating-point. If the solution space is not full dimensional, they tend to provide solutions Q that are "almost" positive semidefinite (all eigenvalues greater than -ǫ for some small positive ǫ), but not positive semidefinite.

Regarding problem 1, bounds on degrees only matter for the completeness of the refutation method: we are guaranteed to find the certificate if we look in a large enough space. They are not needed for soundness: if we find a correct certificate by looking in a portion of the huge search space, then that certificate is correct regardless. This means that we can limit the choice of monomials in M j and hope for the best.

Regarding the second and third problems : what is needed is a way to reduce the dimension of the search space, ideally up to the point that the solution set is full dimensional. As recalled by [START_REF] Peyrl | Computing sum of squares decompositions with rational coefficients[END_REF], in a sum-of-square decomposition of a polynomial P , only monomials x α1 1 . . . x αn n such that 2(α 1 , . . . , α n ) lies within the Newton polytope6 of P can appear [START_REF] Reznick | Extremal PSD forms with few terms[END_REF]Th. 1]. This helps reduce the dimension if P is known in advance (as in a sum-of-squares decomposition to prove positivity) but does not help for more general equations.

Kaltofen et al. [START_REF] Kaltofen | Exact certification in global polynomial optimization via sums-of-squares of rational functions with rational coefficients[END_REF] suggest solving the SDP problem numerically and looking for rows with very small values, which indicate useless monomials that can be safely removed from the basis; in other words, they detect "approximate kernel vectors" from the canonical basis. Our method is somehow a generalization of theirs: we detect kernel vectors whether or not they are from the canonical basis.

In the next section, we shall investigate the fourth problem: how to deal with solution sets with empty interior.

How to Deal with Degenerate Cases

In the preceding section, we have shown how to reduce the problem of finding unsatisfiability witnesses to a SDP feasibility problem, but pointed out one crucial difficulty: the possible degeneracy of the solution set. In this section, we explain more about this difficulty and how to work around it.

Let K be the cone of positive semidefinite matrices. We denote by M 0 a positive semidefinite matrix M , by M ≻ 0 a positive definite matrix M . The vector y is decomposed into its coordinates y i . x denotes a floating-point value close to an ideal real value x.

We consider a SDP feasibility problem: given a family of symmetric matrices F 0 , F i , . . . , F m , find (y i ) 1≤i≤m such that

F (y) = -F 0 + m i=1 y i F i 0. (4)
The F i have rational coefficients, and we suppose that there is at least one rational solution for y such that F (y) 0. The problem is how to find such a solution.

If nonempty, the solution set S ⊆ R m for the y, also known as the spectrahedron, is semialgebraic, convex and closed; its boundary consists in y defining singular positive semidefinite matrices, its interior are positive definite matrices. We say that the problem is strictly feasible if the solution set has nonempty interior. Equivalently, this means that the convex S has dimension m.

Interior point methods used for semidefinite feasibility, when the solution set has nonempty interior, tend to find a solution ỹ in the interior away from the boundary. Mathematically speaking, if ỹ is a numerical solution in the interior of the solution set, then there is ǫ > 0 such that for any y such that yỹ ≤ ǫ, y is also a solution. Choose a very close rational approximation y of ỹ, then unless we are unlucky (the problem is almost degenerate and all any suitable ǫ is extremely small), then y is also in the interior of S. Thus, F (y) is a solution of problem 4.

This is why earlier works on sums-of-square methods [START_REF] Peyrl | Computing sum of squares decompositions with rational coefficients[END_REF] have proposed finding rational solutions only when the SDP problem is strictly feasible. In this article, we explain how to do away with the strict feasibility clause. Some problems are not strictly feasible. Geometrically, this means that the linear affine space {-F 0 + m i=1 y i F i | (y 1 , . . . , y m ) ∈ R m } is tangent to the semidefinite positive cone K. Alternatively, this means that the solution set is included in a strict linear affine subspace of R m . Intuitively, this means that we are searching for the solution in "too large a space"; for instance, if m = 2 and y lies in a plane, this happens if the solution set is a point or a segment of a line. In this case, some SDP algorithms may fail to converge if the problem is not strictly feasible, and those that converge, in general, will find a point slightly outside the solution set. The main contribution of this article is a workaround for this problem.

Simplified algorithm

We shall thus now suppose the problem has empty interior.

The following result is crucial but easily proved:

Lemma 2. Let E be a linear affine subspace of the n × n symmetric matrices such that E ∩ K = ∅. F in the relative interior I of E ∩ K. Then it follows:

1. For all F ′ ∈ E ∩ K, ker F ⊆ ker F ′ .

The least affine space containing

E ∩ K is H = {M ∈ E | ker M ⊇ ker F }.
Suppose we have found a numerical solution ỹ, but it is nearly singularmeaning that it has some negative eigenvalues extremely close to zero. This means there is v = 0 such that |v.F (ỹ)| ≤ ǫ v . Suppose that ỹ is very close to a rational solution y and, that v.F (y) = 0, and also that y is in the relative interior of S -that is, the interior of that set relative to the least linear affine space containing S. Then, by lemma 2, all solutions F (y ′ ) also satisfy v.F (y ′ ) = 0. Remark that the same lemma implies that either there is no rational solution in the relative interior, or that rational solutions are dense in S.

How can finding such a v help us? Obviously, if v ∈ m i=0 ker F i , its discovery does not provide any more information than already present in the linear affine system -F 0 + Vect (F 1 , . . . , F m ). We thus need to look for a vector outside that intersection of kernels; then, knowing such a vector will enable us to reduce the dimension of the search space from m to m ′ < m.

Thus, we look for such a vector in the orthogonal complement of m i=0 ker F i , which is the vector space generated by the rows of the symmetric matrices F 0 , . . . , F m . We therefore compute a full rank matrix B whose rows span the exact same space; this can be achieved by echelonizing a matrix obtained by stacking F 0 , . . . , F m . Then, v = wB for some vector w. We thus look for w such that G(y).w = 0, with G(y) = BF (y)B T .

The question is how to find such a w with rational or, equivalently, integer coefficients. Another issue is that this vector should be "reasonable" -it should not involve extremely large coefficients, which would basically amplify the floating-point inaccuracies.

We can reformulate the problem as: find w ∈ Z m \ {0} such that both w and G(ỹ).w are "small", two constraints which can be combined into a single objective to be minimized α 2 G(ỹ).w 2 2 + w 2 2 , where α > 0 is a coefficient for tuning how much we penalize large values of G(ỹ).w 2 in comparison to large values of w 2 . If α is large enough, the difference between αG(ỹ) and its integer rounding M is small. We currently choose α = α 0 / G(ỹ) , with M the Frobenius norm of M (the Euclidean norm for n × n matrices being considered as vectors in R n 2 ), and α 0 = 10 15 .

We therefore try searching for a small (with respect to the Euclidean norm) nonzero vector that is an integer linear combination of the l i = (0, . . . , 1, . . . , 0, m i ) where m i is the i-th row of M and the 1 is at the i-th position. Note that, because of the diagonal of ones, the l i form a free family.

This problem is known as finding a short vector in an integer lattice, and can be solved by the Lenstra-Lenstra-Lovász (LLL) algorithm. This algorithm outputs a free family of vectors s i such that s 1 is very short. Other vectors in the family may also be very short.

Once we have such a small vector w, using exact rational linear algebra, we can compute

F ′ 0 , . . . , F ′ m ′ such that    -F ′ 0 + m ′ i=1 y ′ i F ′ i | (y1, . . . , y m ′ ) ∈ R m ′    = -F0 + m i=1 yiFi | (y1, . . . , ym) ∈ R m ∩ {F | F.v = 0} (5)
The resulting system has lower search space dimension m ′ < m, yet the same solution set dimension. By iterating the method, we eventually reach a search space dimension equal to the dimension of the solution set.

If we find no solution F ′ 0 , then it means that the original problem had no solution (the Positivstellensatz problem has no solution, or the monomial bases were too small), or that a bad vector v was chosen due to lack of numerical precision. This is the only bad possible outcome of our algorithm: it may fail to find a solution that actually exists; in our experience, this happens only on larger problems (search space of dimension 3000 and more), where the result is sensitive to numerical roundoff. In contrast, our algorithm may never provide a wrong result, since it checks for correctness in a final phase.

More Efficient Algorithm

In lieu of performing numerical SDP solving on F = -F 0 + y i F i 0, we can perform it in lower dimension on -(BF 0 B T ) + y i (BF i B T ) 0. Recall that the rows of B span the orthogonal complement of m i=0 ker F i , which is necessarily included in ker F ; we are therefore just leaving out dimensions that always provide null eigenvalues.

The reduction of the sums-of-squares problem (Eq. 3) provides matrices with a fixed block structure, one block for each P j : for a given problem all matrices F 0 , F 1 , . . . , F m are block diagonal with respect to that structure. We therefore perform the test for positive semidefiniteness of the proposed F (y) solution block-wise (see Sec. 3.6 for algorithms). For the blocks not found to be positive semidefinite, the corresponding blocks of the matrices B and F (ỹ) are computed, and LLL is performed.

As described so far, only a single v kernel vector would be supplied by LLL for each block not found to be positive semidefinite. In practice, this tends to lead to too many iterations of the main loop: the dimension of the search space does not decrease quickly enough. We instead always take the first vector v (1) of the LLL-reduced basis, then accept following vectors (1) 2 . For practical uses, we took β = γ = 10.

v (i) if v (i) 1 ≤ β. v (1) 1 and G(ỹ).v (i) 2 ≤ γ. G(ỹ).v
When looking for the next iteration ỹ′ , we use the ỹ from the previous iteration as a hint: instead of starting the SDP search from an arbitrary point, we start it near the solution found by the previous iteration. We perform leastsquare minimization so that -

F ′ 0 + m ′ i=1 y ′ i F ′ i is the best approximation of -F 0 + m i=1 y i F i | (y 1 , . . . , y m ).

Extensions and Alternative Implementation

As seen in §4, our algorithm tends to produce solutions with large numerators and denominators in the sum-of-square decomposition. We experimented with methods to get F (y ′ ) ≈ F (y) such that F (y ′ ) has a smaller common denominator. This reduces to the following problem: given v ∈ f 0 + vect(f 1 , . . . , f n ) a real (floating-point) vector and f 0 , . . . , f n rational vectors, find y ′ 1 , . . . , y n such that v ′ = f 0 + i y ′ i f i ≈ v and the numerators of v ′ have a tunable magnitude (parameter µ). One can obtain such a result by LLL reduction of the rows of:

M =      Z(βµ(f 0 -v)) Z(βf 0 ) 1 0 . . . 0 Z(βµf 1 ) Z(βf 1 ) 0 1 . . . . . . . . . 0 . . . Z(βµf n ) Z(βf n ) 0 1      (6)
where β is a large parameter (say, 10 19 ) and Z(v) stands for the integer rounding of v. After LLL reduction, one of the short vectors in the basis will be a combination n i y i l i where l 0 , . . . , l n are the rows of M , such that y 0 = 0. Because of the large βµ coefficient, y 0 (f 0v) + n i=1 y i f i should be very small, thus f 0 + n i=1 y i f i ≈ v. But among those vectors, the algorithm chooses one such that n i=0 y i f i is not large -and among the suitable v ′ , the vector of numerators is proportional to n i=0 y i f i . After computing such a y ′ , we check whether F (y ′ ) 0; we try this for a geometrically increasing sequence of µ and stop as soon as we find a solution. The matrices Qj then have simpler coefficients than the original ones. Unfortunately, it does not ensue that the sums of square decompositions of these matrices have small coefficients.

An alternative to finding some kernel vectors of a single matrix would be to compute several floating-point matrices, for instance obtained by SDP solving with optimization in multiple directions, and find common kernel vectors using LLL.

Sub-algorithms and Implementation

The reduction from the problem expressed in Eq. 3 to SDP with rational solutions was implemented in Sage. 7Solving the systems of linear equations (S) (Sec. 3.1, over the coefficients of the matrices) and 5, in order to obtain a system -F 0 + vect(F 1 , . . . , F m ) of generators of the solution space, is done by echelonizing the equation system (in homogeneous form) in exact arithmetic, then reading the solution off the echelon form. The dimension of the system is quadratic in the number of monomials (on the problems we experimented with, dimensions up to 7900 were found); thus efficient algorithms should be used. In particular, sparse Gaussian elimination in rational arithmetic, which we initially experimented, is not efficient enough; we thus instead use a sparse multi-modular algorithm [27, ch. 7] from LinBox8 . Multi-modular methods compute the desired result modulo some prime numbers, and then reconstruct the exact rational values.

One can test whether a symmetric rational matrix Q is positive semidefinite by attempting to convert it into its Gaussian decomposition, and fail once one detects a negative diagonal element, or a nonzero row with a zero diagonal element (Appendix. A). We however experimented with three other methods that perform better:

• Compute the minimal polynomial of Q using a multi-modular algorithm [START_REF] Adams | Signature of symmetric rational matrices and the unitary dual of Lie groups[END_REF]. The eigenvalues of Q are its roots; one can test for the presence of negative roots using Descartes' rule of signs. Our experiments seem to show this is the fastest exact method.

• Compute the characteristic polynomial of Q using a multimodular algorithm [START_REF] Adams | Signature of symmetric rational matrices and the unitary dual of Lie groups[END_REF] and do as above. Somewhat slower but more efficient than Gaussian decomposition.

• Given a basis B of the span of Q, compute the Cholesky decomposition of B T QB by a numerical method. This decomposition fails if and only if B T QB is not positive definite (up to numerical errors), thus succeeds if and only if Q is positive semidefinite (up to numerical errors).

For efficiency, instead of computing the exact basis B of the span of Q, we use B from §3.3, whose span includes the span of Q. The only risk is that ker B ker Q while Q is positive semidefinite, in which case B T QB will have nontrivial nullspace and thus will be rejected by the Cholesky decomposition. This is not a problem in our algorithm: it just means that the procedure for finding kernel vectors by LLL will find vectors in ker Q \ ker B.

One problem could be that the Cholesky decomposition will incorrectly conclude that B T QB is not positive definite, while it is but has very small positive eigenvalues. In this case, our algorithm may then find kernel vectors that are not really kernel vectors, leading to an overconstrained system and possibly loss of completeness. We have not encountered such cases.

Another problem could be that a Cholesky decomposition is obtained from a matrix not positive semidefinite, due to extremely bad numerical behavior. At worst, this will lead to rejection of the witness when the allegedly semidefinite positive matrices get converted to sums of squares, at the end of the algorithm.

Numerical SDP solving is performed using DSDP 9 [4, 3], communicating using text files. LLL reduction is performed by fpLLL. 10 Least square projection is performed using Lapack's DGELS.

The implementation is available from the first author's Web page (http://bit.ly/fBNLhR and http://bit.ly/gPXNF8).

Preliminary Reductions

The more coefficients to find there are, the higher the dimension is, the longer computation times grow and the more likely numerical problems become. Thus, any cheap technique that reduces the search space is welcome.

If one looks for witnesses for problems involving only homogeneous polynomials, then one can look for witnesses built out of a homogeneous basis of monomials (this technique is implemented in our positivity checker).

One could also make use of symmetries inside the problem. For instance, if one looks for a nonnegativity witness P = N/D of a polynomial P , and P is symmetric (that is, there exists a substitution group Σ for the variables of P such that P.σ = P for σ ∈ Σ), then one may reduce the search to symmetric N and D. If P = N/D is a witness, then DP = N thus for any σ, (D.σ)P = (N.σ) and thus ( σ D.σ)P = ( σ N.σ), thus D ′ = σ D.σ and N ′ = σ N.σ constitute a symmetric nonnegativity witness.

Examples

The following system of inequalities has no solution (neither Redlog nor QepCad nor Mathematica 5 can prove it; Mathematica 7 can):

       P 1 = x 3 + xy + 3y 2 + z + 1 ≥ 0 P 2 = 5z 3 -2y 2 + x + 2 ≥ 0 P 3 = x 2 + y -z ≥ 0 P 4 = -5x 2 z 3 -50xyz 3 -125y 2 z 3 + 2x 2 y 2 + 20xy 3 + 50y 4 -2x 3 -10x 2 y -25xy 2 -15z 3 -4x 2 -21xy -47y 2 -3x -y -8 ≥ 0 ( 7
)
9 DSDP is a sdp tool available from http://www.mcs.anl.gov/DSDP/ 10 fpLLL is a LLL library from Damien Stehlé et al., available from http://perso.ens-lyon.fr/damien.stehle/ form P 1 ≥ 0 ∧ • • • ∧ P n ≥ 0 ⇒ R ≥ 0. In order to prove such implication, we looked for witnesses consisting of sums of squares (Q 1 , . . . , Q n , Q R ),such that j Q j P j + Q R R = 0 with Q R = 0, and thus R = j Qj Pj QR . In some cases, it was necessary to use the products i P wi i for w ∈ {0, 1} n instead of the P i . We could find witnesses for all those problems, 11 . though for some of them, the witnesses are very large, taking up megabytes. Since these searches were done without making use of symmetries in the problem, it is possible that more clever techniques could find smaller witnesses.

Conclusion and further works

We have described a method for solving SDP problems in rational arithmetic. This method can be used to solve sums-of-squares problems even in geometrically degenerate cases. We illustrated this method with applications to proving the nonnegativity of polynomials, or the unsatisfiability of systems of polynomial (in) equalities. The method then provides easily checkable proof witnesses, in the sense that checking the witness only entails performing polynomial arithmetic and applying a few simple mathematical lemmas. We have implemented the conversion of nonnegativeness witnesses to Coq proofs. A more ambitious implementation, mapping Coq real arithmetic proofs goals to Positivstellensatz problems through the Psatz tactic from the MicroMega package [START_REF] Besson | Fast Reflexive Arithmetic Tactics the Linear Case and Beyond[END_REF], then mapping Positivstellensatz witnesses back to proofs, is underway.

One weakness of the method is that it tends to provide "unnatural" witnesses -they tend to have very large coefficients. These are machine-checkable but provide little insights to the reader. An alternative would be to provide the matrices and some additional data (such as their minimal polynomial) and have the checker verify that they are semidefinite positive; but this requires formally proving, once and for all, some non-trivial results on polynomials, symmetric matrices and eigenvalues (e.g. the Cayley-Hamilton theorem), as well as possibly performing costly computations, e.g. evaluating a matrix polynomial.

A more serious limitation for proofs of unsatisfiability is the very high cost of application of the Positivstellensatz. There is the exponential number of polynomials to consider, and the unknown number of monomials. It would be very interesting if there could be some simple results, similar to the Newton polytope approach, for reducing the dimension of the search space or the number of polynomials to consider. Another question is whether it is possible to define SDP problems from Positivstellensatz equations for which the spectrahedron has rational points only at its relative boundary.

While our method performed well on examples, and is guaranteed to provide a correct answer if it provides one, we have supplied no completeness proofthat is, we have not proved that it necessarily provides a solution if there is one. This is due to the use of floating-point computations. One appreciable result would be that a solution should be found under the assumption that floating-point computations are precise up to ǫ, for a value of ǫ and the various scaling factors in the algorithm depending on the values in the problem or the solution.

It seems possible to combine our reduction method based on LLL with the Newton iterations suggested by [START_REF] Kaltofen | Exact certification in global polynomial optimization via sums-of-squares of rational functions with rational coefficients[END_REF][START_REF] Kaltofen | Exact certification of global optimality of approximate factorizations via rationalizing sums-of-squares with floating point scalars[END_REF], as an improvement over their strategy for detection of useless monomials and reduction of the search space. Again, further experiment is needed. Suppose that the entrance of iteration i, m[i . . . n, i . . . n] is positive semidefinite. If m i,i = 0, then the ith base vector is in the isotropic cone of the matrix, thus of its kernel, and the row i must be zero. Otherwise, m i,i 0. By adding ǫ to the diagonal of the matrix, we would have a positive definite matrix and thus the output of the loop iteration would also be positive definite, as above. By ǫ → 0 and the fact that the set of positive semidefinite matrices is topologically closed, then the output of the loop iteration is also positive semidefinite.

The output variable is then a list of couples (c i , v i ) such that c i > 0 and the original matrix m is equal to i c i v T i v i (with v i row vectors). Otherwise said, for any row vector u, umu T = i c i u, v i 2 .

  throw n o n p o s i t i v e s e m i d e f i n i t e e l s e b e g i n v := m. row ( i ) / m[ i , i ] o utput . append (m[ i , i ] , v ) m := mm[ i , i ] * v . t r a n s p o s e ( ) * v end end

http://pvs.csl.sri.com/

http://www.astree.ens.fr/

There exists a theoretical exact algorithm for computing such a decomposition for homogeneous polynomials of at most 3 variables[START_REF] De Klerk | Products of positive forms, linear matrix inequalities, and Hilbert 17th problem for ternary forms[END_REF]; we know of no implementation of it and no result about its practical usability.

Another result, due to Schmüdgen[START_REF] Schweighofer | On the complexity of Schmüdgen's Positivstellensatz[END_REF], gives simpler witnesses forP 1 ≥ 0 ∧ • • • ∧ Pn ≥ 0 ⇒ C in the case where P 1 ≥ 0 ∧ • • • ∧ Pn ≥ 0 defines a compact set.

There exist non-elementary bounds on the degree of the monomials needed[START_REF] Lombardi | Une borne sur les degrés pour le Théorème des zéros réel effectif[END_REF]. In the case of Schmüdgen's result on compact sets, there are better bounds[START_REF] Schweighofer | On the complexity of Schmüdgen's Positivstellensatz[END_REF].

The Newton polytope of a polynomial P , or in Reznick's terminology, its cage, is the convex hull of the vertices (α 1 , . . . , αn) such that x α 1 1 . . . x αn n is a monomial of P .

Sage is a computer algebra system implemented using the Python programming language, available under the GNU GPL from http://www.sagemath.org.

LinBox is a library for exact linear arithmetic, used by Sage for certain operations. http://www.linalg.org/

A 7z archive is given at http://bit.ly/hM7HW3

* This work was partially supported by ANR project "ASOPT". † CNRS / VERIMAG

This system was concocted by choosing P 1 , P 2 , P 3 somewhat haphazardly and then P 4 = -(P 1 +(3+(x+5y) 2 )P 2 +P 3 +1+x 2 ), which guaranteed the system had no solution. The initial 130 constraints yield a search space of dimension 145, and after four round of numeric solving one gets an unsatisfiability witness (sums of squares Q j such that 4 j=1 P j Q j + Q 5 = 0). Total computation time was 4.4 s. Even though there existed a simple solution (note the above formula for P 4 ), our algorithm provided a lengthy one, with large coefficients (and thus unfit for inclusion here).

Motzkin's polynomial M (Eq. 1) cannot be expressed as a sum of squares, but it can be expressed as a quotient of two sums of squares. We solved M.Q 1 -Q 2 = 0 for sums of squares Q 1 and Q 2 built from homogeneous monomials of respective total degrees 3 and 6 -lesser degrees yield no solutions (Fig. 1). The equality relation over the polynomials yields 66 constraints over the matrix coefficients and a search space of dimension 186. Four cycles of SDP programming and LLL are then needed, total computation time was 4.1 s.

We exhibited witnesses that each of the 8 semidefinite positive forms listed by [START_REF] Reznick | Some Concrete Aspects Of Hilbert's 17th Problem[END_REF], which are not sums of squares of polynomials, are quotients of sums of squares (Motzkin's M , Robinson's R and f , Choi and Lam's F , Q, S, H and Schmüdgen's q). These examples include polynomials with up to 6 variables and search spaces up to dimension 1155. We did likewise with delzell, laxlax and leepstarr2 from [START_REF] Kaltofen | Exact certification in global polynomial optimization via sums-of-squares of rational functions with rational coefficients[END_REF]. The maximal computation time was 7'.

We then converted these witnesses into Coq proofs of nonnegativity using a simple Sage script. These proofs use the Ring tactic, which checks for polynomial identity. Most proofs run within a few seconds, though laxlax takes 7'39" and Robinson's f 5'07"; the witness for leepstarr2 is too large for the parser. We also exhibited a witness that the Vor1 polynomial cited by [START_REF] Safey | Computing the global optimum of a multivariate polynomial over the reals[END_REF] is a sum of squares.

John Harrison kindly provided us with a collection of 14 problems that his system [START_REF] Harrison | Verifying nonlinear real formulas via sums of squares[END_REF] could not find witnesses for. These problems generally have the

A Gaussian Reduction and Positive Semidefiniteness

An algorithm for transforming a semidefinite positive matrix into a "sum of squares" form, also known as Gaussian reduction:

f o r i :=1 t o n do b e g i n i f m[ i , i ] < 0 then throw n o n p o s i t i v e s e m i d e f i n i t e i f m[ i , i ] = 0 then i f m. row ( i ) <> 0


Thanh Phuong Nguyen

Isabelle Debled-Rennesson

Arc segmentation in linear time

. By utilizing this representation, we transform the problem of digital arc detection into a problem of digital straight line recognition. We then develop a linear method for arc segmentation of digital curves.

Introduction

The digital arcs and circles are basic geometric objects of which the recognition is an interesting topic. In the literature, some methods have been proposed for the recognition of digital circles. Nakamura et al. [START_REF] Nakamura | Digital circles[END_REF] proposed a recursive algorithm for determining the center of a digital circle, but its complexity is exponential in the general case. Kim [START_REF] Kim | Digital disks[END_REF][START_REF] Kim | Digital disks and a digital compactness measure[END_REF] proposed several results on digital disks. The first result [START_REF] Kim | Digital disks[END_REF] detects if a set of grid points in a N × N image is a digital disk with complexity O(n 3 ). The second result [START_REF] Kim | Digital disks and a digital compactness measure[END_REF] reduces this task to O(n 2 ). Based on the classical separating arc problem, Kovalevsky [START_REF] Kovalevsky | New definition and fast recognition of digital straight segments and arcs[END_REF] (resp. Fisk [START_REF] Fisk | Separating point sets by circles, and the recognition of digital disks[END_REF]) proposed an algorithm for the recognition of a digital disk in O(n 2 log n) (resp. O(n 2 )) time. Coeurjolly [START_REF] Coeurjolly | An elementary algorithm for digital arc segmentation[END_REF] transformed the problem of circle recognition into a problem of search a 2D point that belongs to the intersection of n 2 half-plane. Sauer [START_REF] Sauer | On the recognition of digital circles in linear time[END_REF] (resp. Damaschke [START_REF] Damaschke | The linear time recognition of digital arcs[END_REF]) presented a linear algorithm to decide if a curve is a digital circle (resp. arc) based on Megiddo's algorithm [START_REF] Megiddo | Linear programming in linear time when the dimension is fixed[END_REF]. Worring [START_REF] Worring | Digitized circular arcs: characterization and parameter estimation[END_REF] introduced a digital circle segmentation method by using a fixed size window process. Roussillon [START_REF] Roussillon | On three constrained versions of the digital circular arc recognition problem[END_REF] proposed a linear algorithm of circle recognition in 3 particular cases.

We present in this paper a linear method for the detection of digital circles or digital arcs based on a discrete geometry approach. Firstly, a polygonalization is applied in linear time on the input curve [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF]. Secondly, we use a transform proposed by Latecki et al. [START_REF] Latecki | Shape similarity measure based on correspondence of visual parts[END_REF] to represent the obtained polygon in a novel space called tangent space. We show that a sequence of chords of a circle will correspond to a sequence of collinear points in the tangent space. So the problem of arc/circle detection can be considered as a problem of digital straight line recognition.

This paper is organized as follows. Section 2 recalls some definitions concerning digital circles and blurred segments. The next section presents a technique to transform an arc into the tangent space and proposes some principal properties of the arc in this representation. Section 4 proposes a linear algorithm for the detection of digital arcs or digital circles. In Section 5, we present a linear method for the segmentation of a curve into arcs and some experimentations.

Discrete circle and blurred segment

Discrete circle: In the literature, there exist several definitions of discrete circle. They are proposed by considering a real circle on the grid digitization. The difference among them is the process of discretization. Nakamura et al. [START_REF] Nakamura | Digital circles[END_REF] considered a discretization of a real circle by the points of Z 2 that are the nearest points of that circle. Kim [START_REF] Kim | Digital disks[END_REF] proposed a definition of discrete circle as a boundary of a digital disk superimposed by a real circle. Andres [START_REF] Andres | Discrete circles, rings and spheres[END_REF] used an arithmetic approach to define a digital circle as a sequence of points superimposed by a ring. Discrete line and blurred segment: The notion of blurred segment [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] was introduced from the notion of arithmetical discrete line. An arithmetical discrete line, noted D(a, b, µ, ω), is the set of points (x, y) ∈ Z 2 that verifies: µ ≤ axby < µ + ω with a main vector (b, a), lower bound µ and thickness ω . A width ν blurred segment (BS) is a set of points (x, y) ∈ R 2 that is optimally bounded (see [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] for more details) by a discrete line D(a, b, µ, ω) verifying ω-1 max(|a|,|b|) ≤ ν. Fig. 1 shows a BS of with 1.25 (the sequence of gray points) whose optimal bounding line is D (5, 8, -8, 11). A linear method for recognition of BS has been also proposed in [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF].

Arc representation in tangent space

Modified tangent space representation: We recall in this section some notions concerning a representation of a polygon in the tangent space. Latecki et al. [START_REF] Latecki | Shape similarity measure based on correspondence of visual parts[END_REF] proposed the tangent space representation as a tool of similarity measure for shape matching. Inspired from this representation, we propose a modified tangent space to represent a polygonal curve. The difference is that we do not normalize the axis 0x in the tangent space. Let C = {Ci} n i=0 be a polygonal curve, αi = ∠(

----→ Ci-1Ci, ----→ CiCi+1
) and li -the length of the line segment CiCi+1, i ∈ {0, . . . , n -1}. If Ci+1 is on the right of ----→ Ci-1Ci then αi > 0, otherwise αi < 0. From now, we denote P.x (resp. P.y) to indicate the x (resp. y)-coordinate of point P . We consider a transformation that associates the polygon C of Z 2 to a polygon of R 2 that is constituted by line segments Ti2T (i+1)1 , T (i+1)1 T (i+1)2 for i from 0 to n -1 with T02 = (0, 0), Ti1 = (T (i-1)2 .x + li-1, T (i-1)2 .y), i from 1 to n, Ti2 = (Ti1.x, Ti1.y + αi), i from 1 to n -1.

y x

Fig. 1. A BS [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF].

α 1 α 2 α 3 C0 C1 C2 C3 C4 0 y x α 2 α 3 T 32 T 41 T 02 T 11 T 22 T 31 T 12 T 21 α 1
Fig. 2. Transformation to the tangent space: on the left, the input polygonal curve and on the right, its tangent space representation.

Properties of arcs in the modified tangent space representation The theorem below allows us to study the properties of a representation in the tangent space of a polygon that corresponds to an arc or a circle.

Theorem 1. Let C = {Ci} n i=0 be a polygon, αi = ∠( ----→ Ci-1Ci, ----→ CiCi+1).
The length of CiCi+1 is li, for i ∈ {0, . . . , n -1}. The vertices of C are on a real arc of radius R and of center O such that ∠CiOCi+1 ≤ π 4 for i ∈ {1, . . . , n -1}. This results below is obtained.

1 R < αi l i +l i+1 2 < 1 0.9742979R Proof. Let us consider figure 3. We have αi = ∠CiOHi-1 + ∠CiOHi. We denote that αi1 = ∠CiOHi-1 and αi2 = ∠CiOHi. Moreover, ∠C1OH0 = ∠C 0 OC 1 2 ≤ π 8 , ∠C1OH1 = ∠C 1 OC 2 2 ≤ π 8 .
In addition, we have sin ∠C1OH0 = l 0 2R , sin ∠C1OH1 = l 1 2R . Therefore, l 0 +l 1 2R = sin α11 + sin α12. Similarly, we have

l i-1 +l i 2R = sin αi1 + sin αi2, for i ∈ {1, . . . , n -1}. Because of x ≥ sin x ≥ x -x 3 6 with x > 0, we have αi1 > sin αi1 > αi1(1 - α 2 i1 6 ) > αi1(1 - π 8 2 6 ) > 0.9742979αi1. Similarly αi2 > sin αi2 > 0.9742979αi2. Therefore, we have αi > l i-1 +l i 2R > 0.9742979αi. This theorem is proved. O C2 C3 H0 H1 C0 C5 C4 C1 α1 H Fig. 3.
Property of a set of sequential chords of a partial circle.

x y 0 This theorem allows to deduce that the corresponding curve of midpoints of T (i-1)2 Ti1, 1 ≤ i ≤ n in the tangent space of the curve C is quasi collinear. From now on, the midpoint curve is called M pC. In addition, the more sin αi closes to αi, 1 ≤ i < n, the more M pC is collinear. Therefore, we can decide if a digital curve approximates an arc of circle by verifying the collinearity of its M pC in the tangent space. A qualitatif study on this approximation will be also considered in Section 5.

T (i+1)1 T (i+2)1 T i1 T i2 T (i-1)2 M i+1 I i I i+1 M i T (i+1)2 M i-1

Arc segmentation

Detection of digital arcs Thanks to theorem 1, we introduce now an heuristic algorithm for (see algo. 1) deciding if a digital curve is an arc. Our main idea is to work on the representation of a digital curve in the modified tangent space. In this representation, the set of midpoints M pC = {Mi} n-1 i=0 (see the above section) will be constructed. And we will use a linear procedure [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] to test the collinearity of these points. If the response is positive, we consider that the input digital curve is a digital arc (a partial circle).

A higher value of the range of vertical ordinate in the tangent space leads to false positive detection as an example; an helix can be detected as an arc. So, to avoid this problem, the maximal difference of vertical ordinate is fixed to 2π for detecting an arc. The input parameter αmax of Algorithm 1 allows to control the obtained error.Parameter ν1 is used for polygonalization by using the recognition of blurred segments. The input parameter ν2 is used as the width in the algorithm for recognition of blurred segments [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] to test the collinearity of the midpoint set in the representation of the tangent space. In practice, αmax (resp. ν1) is chosen as π 4 (resp. 1). In addition, ν2 can be chosen as a fixed value from 0.15 to 0.25 without problem.

This heuristic algorithm detect well circular arc. In Section 5 (see corollary 1), we consider an adaptive estimation of ν2 to guarantee the quality of detected arcs. It is Use algorithm [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] to polygonalize P into BS of width ν1: C = {C} m i=0 ; Represent C in the modified tangent space by T (C); BS = ∅; if there exists i such that |Ti2.y -Ti1.y| > αmax then return FALSE; Determine midpoint set M pC = {Mi} m-1 i=0 of {Ti2T (i+1)1 } m-1 i=0 ; i = 1; while |Mi.y -M0.y| ≤ 2π do BS = BS ∪ Mi; i++; Use algorithm [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] The algorithms of polygonalization and of blurred segment recognition are linear. The complexity of the tangent space transform, and the construction of midpoint curve M pC = {Mi} m-1 i=1 is in O(m). Because m << n, the total complexity of our detection method is then in O(n). Segmentation of curves into digital arcs Based on the above idea for detection of an arc, we then develop a linear method for the segmentation of digital arcs by using a width ν blurred segment [START_REF] Debled-Rennesson | Optimal blurred segments decomposition of noisy shapes in linear time[END_REF] polygonalization on the curve of midpoints. Its main idea, illustrated in figure 5, is based on the polygonalization of the midpoint curve (Fig. 5.e). Contrariwise to Algorithm 1, we polygonalize the midpoint curve M pC in spite of recognizing if M pC is a BS and then each line segment corresponds to a circular arc. Experimental results and application to real images We have implemented this linear method. An example of a curve segmented into arcs is presented in Fig. 5. Firstly, the approximating polygon (see Fig. 5.b) is constructed from the input curve in Fig. 5.a. After that, we transform it into the modified tangent space representation (see Fig. 5.c). Then, by polygonalizing the curve of midpoints in this tangent space (see Fig. 5.d), the corresponding arcs can be detected (see Fig. 5.e). Fig. 6 shows an experimentation on technical drawing images. Figs. 6.a, 6.d are input images. Figs. 6.c, 6.f present the extracted arcs from the borders presented in 6.b, 6.e. Our method gives good results on this type of images which frequently contain arc and circle primitives. Fig. 7 presents our obtained result with a real image.

Study of quasi collinearity property of midpoint curve

Algorithm 1 works well in practice. However, we have a problem to estimate error approximation in general case when the value of ν2 is fixed. In this section, we present a first study concerning the utilization of this algorithm where ν2 is chosen adaptively. Let us suppose that αmax = max{αi} n i=1 . Let us suppose that Ri is the radius of the approximating circle that passes through 3 points Ci-1, Ci, Ci+1; αi1 = ∠Hi-1OCi, αi2 = ∠HiOCi (see Fig. 3). We suppose that αi1, αi2 ≤ π 8 for i = 1, . . . , n -1 to guarantee the condition sin x ≃ x in Theorem 1. It means that we consider the condition sin x ≃ x with x ∈ [0, π 8 ]. Therefore, we have αi ≤ π 4 .

Comparison of radius of local circumcircles to the global radius

Proposition 1. Let C = {Ci} n i=0 be a polygon, αi = ∠( ----→ Ci-1Ci, ----→ CiCi+1).
The length of CiCi+1 is li, for i ∈ {0, . . . , n -1}. We denote Oi (resp. Ri) respectively the center (resp. the radius) of circumcirle that passes through 3 points Ci-1, Ci, Ci+1, Hi the projection of Oi on CiCi+1. Suppose that Ri -OHi ≤ h for i ∈ {1, . . . , n -1}. This results below is obtained. Riαi ≥

l i-1 +l i 2 ≥ Riαi -0.3377hαi
Proof. We denote αi1 = ∠Hi-1OiCi, αi2 = ∠HiOiCi (see Fig. 3). Firstly, cos

αi1 = 1- 2 sin 2 α i1 2 = OH i R i ≥ R i -h R i = 1-h R i .
In addition, thanks to αi1 ≤ π 8 and sin(x)

x is decreas- ing in [0, π 16 ], we have sin x ≥ x sin π 16 π 16 . Therefore h R i ≥ 2 sin 2 α i 2 ≥ 2 • sin π 16 π 16 2 α i1 2 2 > 0.9872 2 α 2 i1 . Similarly, we have h R i > 0.9872 2 α 2 i2
, for i ∈ {1, . . . , n -1} (1). In addition, we have this remark

x ≥ sin x ≥ x -x 3 6 with π 4 ≥ x ≥ 0. So, αi ≥ sin αi1 + sin αi2 > αi1 + αi2 -1 6 (α 3 i1 + α 3 i2 ) = (αi1 + αi2)(1 -1 6 (α 2 i1 -αi1αi2 + α 2 i2 )) = αi(1 -1 6 (α 2 i1 -αi1αi2 + α 2 i2 )) ≥ αi(1 -1 6 (α 2 i1 + α 2 i2 ))), for i ∈ {1, . . . , n -1} (2)
. Thanks to (1) and ( 2), we obtain sin αi1 + sin αi2 > α1(1 -

1 3•0.9872 h R ), for i ∈ {1, . . . , n -1} (3).
Moreover, we have αi = αi1 + αi2 and αi1, αi2 ≤ π 8 . In addition, we have sin αi1

= l i-1
2R , sin αi2 = l i 2R . Therefore, we have

l i-1 +l i 2R
= sin αi1 + sin αi2, for 1 ≤ i < n (4). Thanks to (3) and (4), we obtain Riαi ≥

l i-1 +l i 2 ≥ Riαi(1 - 1 3•0.9872 h R i ) = αi(Ri - h 3•0.9872 ) ⇔ Riαi ≥ l i-1 +l i 2
≥ Riαi -0.3377hαi. Now we consider a set of midpoints M pC is a blurred segment whose horizontal width is ǫ. Let us suppose that the slope of this blurred segment is

1 R , R ∈ R. M i-1 M i A B ǫ
Fig. 8. ≥ Rαiǫ. This double inequations can be rewrited as

C 0 C i-1 C 1 C i O ′ i O ′′ i R R 2.1δ 2.1δ
l i-1 +l i 2 = Rαi + ǫi, where ǫi ∈ [-ǫ, ǫ]. Thank to proposition 1, we have: Riαi > Rαi + ǫi > Riαi -0.3377hαi ⇔ Ri > R + ǫ i α i > Ri -0.3377h. So, ǫ α i < Ri -R < ǫ α i + 0.3377h
. So, we have corollary 1.

Corollary 1 Let C = {Ci} n i=0 be a polygon, αi = ∠( ----→ Ci-1Ci, ----→ CiCi+1).
The length of CiCi+1 is li, for i ∈ {0, . . . , n -1}. The set of midpoints {Mi} n-1 i=0 is a blurred segment whose horizontal width is ǫ. We denote Oi (resp. Ri) respectively the center (resp. the radius) of circumcirle that passes through 3 points Ci-1, Ci, Ci+1, Hi the projection of Oi on CiCi+1. Suppose that Ri -OHi ≤ h for i ∈ {1, . . . , n -1}. This results below is obtained.

0 < Ri -R < ǫ α i + 0.3377h ≤ ǫ min{α i } n i=1 + 0.3377h.

Localization of centers of local circumcircles

In this section, we consider the convergence of local circumcircle centers if this condition below is satisfied: |Ri-R| ≤ δ, R ∈ R, 1 ≤ i, j ≤ n -1. Let us consider Fig. 10. To prove the convergence of centers of local circumcircle, we show this property for an approximation of a half of circle.

Proposition 2. Let us consider a sequence of points {C}

n i=0 . There exist R and δ such that R, δ ∈ R, 0 ≤ Ri -R ≤ δ, i = 1, . . . , n -1. Suppose that ∠C k CjCj+1 > π 2 for k ∈ {0, 1}, k < j < n. Therefore, we have this property 0 ≤ R ′ i -R ≤ δ, 0 ≤ R ′′ i -R ≤ δ, for 1 ≤ i ≤ n -1. Proof. We denote O ′ i (resp. O ′′ i ) and R ′ i (resp. R ′′ i
) the centers and radius of circumcircles that passes through 3 points C0 (resp. C1), Ci, Ci+1. Firstly, we have a trivial remark: the perpendicular bisector of CiC k is between that of CiCj and CjC k . Now, we prove this proposition by induction.

Because of R

Conclusions

We have presented a linear method for the detection of digital circles or digital arcs. A linear method for the segmentation of a curve into digital arcs is also proposed. This method is based on a discrete geometry approach. It is simple, easy and robust to implement. A more complete demonstration of the algorithm is under process.

C 0 C i-1 C i H i O ′ i-1 O i O ′ i C i+1 H ′ i-1 H ′ i H i+1 (a) Outside C 0 C i-1 C i C i+1 H ′ i-1 H ′ i H i H i +1 O ′ i O i O ′ i-1 (b) Inside

Fig. 4 .

 4 Fig. 4. Property of a polygon in our modified tangent space representation.

Algorithm 1 :

 1 Detection of a digital arc/circle Data: P = {Pi} n i=0 digital curve, αmax -maximal admissible angle, ν1-width of BS for polygonalization, ν2-width of BS for collinear test 1 Result: ARC (resp. CIRCLE): C is a digital arc (resp. circle), FALSE if not. begin

Fig. 5 .Fig. 6 .

 56 Fig. 5. Arc segmentation results on a digital curve: (a) Input curve, (b) Approximated polygon, (c) Results of arc segmentation.

Fig. 7 .

 7 Fig. 7. Experimentation on a real image at width 2.

Fig. 9 .

 9 Fig. 9. {O ′ i , O ′′ i } n-1i=1 is in a compact zone.

′ 1 = 1 1 8 ≤ 1 .

 11181 R1, the proposition is true with i = 1. Suppose that |R ′ i-1 -R| ≤ δ. Let us consider Fig. 10. We denote H ′ i , Hi are respectively the midpoint of C0Ci, Ci-1Ci. Thank to the above remark, we have O Hi. We consider now the position of Ci+1 with the circle of centerO ′ i-1 , of radius R ′ i-1 . The proposition if trivial if Ci+1 is on this circle because of R ′ i = R ′ i-1 . If Ci+1 is outside of this circle (see Fig. 10.a), we then deduce O Ci. It means that Ri+1 > R ′ i > R ′ i-1 . Thank to 0 ≤ Ri+1 -R ≤ δ and 0 ≤ R ′ i-1 -R ≤ δ, we have 0 ≤ R ′ i -R ≤ δ.In other case (see Fig.10.b), by applying the same arguments, we obtain Ri+1 < R′ i < R ′ i-1 . Therefore, we have 0 ≤ R ′ i -R ≤ δ, for 1 ≤ i ≤ n -1.By replacing C0 by C1 and using the same argument, we have 0≤ R ′′ i -R ≤ δ, for 1 ≤ i ≤ n -1.We have a simple remark: a triangle ABC that satisfies ∠ABC ≥ 7π 8 have :AC > BC + cos π 8 • AB. It is trivial because of AC 2 = BC 2 + AB 2 -2BC • AB cos ∠ABC > BC 2 + cos 2 π 8 AB 2 + 2BC • AB cos π 8. Therefore, we have: O 1δ, for i = 1, . . . , n -1. Thanks to this result and proposition 2, it is trivial now to show that set of center O ′ i , O ′′ i is in a compact zone (see Fig.9).

Fig. 10 .

 10 Fig. 10. Position between Ci+1 and circumcircle of C0Ci-1Ci.

  to verify if BS is a blurred segment of width ν2; if BS is a blurred segment of with ν2 then if |Mm-1.y -M0.y| == 2 * π then return CIRCLE; else return ARC; else return FALSE; end evident that we can lightly change Algorithm 1 to check the condition of Corollary 1 and Proposition 2 in linear time also.

By default, αmax = π 4 , ν1 = 1 for normal curves and ν2 = 0.2 (see algo. 1).