Alin Bostan
email: alin.bostan@inria.fr

Guillaume Ch Èze

Thomas Cluzeau
email: cluzeau@ensil.unilim.fr

Jacques-Arthur Weil
email: jacques-arthur.weil@unilim.fr

Guillaume Chèze
email: guillaume.cheze@math.univ-toulouse.fr

EFFICIENT ALGORITHMS FOR COMPUTING RATIONAL FIRST INTEGRALS AND DARBOUX POLYNOMIALS OF PLANAR POLYNOMIAL VECTOR FIELDS

We present fast algorithms for computing rational first integrals with bounded degree of a planar polynomial vector field. Our approach builds upon a method proposed by Ferragut and Giacomini ([FG10]) whose main ingredients are the calculation of a power series solution of a first order differential equation and the reconstruction of a bivariate polynomial annihilating this power series. We provide explicit bounds on the number of terms needed in the power series. This enables us to transform their method into a certified algorithm computing rational first integrals via systems of linear equations.We then significantly improve upon this first algorithm by building a probabilistic algorithm with arithmetic complexity Õ(N 2ω ) and a deterministic algorithm solving the problem in at most Õ(d 2 N 2ω+1 ) arithmetic operations, where N denotes the given bound for the degree of the rational first integral, and where d is the degree of the vector field, and ω the exponent of linear algebra. We also provide a fast heuristic variant which computes a rational first integral, or fails, in Õ(N ω+2 ) arithmetic operations. By comparison, the best previously known complexity was d ω+1 N 4ω+4 arithmetic operations using the algorithm given in [Chè11]. We then show how to apply a similar method to the computation of Darboux polynomials. The algorithms are implemented in a Maple package RationalFirstIntegrals which is available to interested readers with examples showing its efficiency.

Introduction

Context. Let K denote an effective field of characteristic zero, i.e, one can perform arithmetic operations and test equality of two elements (typically, K = Q or Q(α), where α is an algebraic number). Given two polynomials A, B in K[x, y], we consider the planar polynomial vector field (S) ẋ = A(x, y), ẏ = B(x, y),

and discuss the problem of computing rational first integrals of (S), i.e., rational functions F ∈ K(x, y) \ K that are constant along the solutions (x(t), y(t)) of (S). More precisely, the present article is concerned with the following algorithmic problem:

(P N ): given a degree bound N ∈ N, either compute a rational first integral F ∈ K(x, y) \ K of (S) of total degree at most N , or prove that no such F exists.

Date: December 23, 2014.

This old problem was already studied by Darboux in 1878 ([Dar78]) and has been the subject of numerous works ever since. The naive approach (by indeterminate coefficients) leads to a polynomial system of quadratic equations in the coefficients of F . Other methods use what is called nowadays Darboux polynomials, in the spirit of the celebrated Prelle-Singer's method [START_REF] Prelle | Elementary first integrals of differential equations[END_REF]; see Subsection 2.3 for a review. These methods also require solving a polynomial system of quadratic equations. Recently, Chèze [Chè11] has shown that problem (P N ) can be solved in polynomial time in N . The importance of this result is mainly theoretical since the exponent in the polynomial complexity estimate is bigger than 10.

To improve upon this current state of affairs, our starting point is the article [START_REF] Ferragut | A new algorithm for finding rational first integrals of polynomial vector fields[END_REF] of Ferragut and Giacomini. The key observation is that (S) has a rational first integral if and only if all power series solutions in K [[x]] of the first order non-linear differential equation (E) dy dx = B(x, y) A(x, y)

are algebraic over K(x). Furthermore, they show how minimal polynomials of these algebraic power series may lead to rational first integrals. The methods presented in [START_REF] Ferragut | A new algorithm for finding rational first integrals of polynomial vector fields[END_REF] involve solving a (simpler) polynomial system of quadratic equations (sometimes linear, as suggested in their Remark 1 and the example that follows).

Our main contributions. In the present article, we push further the observation of Ferragut and Giacomini, so as to give fast algorithms solving Problem (P N ). We prove structure results in Section 2.2, notably Theorem 11, which guarantees that computations may be performed in the base field without loss of information. We provide effective bounds (Lemma 18) on the number of terms needed in the involved power series in order to certify that the guessed minimal polynomial is correct. This leads us to a first algorithm in Section 3.3, showing that solving Problem (P N ) is reduced to linear equations only1 . We then turn it into a faster probabilistic algorithm that uses Õ(N 2ω ) arithmetic operations in K, where ω ∈ [2, 3] is the exponent of linear algebra over K, and the soft-O notation Õ( ) indicates that polylogarithmic factors are neglected. Lemma 14 (and part (2) of Proposition 16) bound the number of remarkable values (i.e. bad cases for our probabilistic approach). The probabilistic algorithm is then turned into a deterministic one, that solves Problem (P N ) in arithmetic complexity Õ(d 2 N 2ω+1 ), where d = max(deg(A), deg(B)) denotes the degree of the polynomial vector field (S). This compares well to the previous polynomial time algorithm given in [Chè11], which uses at least d ω+1 N 4ω+4 arithmetic operations. Note that if we take ω = 3 (i.e., the cost of naive linear algebra), then the above means that the best previously known complexity would be in Õ(d 4 N 16 ) whereas our deterministic algorithm would use at most Õ(d 2 N 7 ) arithmetic operations, and our probabilistic one would use Õ(N 6 ). Lastly, we sketch a heuristic method that uses Õ(N ω+2 ) arithmetic operations (i.e., Õ(N 5 ) using classical linear algebra) which is sub-cubic, given that the output has size O(N 2 ).

We provide algorithmic details, notably precise degree bounds and complexity estimates. The algorithms developed in the article are implemented in a Maple package called RationalFirstIntegrals2 . Using this implementation, we demonstrate the efficiency of our algorithms on some examples. Finally, we show how to apply a similar method to the computation of Darboux polynomials.

Structure of the article. In Section 2, we recall Darboux's approach to the integrability of polynomial vector fields, related works, and existing results about the problem (P N ). We also give useful facts on the so-called spectrum problem. We recall in Section 3 the connection between rational first integrals of the polynomial vector field (S) and algebraic power series solutions of (E). We then propose a first algorithm, based on linear algebra, that solves Problem (P N ). Building on this, we develop in Section 4 an efficient probabilistic algorithm, and then turn it into an efficient deterministic algorithm. In Section 5, we study the arithmetic complexity of the algorithms developed in Section 4, and discuss several algorithmic issues. Then, in Section 6 we present our implementation and display its behavior on various examples. Finally, Section 7 shows how similar ideas can be used for computing the set of all irreducible Darboux polynomials (of a given degree) of planar polynomial vector fields.

Notation. The degree deg(P ) of a bivariate polynomial P ∈ K[x, y] is the total degree of P . A rational function P/Q with P, Q ∈ K[x, y] is said to be reduced when P and Q are coprime. The degree deg(F ) of a reduced rational function F = P/Q is the maximum of deg(P ) and deg(Q). We denote by K an algebraic closure of the field K. We write ḟ := ∂f ∂t for the usual formal derivative of the "function" (polynomial, or power series) f with respect to the variable t. For a set Ω, we denote by |Ω| its cardinality.

Review on first integrals of polynomial vector fields

In this section, we recall several useful facts, mainly to keep the exposition as self-contained as possible, and to clarify the understanding of the algorithms that we develop below. Some results are not original.

2.1. First definitions and classical results. We consider an autonomous planar polynomial vector field (S) ẋ = A(x, y), ẏ = B(x, y), where x and y are unknown "functions" of the time variable t, A and B are polynomials in K[x, y], and d := max(deg(A), deg(B)) denotes the degree of the polynomial vector field. Without any loss of generality, A and B will be assumed to be coprime in the remaining of the article. To (S) is attached the derivation

D := A(x, y) ∂ ∂x + B(x, y) ∂ ∂y ,
acting on the polynomial ring K[x, y]. We thus view K[x, y] as a differential ring endowed with the derivation D. We denote by K(x, y) its field of fractions.

Definition 1. A rational first integral of (S) is a non-constant rational function

F ∈ K(x, y) \ K satisfying D(F ) = 0.
A rational first integral F of (S) is thus a non-trivial constant for the derivation D. Intuitively, this means that if (x(t), y(t)) is a pair of "functions" satisfying (S), then F (x(t), y(t)) is constant when t varies. We explain in Theorem 11 below why no algebraic extension of the base field is necessary in Definition 1.

A starting observation is that the rational function F = P/Q is a first integral for (S) if and only if D(P ) Q = D(Q) P . Therefore, if F is a reduced rational first integral for (S), then P divides D(P ), and

Q divides D(Q) in K[x, y]. This motivates the following definition. Definition 2. A polynomial M ∈ K[x, y] \ K is a Darboux polynomial for D if M divides D(M ) in K[x, y]. Therefore, if M is a Darboux polynomial for D, then there exists a polynomial Λ ∈ K[x, y] such that D(M ) = Λ M . Such a polynomial Λ ∈ K[x, y] is called a cofactor associated with the Darboux polynomial M .
Darboux polynomials were introduced by G. Darboux in [Dar78]. These polynomials correspond to algebraic curves invariant under the vector field. The following lemma will be used in the sequel: it means that if we have a non-singular initial condition, then there is a unique irreducible invariant algebraic curve satisfying this initial condition, see [Sin92, Lemma A.1].

Lemma 3. Let D = A(x, y) ∂ ∂x + B(x, y) ∂ ∂y be the derivation attached to (S) and let (x 0 , y 0 ) be a non-singular point of D, i.e., A(x 0 , y 0 ) = 0 or B(x 0 , y 0 ) = 0. If M 1 and M 2 are two Darboux polynomials for D such that M 1 (x 0 , y 0 ) = M 2 (x 0 , y 0 ) = 0 and if M 1 is irreducible, then M 1 divides M 2 .

Darboux polynomials are sometimes called partial first integrals in the literature. The reason is that rational first integrals and Darboux polynomials are intimately related notions: as sketched above, numerators and denominators of reduced rational first integrals are Darboux polynomials. The converse is also true, see Corollary 5 below.

A fundamental property of Darboux polynomials is given in the following lemma (see, e.g., [DLA06, Lemma 8.3, p. 216]) that can be proved by a straightforward calculation.

Lemma 4. Let M ∈ K[x, y] and let M = M 1 M 2 be a factorization of M in K[x, y],
with M 1 and M 2 coprime. Then, M is a Darboux polynomial for D if and only if M 1 and M 2 are Darboux polynomials for D. Furthermore, if Λ M , Λ M1 and Λ M2 denote respectively the cofactors of M, M 1 and M 2 , then

Λ M = Λ M1 + Λ M2 .
As a corollary, one obtains: Corollary 5. Let F = P/Q be a reduced rational function in K(x, y). Then F is a rational first integral of (S) if and only if P and Q are Darboux polynomials for D with the same cofactor.

The previous corollary gives a relation between Darboux polynomials and rational first integrals. The next theorem shows that if we have enough Darboux polynomials, then we have a rational first integral, see [START_REF] Michael | Liouvillian first integrals of differential equations[END_REF]Appendix] for a modern proof. Theorem 6. Jou79,[START_REF] Michael | Liouvillian first integrals of differential equations[END_REF]] If d = max(deg(A), deg(B)), then the polynomial vector field (S) has a reduced rational first integral P/Q if and only if D has at least d (d + 1)/2 + 2 irreducible Darboux polynomials. In this case, D has infinitely many irreducible Darboux polynomials and any of them divides a linear combination λ Pµ Q, for some λ, µ ∈ K not both zero. Moreover, all but finitely many irreducible Darboux polynomials are of the form λ Pµ Q and have the same degree.

A useful corollary of Theorem 6 is the following, see [START_REF] Michael | Liouvillian first integrals of differential equations[END_REF]:

Corollary 7. For each planar polynomial vector field (S), there exists a nonnegative integer N (S) such that any irreducible Darboux polynomial for the derivation D attached to (S) has degree at most N (S) .

Given a derivation D, the problem of finding a bound for the degree of irreducible Darboux polynomials is known to be difficult: this is the so-called Poincaré problem. It has been deeply studied in the literature and many partial results exist ([Poi91, CLN91, Car94, Per02, Wal00, CG06, LY05] and others) though the question is not fully solved yet. The fact that the derivation D = n x ∂ ∂x + y ∂ ∂y with n ∈ N * admits xy n as an irreducible Darboux polynomial shows that a bound depending only on the degrees of the entries cannot exist: arithmetic conditions on the coefficients of D have to be taken into account as well.

Consequently, given a planar polynomial vector field (S), or equivalently a derivation D, two distinct problems occur when we want to compute rational first integrals:

(1) Find a bound on the degree of the numerator and denominator of a rational first integral, that is a bound on the degree of irreducible Darboux polynomials;

(2) (P N ): given a degree bound N ∈ N, either compute a rational first integral F ∈ K(x, y) \ K of (S) of total degree at most N , or prove that no such F exists.

Our aim is to give an efficient algorithm to handle the second problem (P N ).

In this article we suppose that N ≥ d + 1. This hypothesis is natural because if a derivation has a polynomial first integral of degree N , then one can show that

N ≥ d + 1, see [FL07, Theorem 6] or [Poi91].
2.2. Non-composite rational functions and their spectrum. We recall here the definition of composite rational functions and what is called the spectrum of a rational function. We then use these notions to describe the kernel of the derivation D and to give some of its properties.

Definition 8. A rational function F (x, y) ∈ K(x, y) is composite if it can be written F = u • h, i.e., F = u(h), where h ∈ K(x, y) \ K and u ∈ K(T ) with deg(u) ≥ 2.
Otherwise F is said to be non-composite.

In [START_REF] Moulin | Algebraic closure of a rational function[END_REF]Chè12b], the authors propose different algorithms for the decomposition of rational functions using properties of Darboux polynomials and rational first integrals of the Jacobian derivation. Lemma 9. The set of all rational first integrals of (S) is a K-algebra. It is closed under composition with rational functions in K(T ), and moreover, F is a rational first integral of (S) if and only if u • F is a rational first integral of (S) for all u ∈ K(T ) \ K.

Proof. The first assertion directly follows from the fact that the derivation D = A(x, y) ∂ ∂x + B(x, y) ∂ ∂y is K-linear and satisfies Leibniz's rule

D(F 1 F 2 ) = F 1 D(F 2 ) + D(F 1 ) F 2 .
The second assertion follows from the equality

D(u • F ) = u ′ (F ) D(F ),
and the fact that u ′ (F ) is zero if and only if u ∈ K.

A more precise version of Lemma 9 is given by the next theorem which completely describes the K-algebra structure of the set of all rational first integrals of (S). This theorem seems to be a folklore result but we have not found a suitable reference. Consequently, a complete proof is provided here.

Theorem 10. Let D be the derivation attached with (S). Then we have:

{G ∈ K(x, y) | D(G) = 0} = K(F ),
for some non-composite reduced rational first integral F of (S). Then any other rational first integral G of (S) is of the form G = u • F for some u ∈ K(T ) \ K. In particular, any two non-composite reduced rational first integrals are equal, up to a homography.

Proof. Let L = {G ∈ K(x, y) | D(G) = 0}. We have K ⊂ L ⊂ K(x, y), so, from [Sch00, Theorem 1, p. 12], we deduce that L is finitely generated over K and that L = K(f 1 , f 2 , f 3 ) for some f 1 , f 2 , f 3 ∈ K(x, y). As for i ∈ {1, 2, 3}, A(x, y) ∂fi ∂x + B(x, y) ∂fi ∂y = 0, we get that: ∂f i ∂y ∂f j ∂x -∂f i ∂x ∂f j ∂y = 0, for all i, j ∈ {1, 2, 3}.

The Jacobian criterion implies that f 1 , f 2 , f 3 are algebraically dependent and thus the transcendence degree of L over K is equal to one. By the extended Luroth's theorem, see [Sch00, Theorem 3, p. 15], we get L = K(F ), for F ∈ K(x, y). In particular, F is a rational first integral of (S). Now, if F is composite, F = u(H), with deg(u) ≥ 2, then K(F ) K(H), see, e.g., [GRS02, Proposition 2.2]. By Lemma 9, H is also a rational first integral of (S) so that H ∈ L and K(H) ⊂ L. This yields L = K(F ) K(H) ⊂ L, which is absurd. Thus F is non-composite which gives the desired result.

As a consequence of Theorem 10, non-composite reduced rational first integrals coincide with rational first integrals with minimal degree; they will play a key role in the remaining of this text.

In Definition 1, we have defined rational first integrals as elements of K(x, y). The next theorem explains why it is in general not necessary to consider rational first integrals in K(x, y). To our knowledge this result is proved here for the first time. In [MM97], the authors show that if there exists a rational first integral in K(x, y), then there also exists a rational first integral in K(x, y). We improve this result by taking into account the degrees of these rational first integrals.

Theorem 11. If (S) admits a non-composite rational first integral in K(x, y), then it admits a non-composite rational first integral in K(x, y) with the same degree.

Proof. Let f ∈ K(x, y) be a non-composite rational first integral of (S). Let G denote the Galois group over K of the smallest Galois extension containing all the coefficients of f . Given an indeterminate Z, we denote by N (f -Z) the product

N (f -Z) = σi∈G (σ i (f ) -Z) .
By construction, we have

N (f -Z) ∈ K(x, y)[Z]. Indeed, N (f -Z) is the minimal polynomial of f over K(x, y) so, as f is non-constant, so is N (f -Z).
It follows that, for almost any constant c ∈ K, N (fc) is non-constant; so, modulo a shift by some constant c, we may assume that N (f ) is non-constant 3 . As A, B ∈ K(x, y), it follows that N (f ) is a non-trivial rational first integral of (S) in K(x, y). Thus, by Lemma 9, there exists a non-composite rational first integral F ∈ K(x, y) of (S). Now, applying Theorem 10 with ground field K instead of K, we get that F = u(f ), with u ∈ K(T ). Furthermore, by [BCN11, Theorem 13], F is non-composite in K(x, y) implies that F is non-composite in K(x, y). It thus follows that deg(u) = 1 so that deg(F ) = deg(f ). Now, we introduce the spectrum of a rational function which will play a crucial role in our algorithms.

Definition 12. Let P/Q ∈ K(x, y) be a reduced rational function of degree N . The set

σ(P, Q) = {(λ : µ) ∈ P 1 K | λ P -µ Q is reducible in K[x, y], or deg(λ P -µ Q) < N }
is called the spectrum of P/Q.

In the context of rational first integrals of polynomial vector fields, the elements of the spectrum are sometimes called remarkable values, see, e.g., [START_REF] Ferragut | On the remarkable values of the rational first integrals of polynomial vector fields[END_REF]. There exists a vast bibliography about the spectrum, see for example [Rup86, Lor93, Vis93b, Vis93a, AHS03, Bod08, BC11].

The spectrum σ(P, Q) is finite if and only if P/Q is non-composite and if and only if the pencil of algebraic curves λ Pµ Q = 0 has an irreducible general element, see for instance [Jou79, Chapitre 2, Théorème 3.4.6] or [Bod08, Theorem 2.2] for detailed proofs.

To the authors' knowledge, the first effective result on the spectrum is due to Poincaré. In [Poi91], he establishes a relation between the number of saddle points and the cardinal of the associated spectrum, in the case where all the singular points of the polynomial vector field are distinct. In particular this yields the bound |σ(P, Q)| ≤ d 2 on the cardinality of the spectrum. This bound was improved recently in [Chè12a]:

3 We thank the referees for correcting an early mistake in this proof; this part of the proof is also found in [MM97].

Theorem 13. Let D be the derivation attached with (S) and d denotes the degree of (S). If P/Q is a reduced non-composite rational first integral of (S), then:

|σ(P, Q)| ≤ B(d) + 1, where B(d) = d(d + 1) 2 .
As a consequence of Theorem 13 and Corollary 5, if P/Q is a reduced noncomposite rational first integral of (S), then for all but B(d) + 1 constants σ ∈ K, the polynomial Pσ Q has degree N and is an irreducible Darboux polynomial for D with the same cofactor as P and Q. This means that if (S) has a rational first integral, there exist an infinite number of irreducible Darboux polynomials which all have the same degree (and the same cofactor).

The following lemma will be useful in Section 4 for the study of our probabilistic and deterministic algorithms.

Lemma 14. If P/Q ∈ K(x, y) is a reduced non-composite rational function of degree at most N , then the number of values of c ∈ K for which (Q(0, c) : P (0, c)) belongs to σ(P, Q) is bounded by N (B(d) + 1).

Proof. Let c ∈ K be such that (Q(0, c) : P (0, c)) ∈ σ(P, Q). By Theorem 13, σ(P, Q) contains at most B(d) + 1 elements. Now, for each (λ : µ) ∈ σ(P, Q), as P and Q are of degree at most N , there exist at most N values of c such that λ Q(0, c)µ P (0, c) = 0. This ends the proof.

Existing algorithms for computing rational first integrals and Darboux polynomials of bounded degree.

There is a vast literature regarding the computation of rational or elementary first integrals (see for example [FG10, MM97, Chè11, PS83, SGR90, Sin92, LZ10, DDdMS01, Poi91]) and Darboux polynomials (see for example [CMS09, CMS06, CGG05, Wei95, Dar78]). Note that, among these articles, very few restrict to the specific question of rational first integrals. Surveys on computing first integrals (not restricted to planar systems) can be found for example in [Gor01, Sch93, DLA06] and [Olv93] for symmetry methods which we do not address here.

Given a degree bound N , the naive approach to solve (P N ) consists in using the method of undetermined coefficients. This leads to a system of polynomial (quadratic) equations in the unknown coefficients of the rational first integral, see [Chè11] for a complexity estimate of this approach.

Interest in Darboux polynomials has been revived by the appearance of the Prelle-Singer's method, [START_REF] Prelle | Elementary first integrals of differential equations[END_REF]MM97,[START_REF] Duarte | Solving second-order ordinary differential equations by extending the Prelle-Singer method[END_REF]. In [CMS06, CMS09], Coutinho and Menasché Schechter give necessary conditions for the existence of Darboux polynomials. Other necessary conditions are contained in [CGGL03, CGG05] and also in works on inverse integrating factors [CGGL03, CFL10]. The bottleneck of the Prelle-Singer's method and all of its variants is the computation by undetermined coefficients of all irreducible Darboux polynomials of bounded degree, which leads again to solving a polynomial system. This yields an exponential complexity algorithm, see [Chè11].

In [Chè11], Chèze shows that if the derivation D admits only finitely many irreducible Darboux polynomials of degree at most N , then it is possible to compute all of them by using the so-called ecstatic curve from [Per01] within a number of binary operations that is polynomial in the bound N , in the degree d = max(deg(A), deg(B)) of D and in the logarithm of the height of A and B. A nontrivial modification of this algorithm provides a polynomial-time method to solve (P N ), see again [Chè11]. To our knowledge, this is the first algorithm solving (P N ) in polynomial-time. Unfortunately, the exponent is quite large, making the algorithm unpractical even for moderate values of N . This drawback is due to the fact that algorithm Rat-First-Int in [Chè11] needs to compute the irreducible factors of a bivariate polynomial of degree ≈ dN 4 , and the best known algorithms for solving this subtask have arithmetic complexity, e.g., Õ(

d ω+1 N 4ω+4 ), see [BLS + 04, Lec06].
Last, we mention the article [FG10] of Ferragut and Giacomini, which is the starting point of this work. They propose two methods to compute rational first integrals, together with several useful local filters. The key idea is to try to recover a rational first integral H = N i,j=0 fi,j x i y j N i,j=0 gi,j x i y j from the minimal polynomial F N (x, y) := N i+j=0 (f i,j + cg i,j )x i y j of an algebraic power series solution of the differential equation (E) (see next section for justifications of this). Their first method (section 2, page 93) consists in computing a power series solution with generic initial condition y(0) = a 0 , finding c (from equation F N (0, a 0 ) = 0), and then computing the f i,j and g i,j from relation F N (x, y(x)) = 0. This yields equations that are quadratic in the f 0,j , g 0,j and linear in the other variables. The second method in [START_REF] Ferragut | A new algorithm for finding rational first integrals of polynomial vector fields[END_REF] (Remark 1) uses Puiseux series at a resonant node to obtain linear equations (though, strangely, the authors do not mention this). Both methods are designed for hand-driven computer-aided calculations and are remarkably elegant. However, if we ask for an algorithm which either computes a rational first integral of degree ≤ N or proves that no such first integral exists, then several ingredients are missing. For example, a bound is necessary on the number of terms needed in a power series to certify that there is no mistake in the minimal polynomial (we give such tight bounds in Lemma 18). One also needs to specify whether intermediate algebraic extensions will be needed. This is potentially the case in their second method. We prove it to be unnecessary in Theorem 11 and use power series at regular points to ensure that no such undesired extension appears. Our initial goal was to transform their methods into implementable algorithms. While doing this, it turned out that we could design different algorithms that were not only implementable (and implemented) but also much faster than their methods (see Section 6 for timings and comparisons). Our first strategy is to bypass their Step 4: we rather compute the minimal polynomial of a generic power series solution y(x) directly and then show how to obtain H as a quotient of two specializations of this generic object. This simple idea seems to be new.

This strategy allows to design a much faster probabilistic algorithm (based on the fact that, in the above process, almost any specialization works -with a bound on this "almost" in Lemma 14 and second part of Proposition 16). We also get a fast deterministic algorithm based on this "good specialization" property. We study carefully their (polynomial) complexities. Our implementations show that our algorithms are more efficient both in theory and in practice than all previous algorithms, see Section 6.

3. Rational first integrals, differential equations and algebraic power series 3.1. Algebraic power series and rational first integrals.

Definition 15. A formal power series y(x

) ∈ K[[x]] is said to be algebraic if it is algebraic over K(x), that is, if there exists a non-zero polynomial M ∈ K[x, y] such that M (x, y(x)) = 0. An irreducible polynomial M ∈ K[x, y] satisfying M (x, y(x)) = 0 is called a minimal polynomial of y(x) in K[x, y].
With the planar polynomial vector field (S), we associate the first order nonlinear differential equation:

(E) dy dx = B(x, y) A(x, y)
We may assume without any loss of generality that x does not divide A, i.e., A(0, y) ≡ 0. We will explain how we can reduce to this situation and study the complexity of this reduction in Subsection 5.1.

Then, the formal version of the Cauchy-Lipschitz theorem for non-linear (firstorder) differential equations ensures that for any c ∈ K such that A(0, c) = 0, the equation (E) admits a unique power series solution y c (x) ∈ K[[x]] satisfying y c (0) = c. Note that high-order truncations of the power series y c (x) can be computed efficiently using the algorithm of Brent and Kung [START_REF] Brent | Fast algorithms for manipulating formal power series[END_REF].

The following standard result is fundamental to both our method and the one of Ferragut and Giacomini in [START_REF] Ferragut | A new algorithm for finding rational first integrals of polynomial vector fields[END_REF].

Proposition 16. Consider the planar polynomial vector field (S) and assume that A(0, y) ≡ 0. Let c ∈ K satisfy A(0, c) = 0 and y c (x) ∈ K[[x]] be the unique power series solution of (E) such that y c (0) = c.

(1) If (S) admits a non-composite rational first integral P/Q, then the power series y c (x) is algebraic. More precisely, y c (x) is a root of the non-zero polynomial λ Pµ Q, where λ = Q(0, c) and µ = P (0, c).

(2) If P/Q is a reduced non-composite rational first integral of (S) of degree at most N , then, for all but at most N (B(d) + 1) values of c ∈ K, the polynomial λ Pµ Q, where λ = Q(0, c) and µ = P (0, c) is a minimal polynomial of y c (x).

Proof. Let F = P/Q be a non-composite rational first integral of (S). Since the spectrum σ(P, Q) is finite, we can suppose that P and Q are irreducible and coprime. Let us first show that P (0, c) = 0 or Q(0, c) = 0. As (0, c) is a non-singular point of D, if P (0, c) = Q(0, c) = 0, then Lemma 3 implies that P = α Q with α ∈ K.

As P and Q are coprime, we get a contradiction so that necessarily P (0, c) = 0 or Q(0, c) = 0. We thus suppose Q(0, c) = 0 else we consider the rational first integral Q/P instead of P/Q. As Q(0, c) = 0, the power series Q(x, y c (x)) is invertible so that D(F ) = 0 yields D(F )(x, y c (x)) = 0. The latter equality can be written

A(x, y c (x)) ∂F ∂x (x, y c (x)) + B(x, y c (x)) ∂F ∂y (x, y c (x)) = 0.
Dividing the equality by the invertible power series A(x, y c (x)) and using the fact that y c (x) is a solution of (E), we obtain

∂F ∂x (x, y c (x)) + B(x, y c (x)) A(x, y c (x)) ∂F ∂y (x, y c (x)) = 0, ∂F ∂x (x, y c (x)) + dy c (x) dx ∂F ∂y (x, y c (x)) = 0, d F (x, y c (x)) dx = 0.
It follows that F (x, y c (x)) = σ c , for some σ c ∈ K. Consequently, we have

P (x, y c (x)) -σ c Q(x, y c (x)) = 0,
with σ c = P (0, c)/Q(0, c) which proves (1). If, in addition, F = P/Q is reduced non-composite of degree at most N , then (2) follows directly from Lemma 14.

Proposition 16 shows in particular that if (S) has a rational first integral P/Q, then all power series solutions of (E) are algebraic. The next proposition which is well known (see [START_REF] Ferragut | A new algorithm for finding rational first integrals of polynomial vector fields[END_REF]Wei95]) asserts that the converse is also true.

Proposition 17. Let (S) be a planar polynomial vector field, D the associated derivation, and (E) be the associated differential equation.

(1) If M ∈ K[x, y] is an irreducible Darboux polynomial for D, then all roots y(x) ∈ K(x) of M such that A(0, y(0)) = 0 are power series solutions of (E).

(2) The minimal polynomial of an algebraic solution y(x) ∈ K[[x]] of (E) such that A(0, y(0)) = 0 is a Darboux polynomial for D.

(3) (S) admits a rational first integral if and only if all the power series solutions of (E) are algebraic.

Proof. Assume first that M ∈ K[x, y] is an irreducible Darboux polynomial for D, and that y(x) ∈ K(x) is a root of M , i.e., M (x, y(x)) = 0. Since M divides D(M ), it follows that y(x) is also a root of D(M ). This implies

∂M ∂x (x, y(x)) = - B(x, y(x)) A(x, y(x)) ∂M ∂y (x, y(x)).
On the other hand, M (x, y(x)) = 0 implies by differentiation with respect to x that

∂M ∂x (x, y(x)) = - dy(x) dx ∂M ∂y (x, y(x)).
These two equalities provides

∂M ∂y (x, y(x)) dy(x) dx - B(x, y(x)) A(x, y(x)) = 0.
As M is irreducible, ∂M ∂y (x, y(x)) = 0 so that y(x) is a solution of (E). As mentioned before A(0, y(0)) = 0 implies y(x) ∈ K[[x]] (Cauchy-Lipschitz theorem) which proves (1).

Assume now that y(x) ∈ K[[x]

] is an algebraic solution of (E) such that A(0, y(0)) = 0, and let M be its minimal polynomial. Then

0 = d M (x, y(x)) dx = ∂M ∂x (x, y(x)) + dy(x) dx ∂M ∂y (x, y(x)),
and since y(x) is a solution of (E), the latter equality implies that y(x) is also a root of D(M ). Now as M is the minimal polynomial of y(x), it follows that M divides D(M ), i.e., M is a Darboux polynomial for D and we have proved (2).

Let us now prove (3). If all the power series solutions of (E) are algebraic, then by (2), D admits infinitely many Darboux polynomials. Then Theorem 6 shows that (S) admits a rational first integral. The proof ends here since the other implication of (3) has been proved in Proposition 16.

3.2. Algebraic power series solutions of (E). We have seen in Proposition 16 that if P/Q is a reduced non-composite rational first integral of (S), then a minimal polynomial of a power series solution of (E) is generically of the form λ Pµ Q for some constants λ and µ. Thus, if we are able to compute such a minimal polynomial, we can deduce the rational first integral P/Q. In practice, we do not compute a power series y c (x) ∈ K[[x]] solution of (E) but only a truncation of y c (x), i.e., a finite number of terms of its expansion in the monomial basis. Given a degree bound N for the rational first integral that we are searching for, the following lemma shows that computing y c (x) mod x N 2 +1 , i.e., the first N 2 + 1 terms of its expansion, is enough for our purposes.

Such an analysis of the needed precision for the power series solutions of (E) that we compute is not included in [START_REF] Ferragut | A new algorithm for finding rational first integrals of polynomial vector fields[END_REF]. Note that this kind of strategy was already used in a polynomial factorization setting (see, e.g., [START_REF] Kaltofen | Polynomial-time reductions from multivariate to bi-and univariate integral polynomial factorization[END_REF]) and in a differential equations setting (see [ACFG05, BCS + 07]). The next lemma is a small improvement of [ACFG05, Lemma 2.4].

Lemma 18. Let L be a field of characteristic 0 such that K ⊆ L. Let ŷ(x) ∈ L[[x]] denote an algebraic power series whose minimal polynomial M ∈ L[x, y] has degree at most N . If M ∈ L[x, y] is a polynomial of degree at most N satisfying (⋆) : M (x, ŷ(x)) ≡ 0 mod x N 2 +1 , then M (x, ŷ(x)) = 0. Moreover, if M has minimal degree in y among polynomials satisfying (⋆), then M = f M for some f ∈ L[x].
Proof. By definition, M satisfies (⋆) so there exists M ∈ L[x, y] of degree at most N satisfying (⋆). Let M be such a solution of (⋆) and consider R(x) := Res y (M (x, y), M (x, y)), the resultant of M and M with respect to y. As there exist polynomials S and T

in K[x, y] such that S M + T M = R, Relation (⋆) yields R(x) ≡ 0 mod x N 2 +1 . By Bézout's theorem, we have deg(R) ≤ deg(M ) deg( M ) ≤ N 2 , thus R = 0.
This implies that M and M have a non-trivial common factor. Now, as M is irreducible, necessarily M divides M and thus M (x, ŷ(x)) = 0. Finally, if M is supposed to have minimal degree in y among polynomials satisfying (⋆), we have necessarily M = f M for some f ∈ L[x] which ends the proof.

Note that in (⋆), the power series ŷ(x) can be replaced by its truncation ŷ(x) mod x N 2 +1 .

Remark 19. For a given power series ŷ(x), computing all the polynomials M of degree at most N satisfying (⋆) can be done by taking an ansatz for M and performing linear algebra calculations (e.g., solving a system of linear equations). Consequently, computing "all" solutions of (⋆) means computing "a basis" of solutions of the linear algebra problem associated with (⋆). An efficient method to address this problem and to get, via a row-echelon form, a solution of (⋆) with minimal degree in y is given in Subsection 5.3 where a complexity analysis is provided.

In the sequel, we say that M is a minimal solution of (⋆) if it is a solution of (⋆) with minimal degree in y.

3.3.

A first algorithm for computing rational first integrals. We now propose a first algorithm, based on linear algebra, for solving (P N ). More efficient algorithms, based on this one, are given in Section 4. The strategy of this algorithm is then used in Section 7 for computing Darboux polynomials.

Algorithm GenericRationalFirstIntegral

Input:A, B ∈ K[x, y] s.t. A(0, y) ≡ 0 and a bound N ∈ N. Output: A non-composite rational first integral of (S) of degree at most N , or "None".

(1) For an indeterminate c, compute the polynomial y c ∈ K(c)[x] of degree at most (N 2 + 1) s.t. y c (0) = c and dyc dx ≡ B(x,yc) A(x,yc) mod x N 2 +1 .

(2) Compute all4 non-trivial polynomials M ∈ K(c)[x, y] of degree ≤ N s.t.

(⋆) : M (c, x, y c (x)) ≡ 0 mod x N 2 +1 .

If no such M exists, then Return "None". Else, among the solutions of (⋆), pick a minimal solution

M ∈ K[c][x, y].
(3) Let M denote the primitive part of M relatively to y. Set P (x, y) := M (0, x, y).

Pick any c 1 ∈ K s.t. M(c1,x,y) P (x,y) ∈ K and set Q(x, y) := M (c 1 , x, y).

(4) If D(P/Q) = 0, then Return P/Q. Else Return "None".

In the above algorithm, the output "None" means that there is no rational first integral of degree at most N but there may exist a rational first integral of degree strictly greater than N .

Theorem 20. Algorithm GenericRationalFirstIntegral is correct: either it finds a non-composite rational first integral of (S) of degree at most N if it exists, or it proves that no such rational first integral exists.

To prove Theorem 20, we shall need the following lemma.

Lemma 21. Consider the planar polynomial vector field (S) and assume that A(0, y) ≡ 0. If F is a reduced rational first integral of (S), then F (0, y) ∈ K(y) \ K.

Proof. Let F = P/Q be a reduced rational first integral of (S). Proceeding by contradiction, we assume F (0, y) = c 0 ∈ K. Then P (0, y)c 0 Q(0, y) = 0 so that x divides P (x, y)c 0 Q(x, y). Now, from Corollary 5, P (x, y)c 0 Q(x, y) is a Darboux polynomial for D and thus, by Lemma 4, x is also a Darboux polynomial for D. Consequently, we get that x divides A(x, y) and thus A(0, y) = 0. This is absurd so we conclude F (0, y) ∈ K(y) \ K.

Proof of Theorem 20. Suppose first that there exists a rational first integral of (S) of degree at most N . Then, without loss of generality, we can consider a reduced non-composite one P 0 /Q 0 , see Lemma 9. Let y c ∈ K(c) [[x]] be the power series solution of (E) satisfying y c (0) = c. By Proposition 16, y c is a root of λ P 0µ Q 0 , where λ = Q 0 (0, c) and µ = P 0 (0, c).

As P 0 /Q 0 is non-composite, it follows that λ P 0 -µ Q 0 is irreducible in K(c)[x, y]
. Indeed, by Lemma 21, the constant µ/λ belongs to K(c) \ K and thus, from Lemma 14, we can find c 0 ∈ K such that (Q 0 (0, c 0 ) : P 0 (0, c 0 )) ∈ σ(P 0 , Q 0 ). Consequently λ P 0µ Q 0 is a minimal polynomial of y c . In Step (1), we compute the first N 2 + 1 terms of y c . Now, in Step (2), if there exists a solution M ∈ K(c)[x, y] of (⋆) of degree at most N , then, Lemma 18 applied with

L = K(c) implies that M = f (λ P 0 -µ Q 0 ) with f ∈ K(c)[x], where M ∈ K[c][x, y] is defined in Step (2)
. Therefore, taking the primitive part of M with respect to y, in Step (3), we have M = g (λ P 0µ Q 0 ) for some g ∈ K[c]. Now, if P and Q denote the polynomials defined in Step (3) of the algorithm, we necessarily have:

P Q = α P 0 + β Q 0 δ P 0 + γ Q 0 ∈ K(x, y) \ K,
where α, β, δ, γ ∈ K.

As P 0 /Q 0 is a non-composite rational first integral, we deduce that P/Q is also a non-composite rational first integral. Thus, we have D(P/Q) = 0 in Step (4) and the algorithm returns a correct output. Now suppose that (S) has no rational first integral of degree at most N . In Step (4), the test D(P/Q) = 0 guarantees to return a correct output. In Step (2), we can have an early detection of this situation. Indeed by Proposition 16, if (⋆) has no non-trivial solution, then we deduce that (S) has no rational first integral of degree at most N .

This algorithm fits the first part of our goal as it is entirely based on linear operations: we do not need to solve quadratic equations (see Section 5). However, it is not yet very efficient in practice because computations are done over K(c). For example, in the first step, a direct calculation shows that, for n ≥ 1, the coefficient of x n in the power series solution y c of (E) satisfying y c (0) = c is generically a rational function in c of degree (2n -1) d, whose denominator is generically A(0, c) 2 n-1 . In what follows, we accelerate things by using only computations over K instead of computations in K(c).

Efficient algorithms for computing rational first integrals

4.1. A probabilistic algorithm. In this section, we present an efficient probabilistic algorithm of Las Vegas type for solving (P N ). The approach is similar to the one used in the previous section.

Algorithm ProbabilisticRationalFirstIntegral

Input:A, B ∈ K[x, y] s.t. A(0, y) ≡ 0, two elements c 1 , c 2 ∈ K s.t. A(0, c i ) = 0 for i = 1, 2, and a bound N ∈ N. Output: A non-composite rational first integral of (S) of degree at most N , "None" or "I don't know".

(1) For i = 1, 2 do: (1a) Compute y ci ∈ K[x] of degree at most (N 2 + 1) s.t. y ci (0) = c i , and

dyc i dx ≡ B(x,yc i ) A(x,yc i ) mod x N 2 +1 . (1b) Compute all non-trivial polynomials Mi ∈ K[x, y] of degree ≤ N s.t. (⋆) : Mi (x, y ci (x)) ≡ 0 mod x N 2 +1 .
(1c) If no such Mi exists, then Return "None".

Else let M i ∈ K[x, y] be the primitive part relatively to y of a minimal solution of (⋆).

(

) If i = 1, then while (M 1 (0, c 2 ) = 0 or A(0, c 2 ) = 0) do c 2 = c 2 + 1. (2) If D(M 1 /M 2 ) = 0, then Return M 1 /M 2 . Else Return ["I don't know",[c 2 ]]. 1d
Theorem 22. Algorithm ProbabilisticRationalFirstIntegral terminates and satisfies the following properties:

• If it returns M 1 /M 2 , then it is a non-composite rational first integral of (S) of degree at most N . • If it returns "None", then there is no rational first integral of (S) of degree at most N . • If (S) admits a non-composite rational first integral P/Q of degree at most N and (Q(0, c i ) : P (0, c i )) ∈ σ(P, Q) for i = 1, 2, then the algorithm returns a non-composite rational first integral of (S) of degree at most N .

Proof. Let us first prove that the algorithm terminates. This follows directly from the fact that the while loop in Step (1d) terminates after at most N + d + 1 steps. Indeed, we just have to avoid the roots of the product M 1 (0, y) A(0, y) which a univariate polynomial of degree less than N + d. It thus remains to check that it is a non-zero polynomial, i.e., M 1 (0, y) ≡ 0. If M 1 (0, y) ≡ 0, then x divides M 1 .

As M 1 is the primitive part with respect to y of a minimal solution of (⋆), this would imply that M 1 (x, y) = x and thus M 1 (x, y c1 (x)) ≡ 0 mod x N 2 +1 which is a contradiction. Now, if the algorithm returns M 1 /M 2 , then the test in Step (2) ensures that D(M 1 /M 2 ) = 0 and, by construction, M 1 /M 2 is clearly of degree at most N . Furthermore, M 1 /M 2 is non-composite. Indeed, if M 1 /M 2 is composite, then at least one of the M i 's is reducible and thus it can not be the primitive part with respect to y of a minimal solution of (⋆). Finally Step (1d) certifies that M 1 /M 2 ∈ K. Indeed, M 2 satisfies M 2 (0, c 2 ) = 0, thus if M 2 = k M 1 with k ∈ K, then either k = 0 or M 1 (0, c 2 ) = 0 which is not possible thanks to Step (1d). We have then proved that M 1 /M 2 is a non-composite rational first integral of (S) of degree at most N . If the algorithm returns "None" in Step (1c), then by Proposition 16, (S) has no rational first integral of degree at most N . Assume finally that (S) admits a non-composite rational first integral P/Q of degree at most N and that (Q(0, c i ) : P (0, c i )) ∈ σ(P, Q) for i = 1, 2. Then the same strategy as the one used in the proof of Theorem 20 shows that our algorithm returns a non-composite rational first integral of (S).

Proposition 23. Let Ω be a (finite) subset of K of cardinal |Ω| greater than N (B(d)+1) and assume that, in Algorithm ProbabilisticRationalFirstIntegral, c 1 and c 2 are chosen independently and uniformly at random in Ω. Then, if (S) admits a rational first integral of degree at most N , Algorithm ProbabilisticRationalFirstIntegral returns a non-composite rational first integral of (S) of degree at most N with probability at least 1 -N (B(d)+1)

|Ω|

.

Proof. It is a straightforward application of Lemma 14, Theorem 22 and Zippel-Schwartz's lemma (see [gGG99, Lemma 6.44]).

In fact, the "practical" probability will be much better. Indeed, the elements (λ : µ) of the spectrum may be rational or algebraic and hence, the constants c such that (Q(0, c) : P (0, c)) ∈ σ(P, Q) will generally be algebraic. So, if the c i 's are chosen to be rational in the input, then the "bad" values of the c i 's will generally be in very small number. This fact is widely confirmed by experiments. Now, we study all the different situations that can occur and the corresponding output given by the algorithm ProbabilisticRationalFirstIntegral:

(1) (S) has a non-composite rational first integral P/Q of degree at most N . (a) If (Q(0, c 1 ) : P (0, c 1 )) ∈ σ(P, Q), and (Q(0, c 2 ) : P (0, c 2 )) ∈ σ(P, Q) then in this situation the algorithm returns a non-composite rational first integral. (b) Now, we study the opposite situation: (Q(0, c 1 ) : P (0, c 1 )) ∈ σ(P, Q)

or (Q(0, c 2 ) : P (0, c 2 )) ∈ σ(P, Q). If the algorithm computes M 1 and M 2 but M 1 /M 2 is not a rational first integral, then it returns "I don't know". A first example where this case is encountered is given in Subsection 6.3. Furthermore, we may be unlucky enough to choose two bad values of the c i 's, i.e., (Q(0, c i ) :

P (0, c i )) ∈ σ(P, Q) for i = 1, 2.
For example if we consider

A(x, y) = -4 x 3 + 4 xy 2 + 6 x 2 -2 y 2 -2 x, B(x, y) = -4 x 2 y + 4 y 3 + 4 x y -2 y,
then (S) has a non-composite rational first integral P/Q of degree 2, where P (x, y) = (y -x) (y -x+1) and Q(x, y) = (y +x) (y +x-1). But if we choose c 1 = -1 and c 2 = 1, then we will construct two Darboux polynomials M 1 (x, y) = yx + 1 and M 2 (x, y) = y + x -1 of degree only 1 that are minimal polynomials of y c1 and y c2 . As deg(M 1 /M 2 ) is strictly smaller than deg(P/Q), we obtain D(M 1 /M 2 ) = 0 and the algorithm returns "I don't know". (2) (S) does not have a rational first integral with degree at most N .

(a) If (⋆) has no non-trivial solutions, then the algorithm returns "None".

(b) If (⋆) has non-trivial solutions, then the algorithm returns "I don't know". This situation can occur for example when:

• (S) has no rational first integral but it has Darboux polynomials and the choice of c 1 and c 2 gives two Darboux polynomials. For an example of a derivation without rational first integral but with Darboux polynomials, see [Chè11, Remark 15]. • (S) has a rational first integral with degree bigger than the given bound N . For example, consider the derivation D = (x + 1) ∂ ∂xy ∂ ∂y and the degree bound N = 1. In this situation, the differential equation is (E) : dy dx = -y x+1 which admits y c (x) = c 1+x as solution. We set M (x, y) = α + βx + γy, and then M (x, y c (x)) = 0 mod x 2 gives M (x, y) = γ(-c + cx + y). However, M (x, y c (x)) = x 2 mod x 3 , thus y c (x) is not a root of M . Here D admits the rational first integral y (x + 1) so if we set N = 2 in the input, our algorithm returns a non-composite rational first integral of degree 2. In this case we compute y c (x) mod x 5 .

Remark 24. The bivariate polynomials Mi 's computed in Step (1b) have total degree at most N so they have (N + 1)(N + 2)/2 coefficients. Note that, if we assume N ≥ 3, then we have N 2 + 1 ≥ (N + 1)(N + 2)/2. It is tempting to try to compute the Mi 's using only, say, (N + 1)(N + 2)/2 + 2 terms of the power series. This will make the computation a little bit faster, but then the method becomes only a nice heuristic and may fail.

A deterministic algorithm. Algorithm

ProbabilisticRationalFirstIntegral is now turned into a deterministic algorithm. The idea is that if a rational first integral with degree at most N exists then, if we run at most N (B(d) + 1) + 1 times ProbabilisticRationalFirstIntegral, we will get a non-composite rational first integral of degree at most N .

Algorithm DeterministicRationalFirstIntegral

Input:A, B ∈ K[x, y] s.t. A(0, y) ≡ 0 and a bound N ∈ N. Output: A non-composite rational first integral of (S) of degree ≤ N or "None".

(1) Let Ω := ∅. Step (2). (e) Else Return F .

(3) Return "None".

Theorem 25. Algorithm DeterministicRationalFirstIntegral is correct: it returns a rational first integral of degree at most N if and only if it exists, and it returns "None" if and only if there is no rational first integral of degree at most N .

Proof. Assume that (S) has a non-composite rational first integral P/Q with degree at most N . If F ="I don't know" in Step (2), then from Theorem 22, at least one of the c i 's satisfies (Q(0, c i ) : P (0, c i )) ∈ σ(P, Q). The number of such "bad" values of the c i 's is bounded by N (B(d) + 1) by Lemma 14. Hence if we repeat ProbabilisticRationalFirstIntegral at least N (B(d) + 1) + 1 times, then we will get a good pair (c 1 , c 2 ) and by Theorem 22, the probabilistic algorithm will then return a non-composite rational first integral of degree at most N . Now assume that (S) has no rational first integral of degree at most N . Then by Theorem 22, ProbabilisticRationalFirstIntegral returns "None" or "I don't know". If in Step (2), F ="None", then we have a correct output. Now if F ="I don't know", then the algorithm uses again ProbabilisticRationalFirstIntegral with new values of the c i 's and, after at most N (B(d) + 1) + 1 trials, it returns "None" which is the correct output.

It is possible to improve the algorithm DeterministicRationalFirstIntegral in the following way: if the probabilistic algorithm returns "I don't know", then we test if we have computed Darboux polynomials with degree strictly smaller than N . In this case, we keep these polynomials and before running again ProbabilisticRa-tionalFirstIntegral we check if we can get a rational first integral with the help of all the Darboux polynomials computed so far. This can be done easily because it just requires testing if we can construct a cofactor equal to zero which, from Lemma 4, reduces to a linear algebra problem. This approach can sometimes decrease the number of uses of ProbabilisticRational-FirstIntegral in our deterministic algorithm. However, the worst case complexity analysis of these two approaches is the same. Indeed, we have a Darboux polynomial with degree smaller than N when (Q(0, c i ); P (0, c i )) ∈ σ(P, Q) and we can certify the existence of a rational first integral, thanks to Darboux-Jouanolou's theorem when we have B(d) + 2 irreducible Darboux polynomials (see Theorem 6). Thus, it yields (N B(d) + 1) + 1 bad values for c i as before.

Complexity analysis and algorithmic issues

In this section, we describe how the different steps of algorithms ProbabilisticRa-tionalFirstIntegral and DeterministicRationalFirstIntegral can be performed efficiently and we study their arithmetic complexities.

For the complexity issues, we focus on the dependency on the degree bound N and we recall that we assume that N ≥ d, where d = max(deg(A), deg(B)) denotes the degree of the polynomial vector field. More precisely, we suppose that d is fixed and N tends to infinity.

All the complexity estimates are given in terms of arithmetic operations in K. We use the notation f ∈ Õ(g): roughly speaking, it means that f is in O(g log m (g)) for some m ≥ 1. For a precise definition, see [START_REF] Von | Modern computer algebra[END_REF]Definition 25.8]. We suppose that the Fast Fourier Transform can be used so that two univariate polynomials with coefficients in K and degree bounded by D can be multiplied in Õ(D), see [START_REF] Von | Modern computer algebra[END_REF]. We further assume that two matrices of size n with entries in K can be multiplied using O(n ω ), where 2 ≤ ω ≤ 3 is the matrix multiplication exponent, see [START_REF] Von | Modern computer algebra[END_REF]Ch. 12]. We also recall that a basis of solutions of a linear system composed of m equations and n ≤ m unknowns over K can be computed using O(m n ω-1 ) operations in K, see [BP94, Chapter 2]. 5.1. Computation of a regular point. In the algorithms given in the previous sections, we have to choose a regular point for the differential equation (E), i.e., a point x 0 satisfying A(x 0 , y) ≡ 0. To achieve this, we can start from the point x 0 = 0, evaluate A(x, y) at x = x 0 . If A(x 0 , y) ≡ 0, then we are done. Else, we shift x 0 by one to get x 0 = 1 and we iterate the process. Note that the number of iterations is at most d. Consequently, this step can be performed by evaluating d polynomials (namely the coefficients of A(x, y) viewed as polynomials in the variable y) of degree bounded by d at d points (x 0 = 0, 1, 2, . . . , d -1). This can thus be done in Õ(d 2 ) arithmetic operations, see [START_REF] Von | Modern computer algebra[END_REF]Corollary 10.8]. This is why, in our algorithms, we always suppose, at neglectible cost and without loss of generality, that A(0, y) ≡ 0.

Power series solutions of (E). In

Step (1) of the algorithm ProbabilisticRa-tionalFirstIntegral, we compute the N 2 + 1 first terms of the power series solution of (E) satisfying a given initial condition. Using the result of Brent and Kung (see [START_REF] Brent | Fast algorithms for manipulating formal power series[END_REF]Theorem 5.1]) based on formal Newton iteration, this can be done using Õ(d N 2 ) arithmetic operations, see also [BCO + 07].

5.3.

Guessing the minimal polynomial of an algebraic power series. We shall now give a method for solving Problem (⋆) in Step (1b) of Algorithm Prob-abilisticRationalFirstIntegral. The problem is the following: given the first N 2 + 1 terms of a power series ŷ(x), find (if it exists), a bivariate polynomial M ∈ K[x, y], with minimal degree in y, such that M (x, ŷ(x)) ≡ 0 mod x N 2 +1 . This can be handled by an undetermined coefficients approach as follows: Algorithm GuessMinimalPolynomial

Input: A polynomial ŷ ∈ K[x] s.t. deg(ŷ) ≤ (N 2 + 1), with N ∈ N.
Output: A minimal solution of (⋆) with degree ≤ N or "None".

(1) Let M (x, y) = N i=0 N -i j=0 m i,j x j y i be an ansatz for the bivariate polynomial that we are searching for.

(2) Construct the linear system (L) for the m i,j 's given by:

M (x, ŷ(x)) = N i=0   N -i j=0 m i,j x j   ŷ(x) i ≡ 0 mod x N 2 +1 .
(3) If (L) does not have a non-trivial solution, then Return "None".

(4) Else compute a row-echelon form of a basis of solutions of (L) to find a solution M (x, y) of minimal degree in y and Return it.

Proposition 26. Algorithm GuessMinimalPolynomial is correct. If we suppose that N ≥ 3, then it uses at most Õ(N 2ω ) arithmetic operations in K.

Proof. The correctness of the algorithm is straightforward. Let us study its arithmetic complexity. To construct the linear system (L), we have to compute ŷ(x) i mod x N 2 +1 for i = 0, . . . , N . This can be done in Õ(N 3 ) arithmetic operations.

The linear system (L) has N 2 + 1 equations and (N + 1) (N + 2)/2 = O(N 2 ) unknowns m i,j 's. Note that we assume N ≥ 3 so that N 2 + 1 ≥ (N + 1) (N + 2)/2. It can thus be solved using O(N 2 (N 2 ) ω-1 ) operations. Finally, in Step (4), the row-echelon form can be computed using at most Õ(N 2 ω ) arithmetic operations (see [BP94, Chapter 3]) since the dimension of a basis of solutions of (L) does not exceed O(N 2 ), which ends the proof.

5.4. Total cost of our algorithms.

Theorem 27. Algorithm ProbabilisticRationalFirstIntegral uses at most Õ(N 2ω ) arithmetic operations in K, when N tends to infinity and d is fixed.

Proof. In Subsection 5.2, we have seen that Step (1a) can be performed in at most Õ(d N 2 ) arithmetic operations. Then, using Algorithm GuessMinimalPolynomial,

Step (1b) can be performed in Õ(N 2 ω ) operations in K, see Subsection 5.3. In

Step (1c), we have to compute the primitive part relatively to y of a minimal solution of (⋆). This reduces to computing N gcd's of univariate polynomials of degree at most N which can be done in O(N 3 ) operations in K (and even faster using half-gcd techniques). In Step (1d), we must avoid the roots of M 1 (0, y) A(0, y) thus we need to run the while loop at most d + N + 1 times. In this loop we evaluate univariate polynomials with degree at most d and N , thus it uses at most Õ((d + N ) 2 ) arithmetic operations. Finally, we test if D(M 1 /M 2 ) = 0 which costs Õ((d + N ) 2 ) arithmetic operations since N ≥ d. Indeed, we multiply bivariate polynomials of degree at most N and we add bivariate polynomials of degree at most d + 2N -1.

Corollary 28. The deterministic algorithm DeterministicRationalFirstIntegral can be done using at most Õ(d 2 N 2ω+1 ) arithmetic operations, when N tends to infinity and d is fixed.

In the previous statement, even if d is fixed, we mention it in the complexity in order to emphasize on the number of iterations of the probabilistic algorithm.

Proof. This estimate is straightforward from Theorem 27 since Algorithm Deter-ministicRationalFirstIntegral calls at most N (B(d) + 1) + 1 times the algorithm Prob-abilisticRationalFirstIntegral.

5.5. Faster heuristic using Padé-Hermite approximation. The algorithm GuessMinimalPolynomial developed in Subsection 5.3 uses an undetermined coefficients method to compute a minimal solution of (⋆) in Step (1b) of Algorithm ProbabilisticRationalFirstIntegral. It consists in finding (if it exists) the minimal polynomial of a power series. In the present section, we propose another approach to solve that problem using Padé-Hermite approximation, see [START_REF] Beckermann | A uniform approach for the fast computation of matrix-type Padé approximants[END_REF].

Indeed, the problem of computing a bivariate polynomial annihilating a power series can be handled by means of computing a Padé-Hermite approximant, see [Sha74, Sha78]. More precisely, given a power series ŷ(x), if there exists a bivariate polynomial M of degree N such that M (x, ŷ(x)) = 0, then the coefficients of the powers of y are a Padé-Hermite approximant of type (N, N -1, . . . , 0) of the vector of power series (1, ŷ(x), . . . , ŷ(x) N ) T . Computing such a Padé-Hermite approximant provides a polynomial M satisfying M (x, ŷ(x)) ≡ 0 mod x σ where σ = N (N + 1)/2 + N -1. Unfortunately σ < N 2 + 1 so that we have no way to ensure, using Lemma 18, that the Padé-Hermite approximant computed satisfies M (x, ŷ(x)) = 0. Consequently, using this method to compute the M i 's in Step (1b) of Algorithm ProbabilisticRationalFirstIntegral only provides a heuristic.

Proposition 29. Using Padé-Hermite approximation in Step (1b), Algorithm Prob-abilisticRationalFirstIntegral becomes a heuristic for computing a non-composite rational first integral of (S) of degree at most N using only Õ(N ω+2 ) arithmetic operations.

Proof. Beckermann-Labahn's algorithm (see [START_REF] Beckermann | A uniform approach for the fast computation of matrix-type Padé approximants[END_REF]) computes a Padé-Hermite approximant of type (N, N -1, . . . , 1) of the vector of power series (1, ŷ(x), . . . , ŷ(x) N ) T in Õ(N ω σ) arithmetic operations, where σ = N (N + 1)/2 + N -1. Using the proof of Theorem 27, we obtain the desired complexity estimate.

Implementation and experiments

The algorithms developed in the previous sections have been implemented in a Maple package called RationalFirstIntegrals5 .

Our implementation of the heuristic proposed in Subsection 5.5 is called Heuris-ticRationalFirstIntegral. It uses the gfun package [START_REF] Salvy | GFUN: a Maple package for the manipulation of generating and holonomic functions in one variable[END_REF] 6 and more precisely its seriestoalgeq command to search for a bivariate polynomial annihilating the power series computed using Padé-Hermite approximation.

We shall now illustrate our implementation and give some timings7 . 6.1. Comparison to previous methods. We start by comparing our implementation DeterministicRationalFirstIntegral to two previous methods, namely:

(1) the naive approach which consists in using the method of undetermined coefficients to search for two polynomials P and Q of degree at most N satisfying D(P ) Q-P D(Q) = 0. This implies solving a system of quadratic equations in the coefficients of P and Q. In our implementation, we use the solve command of Maple to solve the quadratic system, (2) the approach developed in [Chè11] based on the ecstatic curve.

Consider the planar polynomial vector field given by A(x, y) = -7 x + 22 y -55 and B(x, y) = -94 x + 87 y -56 which has no rational first integral of degree less than 6. The following table compares the timings (in seconds) of the different implementations for proving the non-existence of a rational first integral of degree less than N = 2, . . . , 6. If we now consider the vector field given by the polynomials A(x, y) = x + 2 and B(x, y) = -x 2 -2 x yy 2 -2 xy -2 which admits the rational first integral

x 2 +x y-2 x+y+1 of degree 2, we obtain the following timings (in seconds) depending on the degree bound N given in the input: In the latter table, the timings indicated for the "ecstatic curve method" correspond only to the computation of the N th ecstatic curve (which will be zero in all cases as there exists a rational first integral of degree 2) and not to the entire computation of a rational first integral of degree at most N which requires some more computations, see [Chè11, Subsection 5.2] for more details. The output of the two other methods consists in a rational first integral of degree 2.

P
The timings presented in this subsection illustrate that our implementation of DeterministicRationalFirstIntegral is significantly faster than our implementations of the two previous methods considered both in the case where there exists no rational first integral and in the case where there exists a rational first integral. This is coherent with the complexity analysis developed in this article. 6.2. Generic polynomial vector fields. If we choose at random two bivariate polynomials A and B, then the associated planar polynomial vector field has generically no rational first integral. In this subsection, we show that our implementation of DeterministicRationalFirstIntegral detects quickly the non-existence of rational first integral of these generic polynomial vector fields. The following table of timings is constructed as follows: for d = 1, . . . , 10, we generate two randomized bivariate polynomials A and B of degree d using the randpoly command of Maple and we check that A(0, y) ≡ 0. Then, we run DeterministicRationalFirstIntegral with N = d, . . . , 10 and we indicate the timings (in seconds) for detecting the nonexistence of a rational first integral of degree at most N , i.e., for returning the output "None". and Q(x, y) = 513 x 4 -342 x 3 y + 57 x 2 y 2 + 1710 x 3 -1254 x 2 y + 342 x y 2 -38 y 3 -10869 x 2 -1026 x y + 171 y 2 -37224 x + 12408 y -12560.

Two observations allow us to obtain a more compact form for F . First, looking at the syzygy in the leading term in x 4 , we see that 513 P (x, y) + 216 Q(x, y) = 2201580 (x 2 + 3 xy + 1).

Secondly, the discriminant of Pc Q shows that 117 P + 89 Q has a multiple factor, namely 117 P (x, y) + 89 Q(x, y) = 755 3 x 2 + 6 x -2 y + 1 (2 + 3 xy) 2 .

It follows that we have the following "nicer" rational first integral:

F (x, y) = x 2 + 3 x -y + 1 (3 x 2 + 6 x -2 y + 1) (2 + 3 x -y) 2 .
This simplification heuristics (using the spectrum) of the expression of a rational first integral to a more compact form can be obtained automatically by running the command SimplifyRFI of our package RationalFirstIntegrals.

In this example, the generic algorithm GenericRationalFirstIntegral run with N = 4 takes 0.342 seconds to compute a rational first integral; we see that, though it is 15 times slower than HeuristicRationalFirstIntegral (or DeterministicRationalFirstIntegral), it still has good performances on relatively small degrees.

Let us now have a look at the polynomial vector field given by A(x, y) = -18 x 8 y 8 -20 x 6 y 9 -6 x 2 y 12 + 24 x 10 y 3 -6 x 4 y 9 -4 y 13 -3 x 12 -7 x 2 y 10 , B(x, y) = 2 x -16 x 6 y 9 + 8 x 14 -18 x 4 y 10 -2 y 13 + 10 x 8 y 4 -2 x 2 y 10 -2 x 10 y -3 y 11 , considered by A. Ferragut in one of his talks concerning [START_REF] Ferragut | A new algorithm for finding rational first integrals of polynomial vector fields[END_REF]. It admits a rational first integral of degree 18. We have run our implementations of Heurist-icRationalFirstIntegral and ProbabilisticRationalFirstIntegral with the given bounds N = 3, 6, 9, 12, 15, and 18 in the input. The following table presents the outputs and the timings (in seconds) that we have obtained:

X X X X X X X X X X In the latter table, ? means that our implementation returns "I don't know" and F = P/Q is the rational first integral of degree 18 given by P (x, y) = -24 x 2 y 9 + 24 x 10 -24 y 10 , Q(x, y) = 8 x 18 -24 x 12 y 4 + 12 x 14 y + 24 x 6 y 8 -24 x 8 y 5 + 6 x 10 y 2 -8 y 12 + 44 x 2 y 9 Note that we obtain approximatively the same timings if we run the deterministic algorithm DeterministicRationalFirstIntegral instead of ProbabilisticRationalFirstIntegral. We can remark that our implementation of HeuristicRationalFirstIntegral is faster in this example than our implementation of ProbabilisticRationalFirstIntegral when there exists a rational first integral whereas ProbabilisticRationalFirstIntegral is much faster at discarding cases when no rational first integral exists. Moreover, we can see that, in this example, HeuristicRationalFirstIntegral only returns "I don't know" for N = 6, 9, 12, 15 whereas in these cases, ProbabilisticRationalFirstIntegral proves that there is no rational first integral of degree at most N . Note that these two drawbacks of HeuristicRationalFirstIntegral come from our implementation, which uses the command seriestoalgeq of the gfun package, and not from the algorithm itself.

Algorithm

In this example, if we replace Q by P + 3 4 Q, we obtain a new rational first integral F = P P + 3 4 Q which has a "nicer" (more compact) form F (x, y) =

x 2 y 9x 10 + y 10

(2 x 6 -2 y 4 + x 2 y) 3 .

This simplification of the expression of the rational first integral P/Q to a more compact form is obtained with the command SimplifyRFI of our package. 6.5. A hypergeometric example. Consider the family of polynomial vector fields given by A = 4 n 2 (x -1) (x + 1) and B = 1 + -4 n 2 x 2 + 4 n 2 y 2 -4 x y n 2 . For each integer n ∈ N * , it admits a rational first integral of degree N = 4 n + 1. This system is derived from the Riccati equation of a standard hypergeometric equation with a finite dihedral differential Galois group, see [vHW05]. The following table contains the timings (in seconds) for HeuristicRationalFirstIntegral to find a rational first integral of degree N = 4 n + 1 when it is run with N = 4 n + 1. In short, it takes 2 minutes to compute a rational first integral of degree 25 and 52 minutes to compute a rational first integral of degree 41 for this family of examples. 6.6. An Abel equation. We consider the rationally integrable Abel differential equation (3) in the article of Gine and Llibre [GL10]. It corresponds to the polynomial vector field given by A(x, y) = x (8 y -9) and B(x, y) = 3 y 2x -3 y. A rational first integral of degree 12 is computed in 4.142 seconds by HeuristicRa-tionalFirstIntegral and in 31.976 seconds by DeterministicRationalFirstIntegral if they are both run with N = 12. The rational first integral returned by HeuristicRational-FirstIntegral is given by P/Q with P (x, y) = 80 y 12 + 480 xy 10 + 1200 x 2 y 8 -1440 xy 9 + 1600 x 3 y 6 -5760 x 2 y 7 + 1200 x 4 y 4

-8640 x 3 y 5 + 8640 x 2 y 6 + 480 x 5 y 2 -5760 x 4 y 3 + 13248 x 3 y 4 + 80 x 6 -1440 x 5 y + 576 x 4 y 2 -13248 x 3 y 3 -4032 x 5 + 36288 x 4 y -27216 x 4 and Q(x, y) = 3 y 12 + 18 xy 10 + 45 x 2 y 8 -54 xy 9 + 60 x 3 y 6 -216 x 2 y 7 + 45 x 4 y 4 -324 x 3 y 5 + 324 x 2 y 6 + 18 x 5 y 2 -216 x 4 y 3 + 680 x 3 y 4 + 3 x 6 -54 x 5 y + 388 x 4 y 2 -680 x 3 y 3 + 32 x 5 -288 x 4 y + 216 x 4

Using the SimplifyRFI procedure, we find a rational first integral written in a more compact form:

F (x, y) = y 4 + 2 y 2 x + x 2 -6 yx 3 x 3 (4 y 4 + 8 y 2 x -4 y 3 + 4 x 2 -36 yx + 27 x)

Computation of Darboux polynomials

In this section, we show how the approach used above for computing rational first integrals of (S) of degree bounded by a fixed N ∈ N can be slightly modified for computing all irreducible Darboux polynomials for the derivation D associated (S) of degree at most N . In the output of our algorithms, irreducible Darboux polynomials in K[x, y] will be given by M

(c, x, y) ∈ K[c, x, y] and f (c) ∈ K[c]. The univariate polynomial f (c) is irreducible in K[c] and for all roots c i of f (c), we have an irreducible Darboux polynomial M (c i , x, y) ∈ K[x, y].
7.1. A deterministic algorithm. In this section we give a deterministic algorithm for computing all irreducible Darboux polynomials for the derivation D associated with (S) of degree at most N . This algorithm is divided into two steps. First, we compute all irreducible Darboux polynomials M (x, y) such that M (0, y) ∈ K: this is the task of Algorithm IrreducibleDarbouxPolynomialsPartial below applied to A and B. Then, in a second step, we show how we can compute the missing Darboux polynomials (those satisfying M (0, y) ∈ K) by applying IrreducibleDarboux-PolynomialsPartial to relevant polynomials constructed from A and B by a change of coordinates.

In these algorithms we suppose A(0, y) ≡ 0 and A(0, y), B(0, y) coprime. We can easily reduce our study to this situation. We have already explained how we can get A(0, y) ≡ 0. Now, we just have to remark that the second condition corresponds to the choice of an element which is not a root of the resultant Res y (A(x, y), B(x, y)). Thus after a finite number of shifts, we can assume that A(0, y) ≡ 0 and that A(0, y) and B(0, y) are coprime. In particular, this implies that x is not a Darboux polynomial and if M is a Darboux polynomial, then M (0, y) ≡ 0 in K[y]. We also assume that D would have no rational first integral with degree at most N . Indeed, from Theorem 6, in this situation D has an infinite number of irreducible Darboux polynomials. We can check this hypothesis with the previous algorithms.

Algorithm IrreducibleDarbouxPolynomialsPartial

Input: A, B ∈ K[x, y] s.t. A(0, y) ≡ 0, A(0, y), B(0, y) coprime, and a bound N ∈ N such that (S) has no rational first integral of degree at most N . Output: The set of all irreducible Darboux polynomials M for the derivation D order to get a new polynomial M such that M (0, y) has a root. If M (0, y) ∈ K, then M has a root at infinity. Thus we consider the following change of coordinates: we set

A ♯ (x, y, z) = A x z , y z z d , B ♯ (x, y, z) = B x z , y z z d , M ♯ (x, y, z) = M x z , y z z k ,
where k = deg(M ), and we consider the following polynomials:

Ã(y, z) = A ♯ (1, y, z), B(y, z) = B ♯ (1, y, z), M (y, z) = M ♯ (1, y, z).
A straightforward computation shows that:

Lemma 31. With the above notation, if M is a Darboux polynomial for the derivation D = A(x, y) ∂ ∂x +B(x, y) ∂ ∂y , then M is a Darboux polynomial for the derivation

D = -y Ã(y, z) + B(y, z) ∂ ∂y -Ã(y, z) z ∂ ∂z .
Furthermore, if M (0, y) ∈ K \ {0}, then M (0, z) ∈ K.

We deduce the following algorithm:

Algorithm IrreducibleDarbouxPolynomials

Input: A, B ∈ K[x, y] s.t. A(0, y) ≡ 0, A(0, y), B(0, y) coprime, B(0, z) ≡ 0, Ã(0, y) and B(0, y) coprime, and a bound N ∈ N such that (S) has no rational first integral of degree at most N .

Output:

The set of all irreducible Darboux polynomials M for the derivation D such that deg(M ) ≤ N .

(1) E := IrreducibleDarbouxPolynomialsPartial(A, B, N ).

(2) E ′ := IrreducibleDarbouxPolynomialsPartial(-y à + B, -à z, N ).

( For the same reasons as before, using a finite number of shifts we can suppose that the hypotheses "A(0, y) ≡ 0, A(0, y), B(0, y) coprime, B(0, z) ≡ 0, Ã(0, y), B(0, y) coprime" are satisfied so that these conditions are not restrictive.

As a direct consequence of Proposition 30 and Lemma 31, we obtain the following result.

Proposition 32. Algorithm IrreducibleDarbouxPolynomials is correct.

A probabilistic algorithm.

As we have seen in Section 5, the computation of a basis of solutions of a system of linear equations is the most costly step of our algorithms. In IrreducibleDarbouxPolynomials, we have to consider four systems of linear equations. The first reason is that in IrreducibleDarbouxPolynomialsPartial, we need to study two linear systems in order to take into account the situation where x = 0 is a vertical tangent of the curve M (x, y) = 0. Indeed, in this situation we can not get a parametrization x, y(x) of the curve. The second reason is that we need to use a change of coordinates in order to control the situation where M (0, y) has a root at infinity. Of course, for a generic polynomial vector field, these two situations (i.e., a vertical tangent and a root at infinity) do not appear. We then deduce the following probabilistic algorithm.

Algorithm ProbabilisticIrreducibleDarbouxPolynomials

Input: A, B ∈ K[x, y], a bound N ∈ N such that (S) has no rational first integral of degree at most N , and two elements x 0 , α ∈ K. Output: The set of all irreducible Darboux polynomials M for the derivation D such that of deg(M ) ≤ N .

(1) E := ∅.

(2) Set A α (x, y) = A(x + α y, y)α B(x + α y, y), B α (x, y) = B(x + α y, y) and

D α = A α (x, y) ∂ ∂x + B α (x, y) ∂ ∂y .
(3) For an indeterminate c, compute the polynomial y c ∈ K(c)[x] of degree ≤ (N 2 + 1) s.t. y c (x 0 ) = c and dyc dx ≡ Bα(x,yc) Aα(x,yc) mod x N 2 +1 . (4) Let M (x, y) = N i=0 N -i j=0 m i,j x j y i be an ansatz for the Darboux polynomials that we are searching for.

(5) Construct the linear system L(c) for the m i,j 's given by: M (x, y c (c, x)) ≡ 0 mod x N 2 +1 . Proof. First, we remark that M is a Darboux polynomial for D if and only if M α (x, y) = M (x + αy, y) is a Darboux polynomial for D α . Thus the strategy used in this algorithm is to perform a change of coordinates in order to be in a generic position, and then to compute all irreducible Darboux polynomials M α of degree at most N by considering only one linear system.

The proof of Proposition 30 shows that from Step (2) to Step (8), we compute all irreducible Darboux polynomials satisfying:

M α (x 0 , y) ∈ K and Res y (M α (x 0 , y), A α (x 0 , y)) = 0.

Let us study the probability to get M α (x 0 , y) ∈ K. If M (x, y) = 0≤i+j≤N a i,j x i y j , then M α (x 0 , y) = ( i+j=N a i,j α i ) y N + • • • where the other terms have degree relatively to y strictly less than N . Thus, if i+j=N a i,j α i is not equal to zero, then we have M α (x 0 , y) ∈ K. As (S) has no rational first integral of degree at most N , then by Darboux-Jouanolou's theorem (see Theorem 6), we have at most B(d)+1 irreducible Darboux polynomials with degree at most N . Thus, by Zippel-Schwartz's lemma, the probability to reach the situation M α (x 0 , y) ∈ K for all irreducible Darboux polynomials is at least 1 -(B(d) + 1) N/|Ω|. Now we suppose that M α (x 0 , y) ∈ K and we study the probability to have the situation Res y (M (x 0 , y), A(x 0 , y)) = 0. If the polynomial Res y (M (x, y), A(x, y)) is not zero, then, by Zippel-Schwartz's lemma, the probability to reach this situation for all irreducible Darboux polynomials, is at least 1 -(B(d) + 1) N d/|Ω|. If the polynomial Res y (M α (x, y), A α (x, y)) is zero, then M α and A α have a common factor. As we suppose M α irreducible, we deduce that M α divides A α . Thus M α divides B α ∂ y (M α ). As A α and B α are coprime, we get that M α divides ∂ y (M α ). This situation is possible only when deg y (M α ) = 0. This means Res y (M α (x, y), A α (x, y)) ≡ 0 when deg y (M α ) = 0 and M α divides A α (x, 0). We compute this kind of irreducible Darboux polynomials in Step (10) of the algorithm.

In conclusion, the algorithm computes all irreducible Darboux polynomials of degree at most N with the announced probability estimate.

7.3. Implementation and example. We have implemented the algorithm Prob-abilisticIrreducibleDarbouxPolynomials in our package RationalFirstIntegrals8 . Let us illustrate the purpose of this section on an interesting example.

Consider the vector field corresponding to the jacobian derivation associated with f (x, y) = (yx -1) (xy 2 ) (x y -1), namely, A(x, y) := -∂f ∂y (x, y) = -3 x 2 y 2 +4 xy 3 +x 3 -2 x 2 y-3 xy 2 +x 2 +2 xy-3 y 2 +x+2 y, B(x, y) := ∂f ∂x (x, y) = 2 xy 3y 4 -3 x 2 y + 2 xy 2 + y 3 -2 xyy 2 + 2 xy + 1.

By construction, it admits the rational first integral f of degree 4 and the Darboux polynomials M 1 (x, y) = yx -1, M 2 (x, y) = xy 2 , M 3 (x, y) = x y -1 of degree at most 2. Let us consider the computation of all irreducible Darboux polynomials of degree at most N = 2. The first Darboux polynomial M 1 satisfies M 1 (0, y) = y -1 ∈ K and its root c M1 = 1 satisfies A(0, c M1 ) = -1 = 0. Therefore it will be found by considering the linear system L 1 (c) in IrreducibleDarbouxPolynomialsPartial, see the proof of Proposition 30. The Darboux polynomial M 2 satisfies M 2 (0, y) = -y 2 ∈ K and but its root c M2 = 0 satisfies A(0, c M2 ) = 0. Thus it will be missed if we only consider system L 1 (c) in IrreducibleDarbouxPolynomialsPartial. It is the case where the curve M 2 (x, y) = 0 has the vertical tangent x = 0. However, if we consider the second system L 2 (c) in IrreducibleDarbouxPolynomialsPartial, we will find this Darboux polynomial, see the proof of Proposition 30. Finally M 3 satisfies M 3 (0, y) = -1 ∈ K so that M 3 (0, y) has a root at infinity. Considering only the systems L 1 (c) and L 2 (c) in IrreducibleDarbouxPolynomialsPartial will not be enough to find this Darboux polynomial. However, performing the change of coordinates as in IrreducibleDarbouxPolynomials and applying Irre-ducibleDarbouxPolynomialsPartial to -y à + B and -à z instead of A and B will provide this Darboux polynomial.

To summarize, applying IrreducibleDarbouxPolynomialsPartial to A and B, we get M 1 and M 2 but we miss M 3 but either applying IrreducibleDarbouxPolynomials or ProbabilisticIrreducibleDarbouxPolynomials we get the three Darboux polynomials. Note also that applying an algorithm similar to ProbabilisticIrreducibleDarbouxPolynomials but where we skip Step (2), i.e., we do not perform the generic change of coordinate, we would obtain only M 1 and miss both M 2 and M 3 . Our implementation of ProbabilisticIrreducibleDarbouxPolynomials requires computations in K(c) so that as for GenericRationalFirstIntegral it is not very efficient and can not be used in practice for examples with large degrees. To give an idea of timings 9 , on the previous example, running ProbabilisticIrreducibleDarbouxPolynomials without the change of coordinates in Step (2), we obtain {M 1 } in 2.737 seconds but we miss M 2 and M 3 whereas running ProbabilisticIrreducibleDarbouxPolynomials, we get the complete set {M 1 , M 2 , M 3 } in 9.775 seconds.

( 2 )

 2 While |Ω| ≤ 2 N (B(d) + 1) + 2 do (a) Choose two random elements c 1 , c 2 ∈ K\Ω s.t. c 1 = c 2 and A(0, c i ) = 0 for i = 1, 2. (b) F := ProbabilisticRationalFirstIntegral(A, B, (c 1 , c 2 ), N ). (c) If F ="None", then Return "None". (d) Else if F =["I don't know",[e 2 ]], then Ω := Ω ∪ {c 1 , e 2 } and go to

  Heuristic 0.540 12.548 118.804 592.494 3247.325

)

  For all [ M (c, y, z), P(c)] ∈ E ′ do: (a) M (c, x, y) := M (c, y x , 1 x ) x deg(M) . (b) Add [M (c, x, y), P(c)] to E. (4) Return E.

( 6 )

 6 Clear the denominator in L(c). (7) Compute the Smith normal form of L(c). Let P(c) be the last invariant factor of L(c). (8) Factorize P(c) over K: P(c) = s i=1 P i (c). (a) For i from 1 to s do: (i) Set K[c i ] = K[c]/(P i (c)).(ii) Compute a solution of L(c i ) s.t. the corresponding polynomial M i has minimal degree in y and is primitive w.r.t. y.(iii) If gcd(D α (M i ), M i ) = M i , then E := E∪{[M i (c, x-α y, y), P i (c)]}. (9) Factorize A α (x, 0) over K: A α (x, 0) = k i=1 A i (x). (10) For i from 1 to k do: (a) If gcd(D α (A i ), A i ) = A i , then E := E ∪ {[A i (xα y, y), c -1]}. (11) Return E.Proposition 33. The algorithm ProbabilisticIrreducibleDarbouxPolynomials is correct. Furthermore, if x 0 and α are chosen uniformly at random in a finite set Ω ⊂ K such that |Ω| > N d (B(d) + 1), then the probability that this algorithm returns all irreducible Darboux polynomials is at least 1 -N (B(d) + 1) |Ω| 1 -N d (B(d) + 1) |Ω| .

Note that Remark 1 in[START_REF] Ferragut | A new algorithm for finding rational first integrals of polynomial vector fields[END_REF] (and the example that follows it) suggests a strategy to obtain linear equations, via Puiseux series at a resonant node.

It is available with various examples at http://www.unilim.fr/pages_perso/thomas.cluzeau/Packages/RFI/RationalFirstInte

i.e., a basis over K(c), see Remark 19.

It is available with some examples at http://www.unilim.fr/pages_perso/thomas.cluzeau/Packages/RFI/RationalFirstIntegr

http://perso.ens-lyon.fr/bruno.salvy/?page_id=48.

All the computations were made on a 2.7 GHz Intel Core i7.

It is available at http://www.unilim.fr/pages_perso/thomas.cluzeau/Packages/RFI/RationalFirstIntegrals.html.

which admits the rational first integral of degree 4

If we run ProbabilisticRationalFirstIntegral with the bound N = 4 and c 1 = 0 or c 2 = 0 in the input, then we get "I don't know". The reason why our algorithm fails is that (Q(0, 0) : P (0, 0)) = (-1 : 0) ∈ σ(P, Q) since -P (x, y) = -(yx) x 2 + y -1 is a reducible polynomial (and also a polynomial of degree less than N = 4). Of course, running ProbabilisticRationalFirstIntegral with values of c 1 and c 2 such that (Q(0, c i ) : P (0, c i )) ∈ σ(P, Q) for i = 1, 2 provides the correct output, i.e., a rational first integral of degree N = 4; see the explanations at the end of Subsection 4.1. The deterministic algorithm DeterministicRationalFirstIntegral calls recursively Prob-abilisticRationalFirstIntegral and exploits the fact that there only exists a finite number of such bad values of the c i 's. So in this example, it returns correctly a rational first integral of degree N = 4. A first integral of degree 4 was found in 12 seconds using their algorithm (see [START_REF] Ferragut | A new algorithm for finding rational first integrals of polynomial vector fields[END_REF]) which was a notable improvement on previous methods. In comparison, running HeuristicRationalFirstIntegral (or DeterministicRationalFirstIntegral) with N = 4 we get such a rational first integral F = P/Q in 0.022 seconds, where P (x, y) = -216 x 4 + 144 x 3 y -24 x 2 y 2 -720 x 3 + 528 x 2 y -144 x y 2 +16 y 3 + 8868 x 2 + 432 x y -72 y 2 + 28548 x -9516 y + 9580, such that deg(M ) ≤ N and M (0, y) ∈ K.

(1) E := ∅.

(2) For an indeterminate c, compute the polynomial y c ∈ K(c)[x] of degree at most (N 2 + 1) s.t. y c (0) = c and dyc dx ≡ B(x,yc) A(x,yc) mod x N 2 +1 .

(3) For an indeterminate c, compute the polynomial

N -i j=0 m i,j x j y i be an ansatz for the Darboux polynomials that we are searching for.

(5) Construct the linear system L 1 (c) for the m i,j 's given by:

(6) Construct the linear system L 2 (c) for the m i,j 's given by: M (x c (c, y), y) ≡ 0 mod y N 2 +1 . Proof. Let M be an irreducible Darboux polynomial such that M (0, y) ∈ K and c M be a root of M (0, y). Then we have: A(0, c M ) = 0 or B(0, c M ) = 0 because A(0, y) and B(0, y) are assumed to be coprime. If A(0, c M ) = 0 and M (0, c M ) = 0, then M admits a root y cM ∈ K(x) such that y cM (0) = c M . Then, from Proposition 17, y cM is a power series solution of (E). Thus c M is a root of P 1 (c). Then, by Lemma 18, M is constructed in Step (7(d)ii). If for a constant c M , we have B(0, c M ) = 0 and M (0, c M ) = 0, then the previous arguments used with P 2 (c) show that M is also constructed.

In the algorithm IrreducibleDarbouxPolynomialsPartial, we compute irreducible Darboux polynomials M such that M (0, y) ∈ K. Indeed, the algorithm finds a irreducible Darboux polynomial M if and only if the curve M (x, y) = 0 and the line x = 0 have an intersection point. Now, we show how to get irreducible Darboux polynomials such that M (0, y) ∈ K. The idea is to use a change of coordinates in


Laurene Aoun

Stanislas Larnier

Pierre Weiss

Ariane Herbulot

Bernard Ducommun

Valérie Lobjois

Christophe Vieu

PDMS MICRODEVICE ARRAYS FOR MEASURING THE FORCES EXERTED BY GROWING MULTICELLULAR TUMOR SPHEROIDS

Keywords: Micro-tumors, Spheroids, Force, Volumetric imaging, Mechanical simulation

This paper presents a novel technique for measuring the mechanical forces exerted by a living MultiCellular Tumor Spehroid (MCTS) during growth. We used high aspect ratio microfabricated polydimethylsiloxane (PDMS) pillars as force sensors. During growth, spheroids induce a deformation of the micropillars, which is the parameter we used in order to measure the exerted forces. From 3 dimensional imaging of micropillars induced deflection and mechanical simulation, we were able to measure the force exerted by mammary cancer cell (MCF7) spheroids. Using this methodology we found a value of force of the order of 100-300 nanoNewtons.

INTRODUCTION

Tumor cells are known to be sensitive to external mechanical stimuli but little is known about the intrinsic growth-associated mechanical properties of micro-tumors [START_REF] Hirschhaeuser | Multicellular tumor spheroids: an underestimated tool is catching up again[END_REF]. In that context, the changing forces that a micro-tumor exerts on its microenvironment need to be better characterized. A MCTS is a 3D complex system that reproduces cell-cell and cell-matrix interactions and exhibits a cellular architecture that mimics micro-tumor organization. Many approaches have been developed to study the mechanical properties of spheroids. It includes cellular capsules to measure the pressure exerted by a spheroid [START_REF] Alessandri | Cellular capsules as a tool for multicellular spheroid production and for investigating the mechanics of tumor progression in vitro[END_REF], parallel plate tensiometry for rheological studies [START_REF] Vasilica Stirbat | Multicellular aggregates: a model system for tissue rheology[END_REF] or the micropipette technique to study the response of cellular aggregates to external stresses [START_REF] Guevorkian | Aspiration of biological viscoelastic drops[END_REF]. Our work differs from those mentioned before, firstly, by its high throughput and the ability to investigate a large number of spheroids within the array of microdevices at the same time. Secondly, our system enables to measure the force directly exerted by the spheroid during growth on its surrounding in the nanoNewton scale and for a large period of time without a drastic perturbation of the cell organization within the spheroid under investigation.

EXPERIMENTAL

The produced devices are arrays of high aspect ratio (1:10), round and flexible PDMS pillars, of 300 µm in height, identical in diameter and circularly distributed around a space where the spheroid will be placed. Each pillar is used as a force sensor by a simple recording of its bending. We fabricated devices of different pillar diameters in order to study the effect of the stiffness on the intensity of the generated forces. The technological challenge behind these devices was the optimization of a process capable to produce PDMS beams high enough compared to the MCTs (around 300 µm), soft enough (stiffness around 10 nN/µm) for measuring nano-newton scale forces, while preventing their mechanical collapse [START_REF] Aoun | Microdevice arrays of high aspect ratio poly(dimethylsiloxane) pillars for the investigation of multicellular tumour spheroid mechanical properties[END_REF]. In order to measure forces exerted by growing spheroids, the deformations of the pillars and the contact zone of the spheroid with the microdevice are required. The growth of cancer cells spheroids was recorded by taking 3D confocal images for fluorescently stained microdevices and cells expressing a fluorescent marker (Figure 1).

RESULTS AND DISCUSSION

Spheroid growth induced bending of the micropillars, reflecting the forces exerted by the spheroid. On one hand, in order to characterize the 3 dimensional deformations of the pillars from the 3D confocal images, we developed a homemade MATLAB program that enables the quantification of a large number of devices at one time point (Figure 2).

Figure 2: Image analysis program for a microdevice with slightly displaced pillars induced by MCF7 spheroid (unlabelled), (a) a 3D confocal view of the whole microdevice to analyze the deformation of the pillars, (b) a representation of the detected centers of all pillars for each slice (80 slices with a z-step of 3,75µm) and (c) a 3 D reconstruction of the smoothed pillars with the microdevice's parameters.

On the other hand, we used finite element method through the software COMSOL, to simulate the bending of the PDMS pillars. From various images, we were able to localize the contact region where the spheroid exerts its efforts on the pillars. Then by simulation, we could calculate the displacement field and the stress accumulated in the pillars for a large range of applied forces. The simulation results were then compared to the 3D experimental displacement in order to extract the value of the force responsible for this deformation (Figure 3).

Figure 3: (a) A simulation example of pillar deformation and Von Mises stress after the application of a distributed force in the region of contact with the spheroid, (b) bending profiles of a PDMS pillar showing the closest match between the simulated result (blue) and the experimental one (red).

Using this procedure of force calculation, we recorded the forces exerted by MCF7 spheroids on the microdevices after a period of 4 days of contact between the spheroids and the microdevices. Our results showed that the forces exerted on a device with micropillars of 28 µm in diameter are in the order of 100-300 nN. Moreover we have found that spheroids confronted with stiffer pillars (36 µm in diameter), exerted slightly higher forces as shown in Figure 4.

CONCLUSION

In this study, we offered a new methodology to investigate systematically the mechanical properties of spheroids with a high throughput. We gave for the first time a range of forces exerted by a living MCTS on its surroundings. In previous work we have shown the dynamic interaction between the spheroids and the microdevice [START_REF] Aoun | Microdevice arrays of high aspect ratio poly(dimethylsiloxane) pillars for the investigation of multicellular tumour spheroid mechanical properties[END_REF], which opens a new perspective towards kinetics studies of the evolution of forces within a growing spheroid, in time. This method also opens original perspectives for pharmacological drug screening and for investigating the mechanisms of growth of micro-tumors.

AKNOWLEDGEMENTS

This work is funded by the foundation 2RITC (Toulouse).

Figure 1 :

 1 Figure 1: (left) Scanning Electron microscopy images of an array of high aspect ratio PDMS microdevices and (right) a 3 dimensional confocal image of a fluorescently stained microdevice (green) surrounding a spheroid expressing a fluorescent marker (red).

Figure 4 :

 4 Figure 4: A scattered plot with the standard deviation of the forces exerted by MCF7 spheroids after 4 days of growth within the microdevices for 2 pillars diameters: 28 and 36µm.