Égoire Gr

Lecerf

NEW RECOMBINATION ALGORITHMS FOR BIVARIATE POLYNOMIAL FACTORIZATION BASED ON HENSEL LIFTING

Keywords: Mathematics Subject Classification. Primary 12Y05, 68W30; Secondary 11Y16, 12D05, 13P05 Polynomial factorization, Hensel lifting

We present new faster deterministic and probabilistic recombination algorithms to compute the irreducible decomposition of a bivariate polynomial via the classical Hensel lifting technique. For the dense bi-degree polynomial representation, the costs of our recombination algorithms are essentially sub-quadratic.

Introduction

Throughout this paper K denotes a commutative field of characteristic p, and F is a polynomial in K[x, y] of bi-degree (d x , d y ), which means degree d x in x and d y in y. We are interested in the complexity of the computation of the irreducible factors F 1 , . . . , F r of F together with their respective multiplicities e 1 , . . . , e r for the dense polynomial representation, which means that a polynomial of bi-degree (d x , d y ) is stored as the vector of its coefficients in the basis of the monomials of bi-degree at most (d x , d y ). In this model, the size of F is precisely (d x + 1)(d y + 1).

Main Results. For the cost analysis we use the computation tree model and the classical O and Õ (read "soft Oh") notation in the neighborhood of infinity as defined in [START_REF] Zur Gathen | Modern computer algebra[END_REF]Chapter 25.7]: f (d) ∈ Õ(g(d)) means that f (d) ∈ g(d)(log 2 (3 + g(d))) O(1) . The factorization algorithms presented in this paper work for any characteristic but the coefficient field K must satisfy a few requirements: its cardinality must be sufficiently large, irreducible factorization in K[y] must be computable; when in positive characteristic, pth root extraction in K must be computable, and the F p -algebra structure of K must be effective. Herein F p represents the prime subfield of K. These requirements are detailed in Section 2.1.

For the sake of simplicity we will provide the reader only with a complete cost analysis in two cases: (a) in characteristic 0 in terms of the number of arithmetic operations in K; (b) when K is the finite field F p k in terms of the number of arithmetic operations in F p .

The constant ω used below represents a feasible matrix multiplication exponent. We require that 2 < ω ≤ 3 (see Section 2.1). The first theorem summarizes the cost of our deterministic reduction from two to one variables: Theorem 1. Assume that K has cardinality at least 2d x d y + max(d x , d y ) + 1. Then the computation of the irreducible decomposition (F 1 , e 1 ), . . . , (F r , e r ) of F reduces to the computation of irreducible decompositions of polynomials in K[y] whose degree sum is at most d x + d y , plus a. in characteristic 0, Õ((d x d y ) (ω+1)/2 ) arithmetic operations in K; b. if K := F p k , Õ(k(d x d y ) (ω+1)/2 )) arithmetic operations in F p .

Our second theorem yields a faster reduction but with a probability of failure. Therein the function R, defined in Section 2.1, is related to the cost of random subsets generation.

Theorem 2. Assume that K has cardinality at least 10d x d y . Then the computation of the irreducible decomposition (F 1 , e 1 ), . . . , (F r , e r ) of F can be done with a probabilistic algorithm which performs irreducible factorizations of polynomials in K[y] whose degree sum is at most d x + d y , plus a. in characteristic 0, O((d x d y ) 1.5 ) operations in K and R(O(d x d y )); b. if K := F p k , Õ(k(d x d y ) 1.5 ) operations in F p and R(O(d x d y )). The algorithm outputs either nothing or a correct result with a probability at least 1/2. The proofs of these two theorems are given at the end of Section 3.2. Informally speaking, if we consider that an arithmetic operation in K is not cheaper than k operations in F p , then the costs in parts (a) and (b) of the above theorems are of the same order of magnitude.

The mathematical ingredients of our algorithms are presented in the first section. Then the second section describes the core algorithms. The third section is devoted to the whole factorization algorithm.

Overview of the Algorithms. The algorithms underlying Theorems 1 and 2 both decompose into the two following stages:

1. Reduction to the separable case. Up to exchanging x and y we can assume that d y ≤ d x . We compute the separable decomposition of F in y by means of the algorithm presented in [START_REF] Lecerf | Fast separable factorization and applications[END_REF]. 2. Factorization of the separable components. In this stage F is assumed to be separable in y. Assuming that the cardinality of K is sufficiently large, a suitable shift of the variable x reduces the problem to the normalized case defined as follows:

(H)

       deg(F (0, y)) = d y ,
Res F (0, y), ∂F ∂y (0, y) = 0, F is primitive when seen in K[x][y],

where Res(A, B) denotes the resultant of the two univariate polynomials A and B. In this case we can use the classical Hensel lifting strategy: from the irreducible factors F 1 (0, y), . . . , F s (0, y) of F (0, y), we compute the irreducible analytic decomposition F 1 , . . . ,

F s of F in K[[x]][
y] to a certain finite precision. Then we discover how these analytic factors recombine into F 1 , . . . , F r . The algorithms underlying Theorems 1 and 2 require a precision of the series only (x dx+1 ). In characteristic 0 or sufficiently large, these results are to be compared to ones of [START_REF] Lecerf | Sharp precision in Hensel lifting for bivariate polynomial factorization[END_REF] where the required precision is (x 2d ) with d being the total degree of F . In positive characteristic, we improve the best previously known recombination algorithm designed in [START_REF] Belabas | Factoring polynomials over global fields[END_REF] that requires a precision (x dx(2dy-1)+1 ) in the worst case [BHKS04, Theorem 5.2].

Related Works. The first polynomial time multivariate factorization algorithms are due to Kaltofen in the beginning of the eighties [Kal82b, Kal82c, Kal85]: they were derived from the algebraic approximant technique invented by A. K. Lenstra, H. W. Lenstra, and Lovász in 1982 to factor univariate polynomials over Q in polynomial time. Then Chistov, von zur Gathen, Grigoriev, Kaltofen, and A. K. Lenstra substantially contributed to this subject in order to reach polynomial time for all the usual fields. In 2003, Gao accomplished an important breakthrough with the design of a softly quadratic time probabilistic reduction of the factorization problem from two to one variable whenever the characteristic of the coefficient field is zero or sufficiently large [START_REF] Gao | Factoring multivariate polynomials via partial differential equations[END_REF]. For more historical details we refer the reader to [Kal82a, Gat84, Kal90, Kal92, Zip93, Kal95, Kal03, GG03, Gao03, CG05].

The Hensel lifting strategy was popularized in computer algebra by Zassenhaus in 1969 in the context of factoring in Q[y] via a p-adic completion of Q. For a long time, in the bivariate case, the cost of this approach has been remaining exponential in the number s of the analytic factors to recombine. Yet, as proved by Gao and Lauder, the cost of the recombination process is indeed softly linear in average [START_REF] Gao | Hensel lifting and bivariate polynomial factorisation over finite fields[END_REF] over finite fields, which explains the practical efficiency of this approach.

Of course a polynomial time factorization algorithm solves the recombination problem in polynomial time, but, to the best of our knowledge, the first polynomial time natural solution to this problem is due to van Hoeij for Q[y] [Hoe02], and to Belabas, van Hoeij, Klüners and Steel in the more general setting of global fields, that includes K(x)[y] [BHKS04]: the key idea is to recombine the logarithmic derivatives of the lifted factors. We refer to [START_REF] Lecerf | Sharp precision in Hensel lifting for bivariate polynomial factorization[END_REF] for a detailed historical survey on the recombination algorithms.

In [BLS + 04, Lec06, Lec07b] we presented recombination algorithms with a subquadratic cost in terms of the total degree, for the dense representation, and when the characteristic of K is 0 or sufficiently large. In the present paper we improve and generalize these results: our new algorithms require a smaller lifting precision and work whatever the characteristic of K is. In addition, instead of the total degree we deal with the bi-degree.

Reduction of the Recombination Problem to Linear Algebra

Throughout this section we assume that Hypothesis (H) holds. We explain how the recombination problem efficiently reduces to linear system solving. We separate the case when p is 0 or large enough to the specific case when p > 0. Yet we will see that the latter case makes use of the former case. The separation between these two cases is precisely given by the following condition:

(C) K has characteristic 0 or at least d x (2d y -1) + 1. 1.1. Main Notation. Let c ∈ K[x] (respectively, c i ∈ K[x] for each i ∈ {1, . . . , r}) denote the leading coefficient of F (respectively, of F i ) when seen in K[x][y].
For convenience we assume that c, c 1 , . . . , c r are all monic (i.e. their leading coefficients equal 1), so that we can write

F = F 1 • • • F r . Recall that F 1 , . . . , F s represent the monic irreducible factors of F seen in K[[x]][y].
Then, to each i ∈ {1, . . . , r}, we associate the vector µ i ∈ {0, 1} s , uniquely defined by

F i = c i s j=1 F µi,j j
.

(1)

Since the µ i have entries in {0, 1} and have pairwise distinct supports, up to a unique permutation, they form a reduced echelon basis, which means that the r × s matrix whose ith row is µ i is in the row-reduced echelon form. From now on we assume that µ 1 , . . . , µ r are actually ordered so that they form a reduced echelon basis. In this setting the recombination problem reformulates as follows: from the knowledge of the F i to a sufficient finite precision, compute all the vectors µ i . Then it is classical that the knowledge of the µ i suffices to recover F 1 , . . . , F r efficiently (see Algorithm 3 in Section 2.4).

If A := i,j≥0 a i,j x i y j denotes a polynomial in K[[x]][y] then we write ⌈A⌉ l and

[A] l k respectively for the following projections to K[x, y]:

⌈A⌉ l := 0≤i≤l-1, j≥0 a i,j x i y j , [A] l k := k≤i≤l-1, j≥0 a i,j x i y j .
We also introduce the following partial products:

Fi := r j=1,j =i F j = F F i and Fi := c s j=1,j =i F j = F F i .
Finally, the central objects of our recombination process are the following:

G i := Fi ∂F i ∂y dx+1
, for all i ∈ {1, . . . , s}.

The recombination algorithms will only make use of these objects, so that the only lifting to precision (x dx+1 ) is necessary. Throughout this paper F denotes a subfield of K, and the space of the bivariate polynomials over K of bi-degree at most (k, l) is written K[x, y] k,l . Our goal will be the computation of the following F-vector space:

L F := (ℓ 1 , . . . , ℓ s ) ∈ F s | s i=1 ℓ i G i ∈
Fi ∂F i ∂y = Fi ∂F i ∂y dx+1 =     s j=1 µ i,j Fj ∂F j ∂y     dx+1 = s j=1 µ i,j G j . (2)
Since

µ i ∈ {0, 1} s ⊆ F s , it follows that µ i ∈ L F .
By Hypothesis (H), G 1 (0, y), . . . , G s (0, y) are linearly independent over F, so are G 1 , . . . , G s . By the same argument, F1 ∂F1 ∂y , . . . , Fr ∂Fr ∂y are also linearly independent over F. Therefore L F has dimension r, which concludes the proof. 1.2. Partial Fraction Decomposition. Let (ℓ 1 , . . . , ℓ s ) be a vector in F s , let G := s i=1 ℓ i G i , and let K denote the algebraic closure of K. In order to construct a set of linear equations over F whose solution set is L F , we consider the partial fraction decomposition of G/F in K [[x]]. Let φ 1 , . . . , φ dy represent the roots of F in K[[x]], and let ρ i := G(x, φ i )/ ∂F ∂y (x, φ i ), for all i ∈ {1, . . . , d y }, so that we have:

G F = dy i=1 ρ i y -φ i . (3)
This partial fraction decomposition is well defined over K [[x]] thanks to Hypothesis (H).

Lemma 2. If (ℓ 1 , . . . , ℓ s ) ∈ L F then ρ 1 , . . . , ρ dy all belong to F. Conversely, if ρ 1 , . . . , ρ dy all belong to K then (ℓ 1 , . . . , ℓ s ) ∈ L F .

Proof. By definition, if (ℓ 1 , . . . , ℓ s ) ∈ L F then G is a F-linear combination of the F1 ∂F1 ∂y , . . . , Fr ∂Fr ∂y . Since, for all i, the residues of Fi ∂Fi ∂y /F = ∂Fi ∂y /F i are all in {0, 1}, we obtain that all the ρ i belong to F. Conversely, assume that all the ρ i are in K. By substituting 0 for x in (3) and in the definition of G, we obtain

G(0, y) F (0, y) = dy i=1 ρ i (0) y -φ i (0) = s j=1
ℓ j ∂Fj ∂y (0, y) F j (0, y) , so that Hypothesis (H) implies that ρ i = ρ i (0) = ℓ j whenever F j (0, φ i (0)) = 0, whence

G = s j=1 ℓ j Fi ∂F i ∂y .
Now let i and j in {1, . . . , s} be such that F i and F j both divide the same irreducible factor F k of F for some k ∈ {1, . . . , r}. By construction, F i and

F j respectively divide ℓ i ∂F ∂y -G and ℓ j ∂F ∂y -G in K[[x]][y]. Since the two latter expressions are in K[x, y], they are multiple of F k . It follows that F i divides ℓ j ∂F ∂y -G in K[[x]][y],
hence that ℓ i = ℓ j since ∂F ∂y is invertible modulo F i by Hypothesis (H). Therefore (ℓ 1 , . . . , ℓ s ) is a F-linear combination of the µ i and the conclusion thus follows from Lemma 1. 1.3. Characteristic Zero or Large Enough. We carry on with the notation: (ℓ 1 , . . . , ℓ s ) ∈ F s , G := s i=1 ℓ i G i , and ρ 1 , . . . , ρ dy are the residues of G/F defined in (3). In characteristic 0, in order to ensure that all the ρ i are in K, it suffices to ensure that all the ρ ′ i are zero. Letting D

: K[x, y] dx,dy-1 → K[x, y] 3dx-1,3dy-3 G → ∂G ∂x ∂F ∂y - ∂G ∂y ∂F ∂x ∂F ∂y - ∂ 2 F ∂xy ∂F ∂y - ∂ 2 F ∂y 2 ∂F ∂x G,
it is classical that ρ ′ i can be calculated as follows:

ρ ′ i = d dx G(x, φ i (x)) ∂F ∂y (x, φ i (x)) = ∂G ∂x (x, φ i (x)) + ∂G ∂y (x, φ i (x))φ ′ i (x) ∂F ∂y (x, φ i (x)) - ∂ 2 F ∂xy (x, φ i (x)) + ∂ 2 F ∂y 2 (x, φ i (x))φ ′ i (x) ∂F ∂y (x, φ i (x)) 2 G(x, φ i (x)) = D(G)(x, φ i (x)) ∂F ∂y (x, φ i (x)) 3 , since φ ′ i (x) = -∂F ∂x (x, φ i (x))/ ∂F ∂y (x, φ i (x)).
Testing the vanishing of all the ρ ′ i thus becomes equivalent to testing if

F divides D(G) in K[[x]][y]. Recall that the division is well-defined since c is invertible in K[[x]
] by Hypothesis (H). The following lemma will lead us to an efficient algorithm for the latter test.

(G) divided by F in K[[x]][y], so that we have D(G) = QF + R, with deg y (R) ≤ d y -1. Then F divides D(G) in K[[x]][y] if, and only if, [Q] 3dx 2dx = ⌈R⌉ 3dx = 0. Proof. Let us assume that F divides D(G) in K[[x]][y]. Then is it clear that R = 0 and D(G) = QF hold in K[[x]][y]. Since D(G) and F are polynomials then Q ∈ K(x)[y] and D(G) = QF holds in K(x)[y]. Let a ∈ K[x] and A ∈ K[x, y] be such that Q = A/a, so that we can write a D(G) = AF in K[x][y]. Since F is primitive by Hypothesis (H), the classical Gauss lemma [Lan02, Chapter IV, Theorem 2.1] implies that a divides the content of A. It follows that Q actually belongs to K[x, y], whence deg x (Q) ≤ 2d x -1. Conversely, let us assume that [Q] 3dx 2dx = ⌈R⌉ 3dx = 0. Then D(G) = ⌈Q⌉ 2dx F holds in K[[x]]/(x 3dx )[y]. Since deg x (⌈Q⌉ 2dx F ) ≤ 3d x -1, we deduce that D(G) = ⌈Q⌉ 2dx F holds in K[x, y], hence that F divides D(G) in K[[x]][y].
For any p > 0 we write K[[x p ]] for the series algebra in x p over K. When p = 0, by convention we let K[[x 0 ]] := K. Lemma 3 motivates us to define the following map:

D F : F s → K[x, y] dx-1,2dy-3 × K[x, y] 3dx-1,dy-1 (ℓ 1 , . . . , ℓ s ) → ([Q] 3dx 2dx /x 2dx , ⌈R⌉ 3dx ). Proposition 1. The inclusion µ 1 , . . . , µ r F ⊆ ker(D F ) always holds. Conversely, if (ℓ 1 , . . . , ℓ s ) belongs to ker(D F ) then ρ 1 , . . . , ρ dy all belong to K[[x p ]].
In addition, if (C) holds then µ 1 , . . . , µ r is the reduced echelon basis of ker(D F ).

Proof. Lemma 3 reformulates as follows: (ℓ 1 , . . . , ℓ s ) belongs to ker(D F ) if, and only if, ρ 1 , . . . , ρ dy all belong to K[[x p ]]. If (ℓ 1 , . . . , ℓ s ) = µ j then we have that G = Fj ∂Fj ∂y by (2), hence that the ρ j belong to {0, 1}, whence the inclusion µ 1 , . . . , µ r F ⊆ ker(D F ).

Let (ℓ 1 , . . . , ℓ s ) ∈ ker(D F ) and assume that (C) holds. In order to prove that (ℓ 1 , . . . , ℓ s ) ∈ µ 1 , . . . , µ r F , it is sufficient to prove that all the ρ i belong K, by Lemmas 1 and 2. The case when p = 0 is immediate so that we can now assume that p ≥ d x (2d y -1)+1. Here we could directly invoke [Gao03, Lemma 2.4] in order to conclude the proof, but, for completeness, let us briefly repeat the arguments. Let i ∈ {1, . . . , d y } and let A ∈ K[x, y] denote the unique irreducible factor of F in K[x, y] such that A(x, φ i )=0. By construction the resultant

B := Res y A, ρ i (0) ∂F ∂y -G ∈ K[x]
has degree at most d x (2d y -1) and belongs to (

x p ), hence it is zero since p ≥ deg(B) + 1. It follows that A divides ρ i (0) ∂F ∂y -G in K[[x]][y], whence ρ i (0) = G(x, φ i )/ ∂F ∂y (x, φ i ) = ρ i . 1.4. Small Positive Characteristic.
Proposition 1 reduces the recombination problem to linear algebra under Hypothesis (C). Until the end of this section, we focus on the situation when (C) does not hold, and F exclusively denotes the prime field F p of K. We shall use the following F-linear maps, respectively reminiscent of Berlekamp's and Niederreiter's univariate factorization algorithms (for instance, see [GG03, Chapter 14]):

B : K[x, y] dx,dy-1 → K(x)[y]/(F ) N : K[x, y] dx,dy-1 → K[x, y p ] pdx,dy-1 G → G p - ∂F ∂y p-1 G G → G p + ∂ p-1 ∂y p-1 F p-1 G .
Here K[x, y p ] pdx,dy-1 is to be understood as the space of the polynomials in x and y p of degrees at most pd x in x and at most d y -1 in y p . One can easily see that N is well defined since ∂ ∂y N(G) = 0. For convenience, Appendix A gathers the classical properties needed in the sequel. As for D, we will be interested in the kernels of the following F-linear maps:

B F : F s → K(x)[y]/(F ) N F : F s → K[x, y p ] pdx,dy-1 (ℓ 1 , . . . , ℓ s ) → B s i=1 ℓ i G i (ℓ 1 , . . . , ℓ s ) → N s i=1 ℓ i G i .
We carry on with the notation: (ℓ 1 , . . . , ℓ s ) denotes a vector in F s , G := s i=1 ℓ i G i , and ρ 1 , . . . , ρ dy represent the residues of G/F as defined in (3).

Proposition 2. µ 1 , . . . , µ r is the reduced echelon basis of ker(B F ) = ker(N F ).

Proof. Proposition 11 of Appendix A applied with F seen in K(x)[y] and E := K((x)) (the fraction field of K[[x]]) gives us that ker(B F ) = ker(N F ), and that (ℓ 1 , . . . , ℓ s ) belongs to these kernels if, and only if, all the ρ i belong to F. The conclusion thus follows from Lemmas 1 and 2. Proposition 2 does not directly provide us with a satisfactory recombination algorithm because the linear system to be solved involves a number of equations that grows linearly with p. The key idea to cut down this dependency in p is to combine N F with D F as follows:

Lemma 4. If (ℓ 1 , . . . , ℓ s ) ∈ ker(D F ) then N F (ℓ 1 , . . . , ℓ s ) belongs to K[x p , y p ] dx,dy-1 . Proof. If (ℓ 1 , . . . , ℓ s ) ∈ ker(D F ) then ρ 1 , . . . , ρ dy all belong to K[[x p ]] by Proposi- tion 1. Therefore equality (8) of Appendix A.1 rewrites into N(G) = F p   dy i=1 ρ p i -ρ i (y -φ i ) p   , whence N F (ℓ 1 , . . . , ℓ s ) = N(G) ∈ K[x p , y p ] dx,dy-1 .
Let ∆ y (x) represent the discriminant of F in y, that is Res y F (x, y), ∂F ∂y (x, y) . We are now ready to present the test for the vanishing of N(G) and B(G) that will be used in the algorithms: Proposition 3. Let S be a set containing at least d x + 1 points in K and assume that (ℓ 1 , . . . , ℓ s ) ∈ ker(D F ). Then the following equivalence holds:

N F (ℓ 1 , . . . , ℓ s ) = 0 ⇐⇒ ∀a ∈ S, N F (ℓ 1 , . . . , ℓ s )(a, y) = 0.
In addition, if c(a)∆ y (a) = 0 for all a ∈ S, then the latter equivalence holds if we replace N F by B F .

Proof. Since the map z → z p is injective in K, the vanishing of N(G) can be tested with only d x + 1 distinct specializations of x by Lemma 4. By Proposition 11 of the appendix applied with F seen in K(x)[y], the equality B(G) = 0 is equivalent to N(G) = 0. Finally by the same proposition 11 applied with F (a, y) seen in K[y], the equality B(G)(a, y) = 0 is equivalent to N(G)(a, y) = 0, for any a ∈ K satisfying c(a)∆ y (a) = 0.

Recombination Algorithms

In this section, we present a deterministic recombination algorithm and a faster probabilistic one both derived from the results of the previous section. We still assume that the normalization Hypothesis (H) holds. The subfield F is set to be K if (C) holds and to be the prime field F p of K otherwise.

2.1. Complexity Model. For our cost analysis, we use the computation tree model [BCS97, Chapter 4] from the total complexity point of view. Roughly speaking, this means that complexity estimates charge a constant cost for each arithmetic operation (+, -, ×, ÷) and the equality test. Yet all the constants in the base fields (or rings) of the trees are thought to be freely at our disposal.

2.1.1. Polynomial Arithmetic. A univariate polynomial of degree d is thought to be represented as the vector of its coefficients of size d + 1. For each integer d, we assume that we are given a computation tree that computes the products of two polynomials of degree at most d with at most M(d) operations, independently of the base ring. As in [GG03, Chapter 8.3], for any positive integers d 1 and d 2 , we assume that M satisfies the following properties:

M(d 1 d 2 ) ≤ d 2 1 M(d 2 ) and M(d 1 )/d 1 ≤ M(d 2 )/d 2 if d 1 ≤ d 2 . In particular, this implies the super-additivity of M, that is M(d 1 ) + M(d 2 ) ≤ M(d 1 + d 2 ).
During the cost analyzes we will appeal to the following classical results:

• The resultant and the extended greatest common divisor of two univariate polynomials of degree at most d over a field K can be computed with If A is a commutative ring with unity, and if q ∈ A[z] is a monic polynomial of degree d then an element of A[z]/(q(z)) is represented by its coordinate vector in the usual basis 1, . . . , z d-1 . Each arithmetic ring operation in A[z]/(q(z)) reduces to O(M(d)) operations in A.

O(M(d) log(d)) operations in K [GG03,
2.1.2. Linear Algebra. Concerning linear algebra, we assume that, for each n, we are given a computation tree that computes the product of two n × n matrices with at most O(n ω ) ring operations (i.e. without inversion nor division), for a fixed constant ω. We require that 2 < ω ≤ 3 in order to use the following statement: Lemma 5. [Sto00, particular case of Proposition 2.11] The computation of the row reduced echelon form of a m × n matrix over K of rank r takes O(mnr ω-2 ) operations in K.

Corollary 1. If m ≥ n then the computation of the reduced echelon basis of the kernel of a m × n matrix over K takes O(mn ω-1 ) operations in K.

Proof. Let M be such a m × n matrix. Let M denote the m × n matrix row mirrored of M , that is Mi,j = M i,n-j+1 for all (i, j) ∈ {1, . . . , m} × {1, . . . , n}. By the previous proposition, the computation of the reduced echelon form of M costs O(mn ω-1 ). Then the reduced echelon basis of the kernel of M can be directly read off from the latter echelon form with O(n 2 ) operations in K.

2.1.3. Extensions for our Algorithms. In order to present a complete factorization algorithm in Section 3, we need to enlarge the computational model with irreducible factorization in K[y]. Separately we also count the number of pth root extractions in K. By root extraction we mean testing if an element in K is a pth root, and returning its casual root.

When (C) does not hold we further need to know the algebra structure of K over its prime field F. Precisely we require that we can access to any component of any element in some basis over F for free. We also assume that any element of F can be sent into K for free. Our computational trees will thus contain operations both in K and F p . These assumptions are indeed not restrictive and cover the case when K is explicitly finitely generated over F p , that is considered in [START_REF] Davenport | Factorization over finitely generated fields[END_REF][START_REF] Steel | Conquering inseparability: primary decomposition and multivariate factorization over algebraic function fields of positive characteristic[END_REF].

For the cost analysis in positive characteristic, we will mainly focus on the case when K is a finite field, that K := F p k for some k ≥ 1. Such a finite field K is supposed to be given as a quotient F p [z]/q(z) with q being an irreducible polynomial of degree k, so that each field operation in K amounts to O(M(k) log(k)) operations in F p . Over finite fields, we will exclusively use computation trees over F p .

2.1.4. Probabilistic Algorithms. In order to modelize probabilistic algorithms, we extend the computation tree model with a function that takes an integer n as input, and that returns a random subset of a fixed set N of cardinality n, assuming that the cardinality N of N is at least n. The cost of this operation is assumed to be bounded by a super-additive function written R(n), that only depends on n. The probability distribution is supposed to be uniform in the space of subsets of N of cardinality n. Let us recall the following probability estimate: Lemma 6. [Lec07a, Lemma 9] Let M be a subset of N of cardinality M . For any n ≤ M , the density of subsets of N of cardinality 2n having at most n elements in M is at least E(M, N ), where: 7. Compute a solution basis ν 1 , . . . , ν r of the following linear system in the unknowns (ℓ 1 , . . . , ℓ t ) ∈ F t : Proof. If (C) holds then Proposition 1 gives us that (μ t , . . . , μt ) = (µ 1 , . . . , µ r ). If (C) does not hold then Proposition 1 gives us that µ 1 , . . . , µ r F ⊆ μ1 , . . . , μt F , so that µ 1 , . . . , µ r F is the kernel of the restriction of N F to μ1 , . . . , μt F by Proposition 2. Since e is chosen so that z e -1 admits e ≥ d x + 1 distinct roots in K, Proposition 3 implies that ν 1 , . . . , ν r is a basis of the kernel of the latter restriction expressed in the coordinates μ1 , . . . , μt . Therefore the reduced echelon basis returned in step 8 is really µ 1 , . . Let us now assume that (C) does not hold. The resolution of the linear system in step 4 now takes O(max(e D , s)s ω-1 ) operations in F by Corollary 1. In step 5, the computation of each Gi can be done by means of the sub-product tree technique with O(M (d x )M(d y ) log(s)) operations in K.

E(M, N ) := 1 1 + M N -M 2 , for N ≥ 2M := 0, for N < 2M.
t i=1 ℓ i N( Gi )(ζ, y) = 0. ( 5
By Proposition 13 the cost of step 6 belongs to

O M(d x ) M(d y )(d y + log(p)) + d 2 y t ω-2 + td y (log(d x ) + log(p)) .
The resolution of the linear system in step 7 costs O(max(e N , t)t ω-1 ) operations in F by Corollary 1. In step 8 the matrix product and the reduced form computation amount to O(s ω ) operations in F thanks to Proposition 5. 2.3. Probabilistic Algorithm. When (C) holds, when s is close to d y , and when ω is close to 3, the first bottleneck of Algorithm 1 is the resolution of the linear system (4) that amounts to O(d x d ω y ) operations in K. When (C) does not hold, and when t is close to d y , the second bottleneck is the evaluation of N is step 6 that also contributes to Õ(d x d ω y ).

In this section we improve Algorithm 1 by means of probabilistic techniques. Concerning the first bottleneck we appeal to the following classical method for overdetermined linear systems. From now on we let m := 5d x d y -2d x and n := (d x + 1)d y . If (u 2 , . . . , u m ) ∈ K m-1 , we write U for the following upper triangular s × m Toeplitz matrix:

U :=      1 u 2 u 3 • • • u m-1 u m 1 u 2 u 3 • • • u m-1 . . . . . . . . . . . . 1 . . . u m-s+1      . Similarly, if (v 2 , . . . , v n ) ∈ K n-1
, we write V for the following upper triangular t × n Toeplitz matrix:

V :=      1 v 2 v 3 • • • v n-1 v n 1 v 2 v 3 • • • v n-1 . . . . . . . . . . . . 1 . . . v n-t+1      .
The following lemma, derived from [KS91, Theorem 2], will provide us with a probabilistic algorithm to solve the overdetermined linear systems (4) and ( 5):

Lemma 7. [Lec07b, Lemma 9] Let A be a m×s matrix over K. For all (u 2 , . . . , u m ) in K m-1 we have ker(A) ⊆ ker(U A). In addition, there exists a nonzero polynomial P ∈ K[z 2 , . . . , z m ] of total degree at most s such that the latter inclusion is an equality whenever P(u 2 , . . . , u m ) = 0.

Concerning the second bottleneck we replace the operator N by B that can be evaluated faster. The probabilistic algorithm thus proceeds as follows:

Algorithm 2. Probabilistic recombination Input: F ∈ K[x, y], F 1 , . . . , F s to precision (x dx+1 ), and a finite subset N of K of cardinality N . Output: µ 1 , . . . , µ r . Do the same computations as in Algorithm 1, but replace steps 4, 6, and 7 respectively by: 4. a. Take (u 2 , . . . , u m ) at random in N m-1 .

b. For all i ∈ {1, . . . , s}, compute D i := U D(G i ), where D(G i ) is seen as a vector in K m . c. Compute the reduced echelon solution basis μ1 , . . . , μt of the following linear system in the unknowns (ℓ 1 , . . . , ℓ s ) ∈ F s :

s i=1 ℓ i D i = 0. ( 6
)
If (C) holds then return μ1 , . . . , μt . 6. a. Take a finite subset of N at random with cardinality 2(d x + 1).

b. If possible select d x + 1 points a 0 , . . . , a dx from the latter set such that c(a i )∆ y (a i ) = 0. If not possible then stop the execution. c. For all i ∈ {1, . . . , t}, compute the vector A i in K n whose jd y + k + 1th entry is the coefficient of y k in B( Gi )(a j , y), for all k ∈ {0, . . . , d y -1} and all j ∈ {0, . . . , d x }. 7. a. Take (v 2 , . . . , v n ) at random in N n-1 . b. For all i ∈ {1, . . . , t}, compute B i := V A i . Compute a solution basis ν 1 , . . . , ν r of the following linear system in the unknowns (ℓ 1 , . . . , ℓ t ) ∈ F t : Proof. By Lemma 7, the solution space of system (6) always contains the solution space of system (4), that is µ 1 , . . . , µ r . If (C) does not hold then, for any i ∈ {1, . . . , r}, we can write µ i = ℓ 1 μ1 + • • • + ℓ t μt with (ℓ 1 , . . . , ℓ t ) ∈ F t so that we have:

t i=1 ℓ i B i = 0. ( 7
t j=1 ℓ j B( Gj ) = B   t j=1 ℓ j s k=1 μj,k G k   = B s k=1 µ i,k G k = B F (µ i ) = 0,
where the last equality follows from Proposition 2. Therefore, by Lemma 7 again, the space spanned by the output of the algorithm always contains µ 1 , . . . , µ r .

Lemma 7 provides us with a polynomial P such that, for any (u 2 , . . . , u m ) outside the zero locus of P, the systems (4) and ( 6) have the same solution sets. The probability of success thus follows from the Schwartz-Zippel lemma.

If (C) does not hold, then the probability that at most d x + 1 points annihilate c∆ y comes from Lemma 6. Then Lemma 7 provides us with a nonzero polynomial Q such that the solution set of (7) coincides with the solution set of the following system:

t i=1 ℓ i A i = 0.
On the other hand, the solution set of the latter system coincides with the solution set of system (5) by Proposition 3, whence the correctness of Algorithm 2. The probability bound in part (b) follows from Z(s, N ) 2 ≥ Z(2s, N ).

If (C) holds then step 4 performs O(sM(d x d y )) operations in K for the computation of all the D i , plus O(s ω ) operations in K for the resolution of system (6), by Corollary 1. We are done with part (a).

Concerning part (b), the evaluation of c∆ y at 2(d x + 1) in step 6 can be done with O(d y M(d x ) log(d x ) + d x M(d y ) log(d y )). Then F (a 0 , y), . . . , F (a dx , y), and each tuple ( Gi (a 0 , y), . . . , Gi (a dx , y)) can be computed with O(d y M(d x ) log(d x )) operations in K. Each B( Gi )(a j , y) is then obtained as the remainder of Gi (a j , y) p -(F (a j , y) ′ ) p-1

Gi (a j , y) p divided by F (a j , y) with O(M(d y ) log(p)) operations in K. Therefore the total cost of step 6 amounts to O(t(d y M(d x ) log(d x )+d x M(d y ) log(p))) operations in K.

In step 7 the computation of all the B i takes O(tM(d x d y )) operations in K, and the resolution of the linear system costs O(max(f B , t)t ω-1 ) operations in F, by Corollary 1. The total cost follows from summing the costs of each step and taking into account that p ∈ O(d x d y ).

Corollary 3. If K := F p k then Algorithm 2 performs Õ(k(d x d 2 y + d ω y )) operations in F p , plus R(O(d x d y )).
The probability to obtain a correct answer is at least E(2d x d y , N )Z(2d y , N ). 2.4. Testing the Recombination and Recovering the Factors. In this subsection we describe how to test whether the output of Algorithm 2 is correct or not. If correct then we also recover the irreducible factors of F . If not correct then the execution is stopped. Let ν 1 , . . . , ν t now represent the output of Algorithm 2. From Proposition 5 we know that µ 1 , . . . , µ r F ⊆ ν 1 , . . . , ν t F . Algorithm 3. Recovering the irreducible factors.

Input: F ∈ K[x, y], F 1 , . . . , F s to precision (x dx+1 ), and a reduced echelon family ν 1 , . . . , ν t in F s such that µ 1 , . . . , µ r F ⊆ ν 1 , . . . , ν t F . Output: F 1 , . . . , F r .

1. If t = 1 then return F . 2. If the entries of the ν i are not all in {0, 1} then stop the execution. 3. If the supports of ν 1 , . . . , ν t do not form a partition of {1, . . . , s} of size t then stop the execution.

4. For all i ∈ {1, . . . , t}, compute Fi := c s j=1 F νi,j j dx+1 .

5. For all i ∈ {1, . . . , t}, compute Fi as the primitive part of Fi in y. Normalize Fi in K[x][y] so that its leading coefficient becomes monic. Proof. Let us first examine the case when t = r. If r = 1 then the output is clearly correct. If r ≥ 2 then, for all i ∈ {1, . . . , r}, we have that Fi

If

= c/c i F i because deg x (c/c i ) + deg x (F i ) ≤ deg x (F 1 ) + • • • + deg x (F r ) = d x .
Finally Fi = F i holds and the algorithm returns a correct answer.

Let us now show that the algorithm always return a correct result. If t = 1 then the algorithm exits at the first step with t = r. If the algorithm exits in step 2 or 3 then we necessary have t > r. If the algorithm reaches the last step then F admits at least t ≥ r irreducible factors, whence t = r.

Steps 1 to 3 take O(st) operations in K. Letting

The Whole Factorization Algorithm

So far we have dealt with the factorization of normalized polynomials, namely under Hypothesis (H). In this section we explain how to handle the general case by means of successive reductions to (H). For the sake of simplicity, we use fast multiplication everywhere, namely M(d) ∈ Õ(d).

3.1. Factorization of Separable Polynomials. Recall that a polynomial F is said to be separable in y if its discriminant ∆ y , is nonzero. The reduction of the factorization of separable polynomials to normalized polynomials simply consists in searching for a suitable shift of the variable x. The following algorithm summarizes the main steps: Proof. The proof is the same as for the deterministic case, except that we appeal to Propositions 5 and Corollary 3.

3.2. The General Case. In order to raise the separability assumption on F , we use the separable factorization algorithm designed in [START_REF] Lecerf | Fast separable factorization and applications[END_REF]. Let B := {1, p, p 2 , p 3 , . . .} be the set of the powers of p. If F is primitive in y, and if p > 0, then the separable decomposition of F in y is defined to be the set

{(G 1 , q 1 , m 1 ), . . . , (G s , q s , m s )} ⊆ (K[x, y] \ K[x]) × B × N
satisfying the following properties:

(S 1 ) F (y) = s i=1 G i (y qi ) mi ; (S 2 ) for all i = j in {1, . . . , s}, G i (y qi ) and G j (y qj ) are coprime; (S 3 ) for all i ∈ {1, . . . , s}, p does not divide m i ; (S 4 ) for all i ∈ {1, . . . , s}, G i is separable, primitive and of positive degree in y; (S 5 ) for all i = j in {1, . . . , s}, (q i , m i ) = (q j , m j ). If p = 0 then the separable decomposition of F is naturally defined to be the set generated by all the triples (G, 1, m) such that G is a proper squarefree factor of F with multiplicity m ≥ 1. For the existence, the uniqueness, and a short history of the separable decomposition we refer the reader to [START_REF] Lecerf | Fast separable factorization and applications[END_REF].

Under the assumption that K has sufficiently many elements we finally obtain the following top level factorization algorithm:

Algorithm 5. Top level factorization algorithm.

Input:

F ∈ K[x, y].
Output: the irreducible factors F 1 , . . . , F r of F together with their respective multiplicities e 1 , . . . , e r . 1. Compute the content and the primitive part of F in y. Replace F by its primitive part and initialize the list L with the irreducible decomposition of the content. 2. Compute the separable decomposition (G 1 , q 1 , m 1 ), . . . , (G s , q s , m s ) of F . 3. Compute the irreducible decomposition of G 1 , . . . , G s . 4. For all i ∈ {1, . . . , s}, and for all irreducible factor E of G i do:

If p = 0 then append (E, m i ) to the list L, else compute (H, h) ∈ K[x, y] × B such that H h (y qi/h ) = E(y qi ), h ≤ q i and h is as large as possible, and append (H(y qi/h ), hm i ) to L. 5. Return L.

By means of deterministic subroutines we obtain: Proposition 9. Assume that K has cardinality at least 2d x d y + d x + 1. Then Algorithm 5 performs irreducible factorizations of univariate polynomials over K whose degree sum is at most d Proof of Theorem 2. We also reduce to the case when d y ≤ d x so that we have Õ(d x d 2 y + d ω y ) ⊆ Õ((d x d y ) 1.5 ). Thanks to the hypothesis on the cardinality of K we can take N with cardinality N = 10d x d y so that E(4d x dy, N )Z(2d y , N ) > 1/2 holds. The conclusion follows from Proposition 10.

x + d y , plus a. if (C) holds, O(d x d ω y ) arithmetic operations in K; b. if K := F p k , Õ(kd x d ω y )

Conclusion

Although the factorization algorithms presented in this paper have good worst case complexity bounds their implementations require some care to make them very efficient in practice. The first bottleneck of the algorithm underlying Theorem 2 is the computation of the polynomials G i . This bottleneck can be avoided by adapting the heuristic presented in [BLS + 04, [START_REF] Lecerf | Sharp precision in Hensel lifting for bivariate polynomial factorization[END_REF][START_REF] Lecerf | Improved dense multivariate polynomial factorization algorithms[END_REF]. Then the next practical bottleneck really becomes the Hensel lifting. Although this lifting can be done in softly optimal time, one interesting speedup consists in solving the recombination problem progressively at each stage of the lifting in order to decrease the number of the lifted factors and stop the lifting as soon as possible. Finally, in a good implementation, the worst case precision bound for the lifting is hardly never attained. These heuristics and others have been implemented by Steel in Magma [Mag] on the top of the algorithms designed in [START_REF] Van Hoeij | Factoring polynomials and the knapsack problem[END_REF][START_REF] Belabas | Factoring polynomials over global fields[END_REF][START_REF] Steel | Conquering inseparability: primary decomposition and multivariate factorization over algebraic function fields of positive characteristic[END_REF]. Our new algorithms are not intended to compete in general with Steel's implementation. Instead they sensibly improve the performances in very particular cases. Such cases are very difficult to build. One example is provided by [BHKS04, Remark 5.5]. b. g p -(f ′ ) p-1 g ∈ (f ) (Berlekamp's operator). c. g p + (f p-1 g) (p-1) = 0 (Niederreiter's operator).

Proof. In short let ρ i represent g(ϕ i )/f ′ (ϕ i ) for each i ∈ {1, . . . , d}, so that the partial fraction decomposition of g/f over E reads as:

g f = d i=1 ρ i x -ϕ i .
For any i ∈ {1, . . . , d} we have that (

g p -(f ′ ) p-1 g)(ϕ i ) = (f ′ ) p (ϕ i )(ρ p i -ρ i )
, which gives the equivalence between parts (a) and (b). The equivalence between parts (a) and (c) follows from the following calculation:

g p + (f p-1 g) (p-1) = g p + f p g f (p-1) = g p + f p g f (p-1) = f p g f p + g f (p-1) = f p d i=1 ρ p i (x -φ i ) p + d i=1 -ρ i (x -φ i ) p , (8)
where the latter equality makes use of Wilson's theorem (namely, p divides (p -1)! + 1).

A.2. Complexity. Let A be any commutative ring with unity of characteristic p, and let f be a polynomial in A[x] of degree d. We are to study the cost for computing g p + (f p-1 g) (p-1) for a given g ∈ A[x] d-1 . In the latter expression, the most difficult part is the computation of (f p-1 g) (p-1) . We thus focus on the construction of the d × d matrix N of the following map expressed in the usual monomial basis:

A[x] d-1 → A[x p ] d-1
g → (f p-1 g) (p-1) .

Recall that A[x p ] d-1 denotes the space of polynomials in x p of degree at most d -1 in x p .

Proposition 12. If f is monic then N can be computed with O(M(d)(d + log(p))) operations in A.

Proof. We write coeff(A, x i ) for the coefficient of x i in A ∈ A[x]. We aim to computecoeff(f p-1 g, x ip+p-1 ) for all i ∈ {0, . . . , d -1}. Let rev(n, .) denote the reversal endomorphism of A[x] n defined by rev(n, A) = n k=0 a n-k x k , for all A := n k=0 a k ∈ A[x] n . In fact we will compute coeff(rev(dp -1, f p-1 g), x ip ), for all i ∈ {0, . . . , d -1}.

Since f is monic of degree d, the constant term of rev(d, f ) is 1. Hence we can define i≥0 u i x i as the power series expansion of rev(d -1, g)/ rev(d, f ), so that we have: . In order to compute u 0 , u p , . . . , u (d-1)p , we are going to use the fact that (u i ) i∈N satisfies the following linear recurrence:

u i := -(f d-1 u i-1 + • • • + f 0 u i-d ), for i ≥ d,
where f j := coeff(f, x j ) (recall that f d = 1). Of course the initial values u 0 , . . . , u d-1 can be computed from the power series expansion of rev(d -1, g)/ rev(d, f ) to precision (x d ) with O(M(d)) operations in A.

For all i ≥ 0, if a d-1 x d-1 + • • • + a 1 x + a 0 represents the remainder rem(x i , f ) in the division of x i by f , then it is classical that u i = a 0 u 0 + • • • + a d-1 u d-1 . The idea of using the latter formula to efficiently compute u i is due to Fiduccia [Fid85]. This formula expresses the fact that linear recurrence sequence extension is transposed to the remainder operation (see [BLS03, Section 5] for details). Let B be the d × d matrix whose ith column is rem(x p(i-1) , f ) expressed in the usual basis 1, x, . . . , x d-1 . This is nothing else but the Petr-Berlekamp matrix arising in univariate polynomial factorization over finite fields. The construction of B takes O(M(d)(d+log(p))) operations in A. Finally we have shown that N can be factored into N = H 2 B t H 1 , where B t is the transpose of B, and H 1 , H 2 are the Hankel matrices of the following maps: Recall that ⌈A⌉ l represents the projection of the series A := i≥0 a i x i ∈ K[[x]] to the polynomial l-1 i=0 a i x i ∈ K[x] l-1 . Since H 1 is symmetric, (B t H 1 ) t = H 1 B can be computed with O(dM(d)) operations in A. The second product H 2 (B t H 1 ) also costs O(dM(d)).

H 1 :A[x] d-1 → A d
The above decomposition N = H 2 B t H 1 was already stated in [Nie94, Theorem 1] (in the more general setting of the Hasse-Teichmüller derivative). Unfortunately the cost analysis made in [Nie94, Theorem 2] was pessimistic because it did not take into account that left and right products by Hankel matrices could be done in softly optimal time.

From now on we assume that A := K[z]/(q(z)), where q is a non-constant monic separable polynomial in K[z] of degree n. Let g 1 , . . . , g t be a sequence of polynomials in A[x] d-1 with 1 ≤ t ≤ d. We are interested in evaluating Niederreiter's operator simultaneously at each g i .

Proposition 13. The computation of g p 1 + (f p-1 g 1 ) (p-1) , . . . , g p t + (f p-1 g t ) (p-1) takes O M(n) M(d)(d + log(p)) + d 2 t ω-2 + td(log(n) + log(p))

operations in K.

Proof. The computation of g p 1 , . . . , g p t can be done with O(td log(p)) ring operations in A, which amounts to O(tdM(n) log(p)) operations in K. We can now focus

Lemma 3 .

 3 Let Q and R respectively denote the quotient and remainder of D

  Chapter 11]. • The product of r univariate polynomials G 1 , . . . , G r over a field K whose degree sum is d takes O(M(d) log(r)) operations in K by means of the subproduct tree technique [GG03, Chapter 10]. • If F ∈ K[y] has degree d then the remainders of F modulo all the G i can also be computed with O(M(d) log(r)) operations in K. This task is usually called the simultaneous reduction. The inverse problem, called Chinese remaindering, also costs O(M(d) log(r))[GG03, Chapter 10].

  Note that E(M, N ) ≥ 1/2 whenever N ≥ 2M . Let us finally recall the classical Schwartz-Zippel lemma: the density of points in N that do not annihilate a given nonzero polynomial in n variables of degree d is at least Z(d, N ) := 1d/N [GG03, Chapter 6, Section 9]. 2.2. Deterministic Algorithm. Propositions 1, 2 and 3 naturally lead to the following recombination algorithm: Algorithm 1. Deterministic recombination. Input: F ∈ K[x, y], and F 1 , . . . , F s to precision (x dx+1 ). Output: µ 1 , . . . , µ r . 1. For each i ∈ {1, . . . , s}, compute Fi as the quotient of F by F i to precision (x dx+1 ). 2. Compute F1 ∂F1 ∂y , . . . , Fs ∂Fs ∂y to precision (x dx+1 ) and deduce G 1 , . . . , G s . 3. Compute D(G 1 ), . . . , D(G s ). 4. Compute the reduced echelon solution basis (μ 1 , . . . , μt ) of the following linear system in the unknowns (ℓ 1 , . . . , ℓ s ) ∈ F s : s i=1 ℓ i D(G i ) = 0. (4) If (C) holds then return μ1 , . . . , μt . 5. For all i ∈ {1, . . . , t}, compute Gi = s j=1 μi,j G j . 6. If p divides d x + 1 then let e := d x + 2 else let e := d x + 1. Let ζ denote the residue class of z in K[z]/(z e -1). Compute N( G1 )(ζ, y), . . . , N( Gt )(ζ, y).

) 8 .

 8 Compute and return the reduced echelon basis of the space t j=1 ν i,j μj | i ∈ {1, . . . , r} F . Proposition 4. Under Hypothesis (H) Algorithm 1 works correctly as specified. a. If (C) holds then Algorithm 1 performs O(d x d y s ω-1 + sM(d x )M(d y )) operations in K. b. If (C) does not hold then Algorithm 1 performs O M(d x ) d 2 y s ω-2 + M(d y )(s log(s) + d y + log(d x d y )) + sd y log(d x d y ) operations in K plus O(max(e D , e N , s)s ω-1 ) operations in F p , where e D and e N represent the number of equations in systems (4) and (5) respectively.

  . , µ r . Let us now analyze the cost of the algorithm. Testing whether (C) holds or not, and obtaining the characteristic p of K when (C) does not hold, takes O(d x d y ) operations in K. Steps 1 to 3 take O(sM(d x )M(d y )) operations in K. If (C) holds then the linear system in step 4 involves s unknowns and O(d x d y ) equations. Therefore its resolution over K takes O(d x d y s ω-1 ) operations in K by Corollary 1. We are done with part (a).

  Finally part (b) follows from summing the costs of all steps, and taking into account that log(p) ∈ O(log(d x d y )). Corollary 2. If K := F p k then Algorithm 1 performs Õ(kd x d ω y ) operations in F p . Proof. This follows from the previous proposition since e D ∈ O(kd x d y ) and e N ∈ O(kd x d y ).

) Proposition 5 .

 5 The space spanned over F by the output of Algorithm 2 always contains µ 1 , . . . , µ r .a.If (C) holds then Algorithm 2 takes O(s(M(d x )M(d y ) + M(d x d y )) + s ω ) operations in K, plus R(m -1). The probability to obtain a correct result is at least Z(s, N ). b. If (C) does not hold then Algorithm 2 takes O(s(M(d x )M(d y ) log(d x d y ) + M(d x d y )) operations in K, plus R(2d x + m + n), plus O(max(f D , f B , s)s ω-1 )operations F p , where f D and f B represent the number of equations in systems (6) and (7) respectively. The probability to obtain a correct result is at least E(2d x d y , N )Z(2s, N ).

Proof.

  It suffices to set f D ∈ O(ks) and f B ∈ O(ks) in Proposition 5.

ti=1

  Fj = F then return F1 , . . . , Ft . Otherwise stop the execution. Proposition 6. Algorithm 3 returns a correct answer if, and only if, t = r. Otherwise the execution is stopped. The algorithm performs O(M(d x )M(d y ) log(d x d y ) + s 2 ) arithmetic operations in K.

  d i := deg y ( Fi ), each Fi can be computed by the sub-product tree technique with O(M(d x )M(d i ) log(d i )) operations, so that the total cost of step 4 amounts to O(M(d x )M(d y ) log(d y )). In step 5 each primitive part computation amounts to O(d i M(d x ) log(d x )) operations, whence a total cost in O(d y M(d x ) log(d x )) for this step. In step 6, one has to check whether t i=1 deg x ( Fi ) = d x holds or not before computing the product. This way the cost of step 6 drops to another O(M(d x )M(d y ) log(d y )).

Algorithm 4 .

 4 Factorization of separable polynomials. Input: F ∈ K[x, y] primitive and separable in y. Output: the irreducible factors F 1 , . . . , F r of F . 1. Find b ∈ K such that F (b, y) is separable of degree d y . Replace F by F (x + b, y) 2. Compute the irreducible decomposition of F (0, y). 3. Compute F 1 , . . . , F s to precision (x dx+1 ) by means of Hensel lifting. 4. Recombine the lifted factors by means of the algorithms presented in the previous section in order to obtain the irreducible factors F 1 , . . . , F r of F over K. 5. For all i ∈ {1, . . . , r} replace F i by F i (xb, y).With deterministic subroutines we obtain the following estimates: Proposition 7. Assume that F is primitive and separable in y, and that K has cardinality at least 2d x d y + 1. Algorithm 4 is correct and performs one irreducible factorization inK[y] in degree d y , plus a. if (C) holds, O(d x d ω y ) operations in K; b. if K := F p k , Õ(kd x d ω y ) operations in F p . Proof.Since deg(c∆ y ) ≤ 2d x d y , the assumption on the cardinality of K ensures us to find a suitable b. This search can be done with Õ(d x d 2 y ) operations by means of fast multi-point evaluation. Then the shift of the variable x can be done in softly optimal time via evaluation/interpolation.In steps 2 to 4 Hypothesis (H) holds. The computation of F 1 , . . . , F s to precision (x dx+1 ) from the irreducible factors of F (0, y) costs Õ(d x d y ) operations in K by[START_REF] Zur Gathen | Modern computer algebra[END_REF] Theorem 15.18] (see also [BLS + 04] for implementation details). Then part (a) follows from Propositions 4 and 6. Then part (b) is completed thanks to Corollary 2.As for the probabilistic version, we call Algorithm 2 in step 4. For doing so we need to provide the algorithm with a subset N of K of cardinality N , used for random element generation. With these modifications we obtain: Proposition 8. Assume that F is primitive and separable in y, and that K has cardinality at least 2d x d y + 1. Algorithm 4 either stop prematurely or returns a correct answer. It takes one irreducible factorization in K[y] in degree d y , plus a. if (C) holds, Õ(d x d 2 y +d ω y ) operations in K and R(O(d x d y )), with a probability of success at least Z(d y , N ); b. if K := F p k , Õ(k(d x d 2 y + d ω y )) operations in F p and R(O(d x d y )), with a probability of success at least E(2d x dy, N )Z(2d y , N ).

  arithmetic operations in F p . Proof. In Algorithm 5, the computation of the primitive part can be done in softly optimal time. Then step 2 takes Õ(d x d 2 y ) by [Lec07a, Proposition 8] thanks to the hypothesis on the cardinality of K. Step 3 can be done by means of Algorithm 4: Proposition 7(a) yields a total cost in O(d x d ω y ) operations in K. We are done with part (a). As for part (b), computing each pair (H, h) costs O(deg x (E) deg y (E)) extractions of pth roots in F p k . One pth root extraction in F p k requires O((k -1) log(p)) arithmetic operations in F p k . The conclusion thus follows from Proposition 7(b).In order to use probabilistic subroutines in Algorithm 5 we add an extra input N that is a subset of K of cardinality N used for the random choices. Proposition 10. Assume that K has cardinality at least 2d x d y + 1. Then Algorithm 5 either returns a correct answer or stops prematurely. It performs irreducible factorizations of univariate polynomials over K whose degree sum is at most d x +d y , plus a. if (C) holds, O(d x d 2 y + d ω y ) arithmetic operations in K, plus R(O(d x d y )) with a probability of success at least E(2d x d y , N )Z(d y , N ); b. if K := F p k , Õ(k(d x d 2 y + d ω y )) arithmetic operations in F p , plus R(O(d x d y )) for random set generation in K, with a probability of success at least E(4d x dy, N )Z(2d y , N ). Proof. Step 2 of Algorithm 5 takes Õ(d x d y ) operations in K plus R(O(d x d y )), with a probability of success at least E(2d x d y , N ) by [Lec07a, Proposition 9]. The proof thus follows as for the previous proposition thanks to Proposition 8. In part (b) the probability bound makes use of E(2d x dy, N ) 2 ≥ E(4d x dy, N ) [Lec07a, Lemma 12].Proof of Theorem 1. Without loss of generality we can assume that d y ≤ d x so that we have Õ(d x d ω y ) ⊆ Õ((d x d y ) (ω+1)/2 ). The conclusion thus follows from Proposition 9.

  rev(dp -1, f p-1 g) = rev(d, f ) p-1 rev(d -1, g) = rev(d, f ) p rev(d -1, g)/ rev(d, f ) = rev(d, f ) p i≥0 u i x i .Therefore we are led to compute u 0 , u p , . . . , u (d-1)p . The computation of rev(d, f ) p takes O(d log(p)) operations in A, and then the product rev(d, f ) p d-1 i=0 u ip x ip costs O(M(d))

g

  → (u 0 , . . . , u d-1 ),H 2 :A d → A[x p ] d-1 (u 0 , . . . , u (d-1)p ) →rev   dp -1, rev(d, f )

This work was partly supported by the French Research Agency via the Gecko project (http://gecko.inria.fr).

Appendix A. Berlekamp's and Niederreiter's Equations

This appendix gathers some classical properties of Berlekamp's and Niederreiter's operators used in Section 1, and contains a cost analysis of the evaluation of Niederreiter's operator used in Section 2. From now on K denotes a field of characteristic p > 0, whose prime field is written F.

A.1. Classical Properties. Let f be a monic separable polynomial of degree d in K[x], and let g ∈ K[x] be of degree at most d -1. Let E be a field extension of K containing all the roots ϕ 1 , . . . , ϕ d of f . The kth derivative of h ∈ K[x] is written h (k) .

Proposition 11. The following assertions are equivalent: a. g(ϕ i )/f ′ (ϕ i ) ∈ F, for all i ∈ {1, . . . , d}.

on the computation of (f p-1 g 1 ) (p-1) , . . . , (f p-1 g t ) (p-1) . If f is monic then this computation can be achieved by means of the matrix N of Proposition 12 that can be built with O(M(d)(d + log(p))) operations in A: it suffices to calculate the product of N with the d × t matrix whose columns are the coefficients of the g i . This product takes O(d 2 t ω-2 ) operations in A. We are done with the monic case. If f is not monic then we decompose

, where q i will be a monic polynomial in K[z] if degree n i ≥ 1. We perform the computations in each A i and recover the result by means of Chinese remaindering. This decomposition is obtained as follows. For all i ∈ {0, . . . , d}, let f i represent the canonical preimage in K[z] of the coefficient of x i in f . We construct A 1 := K[z]/(q 1 (z)) so that the projection of f into A 1 [x] has an invertible leading coefficient. More precisely, q 1 is obtained as q 1 := q/ gcd(f d , q). Since q is separable, q/q 1 is prime with q 1 . It follows that f d is invertible in A 1 . Let k 1 be the largest integer such that f k1 is nonzero modulo q/q 1 . Of course we have k 1 ≤ d -1, and k 1 can be found with (dk 1 )M(n) operations in K. From f k1 and q/q 1 we can compute q 2 such that f k1 is invertible in A 2 = K[z]/(q 2 (z)) and reduces to zero modulo q/(q 1 q 2 ). By iterating this process we obtain the claimed decomposition with O(dM(n) log(n)) operations in K.

By [GG03, Corollary 10.7] any element of A can be sent to

Therefore all the projections of f and the

Now let i ∈ {1, . . . , r} and let us analyze the computation of

. For each j ∈ {1, . . . , t} we introduce π j and ρ j for the quotient and the remainder of g j divided by f in A i , so that g j = π j f + σ j holds in A i . By construction the divisions are well defined since the leading coefficient of f is invertible in A i . The inversion of the leading coefficient of f in A i takes O(M(n i ) log(n i )) operations in K. Then the computation of all the π j and σ j amounts to O(tM(d)) ring operations in A i . For each j we thus have to compute

The computation of all the (f p-1 σ j ) (p-1) can be done with O(d

The sum of these costs over i ∈ {1, . . . , r} leads to O(M(n)(M(d)(d + log(p)) + d 2 t ω-2 + log(n))) thanks to the super-additivity of M. Finally, all the (f p-1 g i ) (p-1) can be lifted into


Erwan Verron

Roger E Khayat

Abdelouahed Derdouri

Bernard Peseux

Dynamic inflation of hyperelastic spherical membranes

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

While the stability of statically inflated nonlinear materials is well understood, that taking place during dynamic inflation is relatively unexplored. The existence of instabilities in statically inflated hyperelastic membranes is well known, and such instabilities were observed experimentally by Hart-Smith ͑1966͒ and Alexander ͑1971͒. The general equations for static inflation of a thin spherical shell were reviewed by Beatty ͑1987͒, and the stability of the static solution was extensively analyzed recently by Shang ͑1991͒ using bifurcation theory.

More recently, there has also been considerable interest in the dynamic inflation, taking into account the inertia of the membrane. Such investigations are reviewed in detail by Jenkins and Leonard ͑1991͒ and Jenkins ͑1996͒. The theoretical analysis of Akkas ͑1978͒ shows the existence of instabilities during the dynamic inflation of spherical membranes. His theory, however, assumes that the membrane thickness remains constant during inflation.

In this paper, we first review the governing equations for the motion of a spherical membrane for a material obeying the Mooney-Rivlin constitutive equation, subjected to an arbitrary and time-dependent driving pressure. A summary of the various regimes in static inflation is then given. The dynamic inflation is extensively investigated to generate both the linear stability picture and the nonlinear behavior. The influence of various control parameters, namely, the material constant͑s͒ and driving pressure are examined for the initially unconstrained spherical membrane. We are particularly interested in examining the conditions for the onset of sustained oscillatory behavior.

II. EQUATION OF MOTION

In this section, the derivation of the governing equation of the motion of a spherical membrane under internal pressure is briefly reviewed for completeness. Although the equations for static inflation are well known ͓Beatty ͑1987͔͒, their dynamic counterparts are not as extensively used. The classical Mooney strain energy function is used to described the rubber-like material.

Consider a spherical membrane of elastic, isotropic, and incompressible material of density 0 . The membrane deformation is described in spherical coordinates (R,⌰,⌽) in the undeformed state, and (r,,) in the deformed state. Let (R 0 ,H 0 ) and (r 0 ,h 0 ) denote the membrane midsurface radius and thickness in the undeformed state and in the deformed state, respectively.

Let rr and ϭ be the diagonal components of the Cauchy stress tensors in the radial and circumferential directions, respectively. Because of the spherical symmetry, all variables depend only on the current radial coordinate r and time t. In this case, the governing equation of interest is the conservation of momentum equation in the radial direction:

0 r ¨ϭ ץ rr ץr ϩ 2 r ͑ rr Ϫ ͒. ͑1͒
Let 1 ϭ h 0 /H 0 denote the principal stretch ratio in the thickness direction, and 2 ϭ 3 ϭ ϭ r 0 /R 0 be the principal stretch ratios in the circumferential directions. In this case, the incompressibility assumption implies that 1 2 ϭ 1. Let ⌬ P denote the imposed pressure difference. Integration of Eq. ͑1͒ through the thickness, and use of the membrane approximation and the expressions for the stretch ratios lead to the following reduced form ͓Eringen and Suhubi ͑1974͔͒:

0 H 0 R 0 ¨ϭ ⌬P 2 ϩ 2H 0 R 0 1 ͑ rr Ϫ ͒. ͑2͒
The expression for the stress difference rr Ϫ will be determined through the constitutive relation between stress and strain ͑principal stretch ratio͒.

In the present work, the material studied is assumed to obey the Mooney-Rivlin constitutive equation ͓Mooney ͑1940͒; Treloar ͑1975͔͒. The corresponding expression for rr Ϫ is given in terms of the principal circumferential stretch ratio and the material constants C and ␣ as ͓Green and Adkins ͑1960͔͒

rr Ϫ ϭ Ϫ2C ͩ 2 Ϫ 1 4 ͪ ͑1ϩ␣ 2 ͒.

͑3͒

The differential equation ͑2͒ becomes

0 H 0 R 0 ¨ϭ ⌬P 2 ϩ 4CH 0 R 0 ͩ 1 5 Ϫ ͪ ͑1ϩ␣ 2 ͒.

͑4͒

To simplify the discussion on the influence of the different parameters in Eq. ͑4͒, we define the reduced time and pressure by

ϭ t 2 R 0 ͱ C 0 , ͑5͒ ⌬p ϭ ⌬P R 0 4CH 0 . ͑6͒
Equation ͑4͒ takes the following nondimensional form:

¨ϭ ⌬p 2 ϩ ͩ 1 5 Ϫ ͪ ͑1ϩ␣ 2 ͒,

͑7͒

in which is a function of the reduced time variable . Now a dot denotes differentiation with respect to .

III. DYNAMIC BEHAVIOR

The dynamic inflation is now examined by solving Eq. ͑7͒ for a constant inflation pressure step. However, it is helpful to consider first the static inflation of a spherical membrane, which will be referred to later when dynamic inflation is analyzed.

A. Static inflation

In Eq. ͑7͒, the acceleration term is set to zero, and the inflation pressure step becomes an explicit function of the material parameter ␣ and the principal circumferential stretch ratio :

⌬p ϭ ͩ 1 Ϫ 1 7 ͪ ͑1ϩ␣ 2 ͒.

͑8͒

The curves ⌬p vs are presented in Fig. 1 for three different values of ␣: 0., 0.1, and 0.25. The three curves are typically representative of three possible behaviors ͓Beatty ͑1987͔͒.

The curve ␣ ϭ 0 represents the behavior of a neo-Hookean membrane. The maximum pressure the membrane can withstand in static inflation is called the critical static pressure or limit pressure, and is denoted ⌬p S C . For 0 Ͻ ⌬p Ͻ ⌬p S C , Eq. ͑8͒ has two distinct real roots, one on the ascending portion of the curve ͑for Ͻ S C , where S C corresponds to ⌬p S C ͒ and another on the descending portion ( Ͼ S C ). For

⌬ p Ͼ ⌬p S C
, there is no real solution for Eq. ͑8͒: the membrane cannot withstand such an inflating pressure. For ⌬ p ϭ ⌬p S C , there is a double root, and for ⌬ p ϭ 0, there are two real distinct roots at 0 and ϩϱ. The curve ␣ ϭ 0.1 is representative of the situation in which the inflation curve always presents two stable and one unstable parts. This is the case in which 0 Ͻ ␣ Ͻ ␣ S lim , with ␣ S lim Ӎ 0.214. Here, the curve ⌬p vs consists of three branches. The curve has two local extrema, one maximum for ⌬p ϭ ⌬p S C ͑this is not a critical value, because the membrane can withstand greater pressures͒ at ϭ S C and one minimum for ⌬p ϭ ⌬p S min at ϭ S min . From Fig. 1, it can be seen that for 0 Ͻ ⌬p Ͻ ⌬p S min , Eq. ͑8͒ has only one real root in the first stable branch. For ⌬p S min Ͻ ⌬p Ͻ ⌬p S C , Eq. ͑8͒ has three real roots, one in the first stable branch, one in the unstable branch, and one in the second stable branch. For greater pressures (⌬p Ͼ ⌬p S C ), Eq. ͑8͒ presents only one real root in the second stable branch. It is to be noted that for ⌬p ϭ ⌬ p S min and ⌬p ϭ ⌬p S C there are only two real roots. In Fig. 1, the curve ␣ ϭ 0.25 stands for all the cases where ␣ Ͼ ␣ S lim . All the curves are monotonically increasing and the function ͑8͒ is a bijection: for every inflating pressure step ⌬p, there is only one real root.

B. Dynamic inflation

Equation ͑7͒ is rewritten as a set of two differential equations of the first order:

˙ϭ v, v ˙ϭ ⌬p 2 ϩ ͩ 1 5 Ϫ ͪ ͑1ϩ␣ 2 ͒.

͑9͒

The initial conditions are: where init and v init will be given when specific examples are covered. With these initial conditions, problem ͑9͒ and ͑10͒ is well posed and has a unique solution. Of course, given the nonlinearities present, the problem can admit a multiple-valued solution. However, the problem is deterministic, and only one solution corresponds to a set of initial conditions. The equations are solved by using a sixth-order Runge-Kutta method ͑subroutine IVPRK from the IMSL library͒.

͑ ϭ 0͒ ϭ init , ͑10͒ v͑ ϭ 0͒ ϭ v init ,
The dynamical system ͑9͒ is codimension 2 since two parameters are required for a bifurcation to occur. In this case, the two parameters are ⌬p and ␣. The coexistence of fixed points is illustrated in Fig. ͑2͒, which shows constant curves in the (⌬p,␣) plane. It is clear from Fig. 2 that a given stretch ratio may be maintained by an infinite number of combinations of the two parameters ͑each curve intersects an infinite number of other curves͒. More importantly, Fig. 2 shows that two different values of the stretch ratio can correspond to the same value of ⌬ p and ␣.

Since ⌬p is independent of time, v can be obtained as an explicit function of by taking the ratio of the two equations system ͑9͒ and integrating from ( init ,v init ) to (,v). We thus obtain the equation of energy conservation for the inflation of the spherical membrane ͓Dym ͑1974͔͒: In Eq. ͑11͒, the left-hand side is the change in kinetic energy, the first term on the right-hand side is the work done by the inflating pressure, and the last terms represent the elastic strain energy stored by the membrane. This equation allows us to directly obtain the phase curves v vs . We next examine the linear stability of the static ͑steady-state͒ solution.

1 2 v 2 Ϫ 1 2 v init 2 ϭ ⌬p ͩ 3 3 Ϫ init 3 3 ͪ Ϫ ͫ 2 2 ϩ 1 4 4 ϩ␣ ͩ 4 4 ϩ 1 2 2 ͪͬ ϩ ͫ init 2 2 ϩ 1 4 init 4 ϩ␣ ͩ init 4 4 ϩ 1 2 init 2 ͪͬ .

Linear stability analysis

Assuming that ␣ and ⌬p are fixed, the equilibrium states are the points that satisfy ˙ϭ 0 and v ˙ϭ 0 in system ͑9͒. These points are the couples ( e ,v e ) of the phase curves where e satisfies Eq. ͑8͒ and v e ϭ 0. The number of such points depends on the values of ␣ and ⌬p and was investigated in detail in Sec. III A.

In order to determine the stability of each equilibrium state, we must study the behavior of system ͑9͒ about this state. The Jacobian matrix ͓J͔ of system ͑9͒ at the equilibrium point ( e ,0) is given by

͓J͔ ϭ e ,v ϭ 0 ϭ ͫ 0 ⍀ 1 0 ͬ , ͑12͒ in which ⍀ ϭ 2⌬p e Ϫ ͩ 5 e 6 ϩ1 ͪ Ϫ␣ ͩ 3 e 4 ϩ3 e 2 ͪ .

͑13͒

Let 1 and 2 be the two eigenvalues of the problem, which dictate the approximate form of and v in the neighborhood of ( e ,0). The two roots 1 and 2 depend on the sign of ⍀:

͑i͒ If ⍀ Ͼ 0, the two roots are real and of opposite signs, 1 ϭ Ϫͱ⍀ and 2 ϭ ͱ⍀.

͑ii͒ If ⍀ Ͻ 0, the two roots are purely imaginary, 1 ϭ ϪiͱϪ⍀ and 2 ϭ iͱϪ⍀.

͑iii͒ If ⍀ ϭ 0, zero is the double root, 1 ϭ 2 ϭ 0.

Let us now examine the dependence of ⍀ ͑more importantly, the sign of ⍀͒ on ⌬p and ␣ ͓recall that e is determined from Eq. ͑8͔͒. If ⍀ Ͼ 0, the equilibrium point ͑ ϭ e , v ϭ 0͒ is unstable and is a saddle point; if ⍀ Ͻ 0, (t) and v(t) are periodic in the neighborhood of ͑ ϭ e , v ϭ 0͒, which is an elliptic center. The number and stability of the equilibrium branches depend on ⌬ p and ␣. If we note that

d⌬p d ͯ ϭ e ϭ Ϫ ⍀ e 2 , ͑14͒
then we can clearly see that if ⍀ Ͻ 0, the ascending branch͑es͒ of a given static curve in Fig. 1 is ͑are͒ stable. If ⍀ Ͼ 0, the descending branch of the static curve is unstable, which confirms previous stability findings ͓Alexander ͑1971͔͒. ͓Figure 1 also shows the span of stable and unstable branches with the limit points in ⌬p corresponding to the extrema in the static curves. Linear stability analysis is limited to small departures from equilibrium. Before analyzing the nonlinear dynamical behavior, it is helpful to examine the solution in the (,v) phase plane for the three values of the parameter ␣ as in Fig. 1. The resulting phase diagrams are shown in Fig. 3. For ␣ ϭ 0, if 0 Ͻ ⌬p Ͻ ⌬p S C , there are two equilibrium points: one elliptic center ( v ,0) and one saddle point ( s ,0) with v Ͻ s . The corre-sponding (,v) trajectories are shown in Fig. 3͑a͒. There are closed curves in the vicinity of the elliptic center, and these curves become open branches in the neighborhood of the saddle point. Depending on the initial conditions ( init ,v init ), the solution may oscillate around the elliptic center or grow exponentially away from the saddle point. Figure 3͑a͒ displays a global homoclinic ͑separatrix͒ solution confining the area of periodic solutions about the stable center and diverging solutions outside; note the similarity with Hamiltonian problems with ''escape from a potential well'' ͓Guckenheimer and Holmes ͑1983͔͒. The homoclinic solution is obtained from Eq. ͑11͒ using the initial conditions of the saddle. Figure 3͑b͒ displays the trajectories corresponding to ⌬p Ͼ ⌬p S C . Regardless of the initial conditions, the solution does not converge to a fixed point. This is, of course, resulting from the fact that there are no real solutions to Eq. ͑8͒. The phase plane trajectories in the presence of three fixed points are presented in Fig. 3͑c͒ for ␣ ϭ 0.1 and ⌬p S min Ͻ ⌬p Ͻ ⌬p S C . In this case, there are two elliptic centers ( v 1 ,0) and ( v 2 ,0) and one saddle point ( s ,0) with v 1 Ͻ s Ͻ v 2 . Depending on the values of the initial conditions, the trajectories can be two closed curves, each one of them around one elliptic center or only one closed curve that surrounds all three fixed points. Figure 3͑c͒ shows also the case of a periodic limit orbit that surrounds the three fixed points. In the former case, the solution tends to oscillate around one elliptic center or the other. In the latter case, the solution remains periodic but oscillates around both elliptic centers over one period of time. The trajectory turns around one of the elliptic centers, then enters the saddle-point neighborhood and goes on around the other elliptic center. For ⌬ p ͓0,⌬p S min ͓U͔⌬p S C ,ϩϱ͓, there is only one fixed point, which is an elliptic center. Thus, the trajectories ͑not shown͒ are closed curves centered on this point, regardless of the initial conditions. Figure 3͑c͒ thus exhibits two homoclinic solutions confining ͑internal͒ periodic solutions about the centers and a set of ͑external͒ periodic solutions about both homoclinic solutions. This situation is similar to asymmetric double-well problems. Note that for ␣ ϭ 0.25 and for any value of ⌬p, there is only one fixed elliptic center.

Nonlinear dynamic inflation

In the previous section, we examined the phase-plane trajectories for general initial conditions. We determined the number of fixed points and analyzed their stability as function of the control parameters ␣ and ⌬p. Consider now the special case of a membrane initially in equilibrium subject to the following initial conditions ͑ init ϭ 1, v init ϭ 0͒. This corresponds to a spherical membrane of a given radius initially in equilibrium.

We first examine the conditions for a periodic response. Obviously, when the motion of the membrane is periodic, the corresponding phase curve is a closed loop in the (,v) plane. Noting that the initial point ͑1,0͒ is an intersection point between the v ϭ 0 axis and the phase trajectory, it can be inferred that this trajectory must intersect the v ϭ 0 axis at least one more time for it to be a closed loop. If v is set equal to zero in Eq. ͑11͒, the relation between ⌬p and becomes

⌬p ϭ ͭ ͩ 3 3 Ϫ1 ͪͫ 2 2 ϩ 1 4 4 ϩ␣ ͩ 4 4 ϩ 1 2 2 ͪ Ϫ 3 4 ͑1ϩ␣͒ ͬ ͑ 1͒ 0 ͑ ϭ 1͒ . ͑15͒
The curves (,⌬p), corresponding to Eq. ͑15͒, are presented in Fig. 4, together with the static curves based on Eq. ͑8͒, for different values of ␣. Thus, the solid curves in Fig. 4 correspond to the loci of solutions with zero velocity (v ϭ 0) and nonzero acceleration (v ˙ 0), whereas the dashed ͑static͒ curves correspond to the loci of solutions with both v ϭ 0 and v ˙ϭ 0. It is clear from Fig. 4 that the overall behavior of the (,⌬p) curves is the same for both static and dynamic inflations. There is a value of the material parameter, ␣ D lim Ӎ 0.145, beyond which, similarly to ␣ S lim , ⌬p increases monotonically with . The ␣ ϭ 0 curve has one maximum corresponding to a critical pressure value, It is useful to closely examine the interplay between the static and dynamic behaviors through the corresponding curves, since the static curves represent equilibrium regions towards which the dynamic solution may eventually be attracted in real systems ͑with damping͒. There are four important ranges of values for ␣. Recall that in static inflation three main ranges of interest were examined: ␣ ϭ 0, 0 Ͻ ␣ Ͻ ␣ S lim , and ␣ Ͼ ␣ S lim . In the dynamic case, we identify four regions of interest, namely, ␣ ϭ 0, 0

⌬p D C . For 0 Ͻ ␣ Ͻ ␣ D lim ,
Ͻ ␣ Ͻ ␣ D lim , ␣ D lim Ͻ ␣ Ͻ ␣ S lim
, and ␣ Ͼ ␣ S lim . We are particularly interested in the conditions under which oscillatory movement occurs; we must then determine the location of the first intersection point i ( i Ͼ 1), between a dynamic curve and a horizontal line, representing the value of a given pressure step ⌬p. Note that the same horizontal line intersects the static curve ͑corresponding to the same ␣ value͒ at e ͑there may be more than one intersection point͒. The solution may thus end up oscillating between ϭ 1 and i around e .

We now focus our discussion on the range ␣ Ͼ 0.1. Figures 5 and6 display the resulting dynamics for two different values of ␣. ␣ ϭ 0.18 and 0.25, respectively, with initial conditions init ϭ 1, v init ϭ 0. Each figure shows ͑a͒ the dynamic and static pressure curves together with the intersecting points i 's and e 's for a typical level of inflating pressure, ͑b͒ the phase-plane trajectories for various levels of inflating pressure, and ͑c͒ the corresponding temporal behaviors for . Typically, there are four different ⌬p ranges that can be identified. For 0 Ͻ ⌬p Ͻ ⌬p D C , there are generally two intersection points, i 1 and i 2 ( i 1 Ͻ i 2 ), between the dynamic curve and the horizontal line, and two intersecting ͑equilibrium͒ points, e 1 and e 2 ( e 1 Ͻ e 2 ), for the static curve. In this case, the membrane oscillates between the unconstrained ͑initial͒ state and i 1 around the stable static equilibrium ͑elliptic͒ point e 1 . For the particular case ⌬p ϭ ⌬p D C , there is only one intersection point i and still two intersection points e 1 and e 2 . This is, in fact, a degenerate case where the two roots of Eq. ͑15͒ are equal to i . Note also that i is equal to e 2 , this latter being a saddle point. In this case, the membrane neither oscillates nor inflates any further but maintains a static ͑unstable͒ equilibrium state. This corresponds to an oscillation with an infinite period. For ure 5 shows the behavior for only three ranges of inflating pressure, since the dynamic curve exhibits a strictly monotonic behavior ͓Fig. 5͑a͔͒. Note that in this case the static curve still exhibits two extrema. The phase-plane trajectory is always a closed cycle but its shape depends on the value of ⌬p relative to the extrema of the static curve. For FIG. 5. ␣ ϭ 0.18: ͑a͒ ⌬p vs for static and dynamic inflations, ͑d͒ stable phase-plane trajectories corresponding to init ϭ 1 and v init ϭ 0, and ͑c͒ temporal behavior of the circumferential stretch ratio , for four ⌬ p ranges:

͑-͒ 0 Ͻ ⌬p Ͻ ⌬ p D C , ͑--͒ ⌬p ϭ ⌬ p D C , ͑¯͒ ⌬p D C Ͻ ⌬ p Ͻ ⌬p S C , and ͑-••-͒ ⌬p у ⌬p S C .
⌬p Ͻ ⌬p S min , there is only one equilibrium point on the first stable branch of the static curve and the membrane oscillates around it. For ⌬p S min Ͻ ⌬p Ͻ ⌬p S C , there are three equilibrium points, and the membrane oscillates around the two vortices as shown in Figs. 5͑b͒ and5͑c͒. Note that the necking is insignificant in this case since the motion amplitude is relatively small. For ␣ ϭ 0.25, both the dynamic and static curves are strictly monotonic. In this case, condition ͑15͒ has only one root and there is always one equilibrium point. The membrane always oscillates. Therefore, the period of motion remains finite. Figure 6 shows the corresponding behavior for three different values of the pressure step. After having determined the ranges of values for ␣ and ⌬p for which oscillatory inflation occurs, we next examine the influence of these two parameters on the period of oscillation. It was found above that the membrane tends to oscillate in most cases of practical interest, and that the amplitude of the motion depends on the value of the pressure step ͑for given ␣͒. Similarly, Figs. 5 and6͑c͒ show that the period of oscillation is closely related to ⌬ p. Figure 7 shows the dimensionless period T as function of ⌬p for the four values of ␣ considered above. Generally, the period of oscillation tends to increase at small inflating pressure, reaches a maximum value and then decreases in the higher ⌬p range. The only case in which the membrane becomes unstable is when ␣ ϭ 0 ͑neo-Hookean material͒ and ⌬p Ͼ ⌬p D C . For ␣ ϭ 0 and ⌬p ϭ ⌬p D C , and ␣ ϭ 0.1

͑representing the case 0 Ͻ ␣ Ͻ ␣ S lim
͒ and ⌬p ϭ ⌬ p D C , the membrane tends to reach an unstable equilibrium ͑saddle͒ point and the period goes to infinity. After this point, for ␣ ϭ 0.1 and ⌬p Ͼ ⌬p D C , the membrane oscillates again and the period of motion decreases monotonically. In all the other situations, the membrane oscillates. But, because of the nonlinearity of the differential system, it is not enough to conclude that the period of motion grows as the pressure step grows ͓see Figs. 5 and 6͑c͔͒: for ␣ ϭ 0.25 ͑representing ␣ Ͼ ␣ S lim ͒, the motion is never unstable but the period does not increase monotonically with pressure. The corresponding curve T vs ⌬p admits a maximum. For ␣ ϭ 0.18, the results are similar to the latter case, because ␣ is also greater than ␣ S lim .

This last result differs from Akkas ͑1978͒ who concluded that, for ␣ Ͼ ␣ S lim , as ⌬p gets larger, the period of oscillations also gets larger. This difference may be due to his hypothesis of constant thickness of the membrane during inflation.

IV. DISCUSSION AND CONCLUSION

There are two major issues regarding the class of materials targeted in this study that we would like to discuss in some detail. The first is the use of the Mooney-Rivlin model as opposed to using other elastic models. The second is the fact that the material in question does not include any viscous damping. Thus, it is important to point out the limitations and advantages of the assumptions adopted here. Hyperelastic models, including the Mooney-Rivlin model, can capture some of the phenomena observed during the ͑primary͒ inflation stage in a wide range of forming processes. Beatty ͑1987͒ compared the normalized pressure-stretch curve obtained from the experimental inflation of a spherical balloon against results based on the Mooney-Rivlin model. He found that, while the normalized pressure grows indefinitely with stretch, the slope of the Mooney-Rivlin curve approaches a constant value for large stretch. Hence, the Mooney-Rivlin model does not predict an indefinitely large, bursting pressure at a finite value of the stretch. This shortcoming of the Mooney-Rivlin model is by now well accepted, as the model can only be accurate for moderately large and not very large deformation. Nonetheless, according to Beatty ͑1987͒, ''the Mooney-Rivlin model does provide a satisfactory qualitative ͑and probably conservative͒ description of the overall primary balloon inflation phenomenon.'' Although the discrepancy between theory and experiment could be in part due to the fact that the balloon was not perfectly spherical during the experiment, and also to the assumption of incompressibility of the material applied throughout the analysis. Obviously, some additional and systematical work is needed to examine the effects of materials properties on inflation dynamics.

On the other hand, the present study should not be regarded as limited to the stability analysis of only Mooney-Rivlin materials, but as a study relevant to a large class of hyperelastic materials in the limit of moderately large strain. The methodology developed in the present study can be applied to other constitutive models, and the type and range of phenomenological behaviors examined are likely to be encountered in other hyperelastic materials. The Mooney-Rivlin model is related to or is the limit of many more realistic models for hyperelastic materials. These include the following models: generalized Rivlin ͓Rivlin ͑1948͔͒, three-term model of Isihara et al. ͑1951͒, Rivlin-Saunders ͑1951͒, Gent-Thomas ͑1958͒, Biderman ͑1958͒, Carmichael-Holdaway ͑1961͒, Hart-Smith ͑1966͒, Valanis-Landel ͑1967͒, Alexander ͑1968͒, Ogden ͑1972͒, and Blatz-Sharda-Tschoegl ͑1974͒. These models are derived on the basis of phenomenological, statistical, or empirical arguments. Our own studies ͓Khayat and co-workers ͑1998͒; ͑1995͒; ͑1994a͒; ͑1994b͒; ͑1993͒; ͑1992͔͒ on the modeling of blow-molding and thermoforming processes showed a remarkable usefulness for the Mooney-Rivlin model, for materials undergoing industrially complex transformations. See, in particular, Khayat andDerdouri ͑1995͒ and[START_REF] Derdouri | Study of the Simultaneous Extension and Inflation of PET Performs[END_REF] for comparison between theory and experiment. Similarly, other works also support this claim ͓see, for instance, Charrier et al. ͑1989͒ and Delorenzi and Nied ͑1990͔͒.

Let us now turn to the second issue regarding the absence of damping in the model. In this regard, one may conclude hastily that the purely elastic model in this case is of no relevance to physical reality. It is true that damping will make the oscillatory solution decay to a fixed point. However, there are instances where transient and not only longterm behavior is of the main physical or practical interest. It is possible that the inflation process in question may be so short that sustained oscillations, which stem from dominant elastic effects, are observed during the whole process. In blow molding and thermoforming, for instance, the inflation stage may last only a few milliseconds, and viscous effects are negligible. Moreover, since the sustenance of oscillations depends on the effects of elasticity relative to viscous effects, it is the type of ͑polymer͒ material used that dictates the time it takes for oscillations to decay ͑relaxation time͒. In a previous study ͓Khayat and Garcia-Rejon ͑1992͔͒ we used a highly nonlinear viscoelastic constitutive model to study the inflation of a spherical shell of a viscoelastic fluid. It was found that beyond a critical Deborah number, oscillatory ͑as opposed to monotonic͒ growth sets in. It was also found that the amplitude and frequency of oscillations depend strongly on the pressure history inside the shell; a sudden initial jump in pressure tends to enhance oscillatory behavior.

In conclusion, the dynamic inflation of a spherical membrane made of Mooney-Rivlin material and subjected to a uniform step pressure was extensively studied. The nature of singular points was determined by a linear stability analysis for different ranges of values of the material parameter and the pressure step. The problem of the dynamic inflation of an initially unconstrained membrane exhibits critical values of the material parameter ␣ D lim and of the pressure step, ⌬p D C ͑counterparts of the classical static critical values ␣ S lim and ⌬p S C ͒ ͓Alexander ͑1971͒; Beatty ͑1987͔͒. We examined, in particular, the conditions for the onset of oscillatory inflation around the ͑static͒ equilibrium points. It was found that oscillatory motion tends to set in for Mooney-Rivlin materials with a period that is maximum at some inflating pressure. It was also found that the period tends to decrease as the Mooney parameter ratio increases.

Although the equations governing the inflation of the spherical membrane are Hamiltonian ͑as they do not account for any damping͒, it is difficult to assess their physical or practical relevance. A more realistic model should include the effects of viscosity, leading to a viscoelastic solid or fluid model. The inflation of a viscoelastic liquid shell obeying the modified ZFD model was earlier examined by Khayat and Garcia-Rejon ͑1992͒. It is found that, for a pressure suddenly imposed initially, the shell grows exponentially with time as long as the Deborah number De is small. When De exceeds a critical level, the growth becomes oscillatory. However, the amplitude and frequency of oscillations is strongly influenced by the rate of internal pressure. Calculations show that the oscillations are essentially absent for a slowly increasing pressure ͑compared to an initial sudden step͒. Since the pressure, in practice, cannot be suddenly imposed, it is likely that oscillations do not emerge in reality as a result of the strong damping.

FIG. 1 .

 1 FIG. 1. Static inflation: normalized pressure step vs circumferential stretch ratio for different values of ␣ shown on each curve.

͑11͒FIG. 2 .

 2 FIG. 2. Codimension-two bifurcation: normalized pressure step vs ␣ for different values of the circumferential stretch ratio .

FIG. 3 .

 3 FIG. 3.Phase-plane trajectories v vs for three typical ranges of ␣ and ⌬ p and general initial conditions ( init ,v init ).

  there are three branches, a local maximum ⌬p D C , and a local minimum ⌬p D min , as in the static case. Note that the values of the critical parameters are highly dependent on the initial conditions. It is interesting to note that ␣ D lim Ͻ ␣ S lim and ⌬ p S min Ͻ ⌬p D min Ͻ ⌬p D C Ͻ ⌬p S C independently of the value of ␣.

FIG. 4 .

 4 FIG. 4. ⌬p vs for static ͑¯͒ and dynamic ͑--͒ inflations ͑the corresponding value of ␣ is shown on each curve͒.

  ⌬p D C Ͻ ⌬ p Ͻ ⌬p S C , the phase trajectory is reduced to an open curve. The membrane tends to inflate indefinitely, but because of the presence of the saddle point the phaseplane trajectory exhibits a local minimum in the neighborhood of the saddle point. For ⌬ p у ⌬p S C , there are no intersection points at all and the phase-plane trajectory increases monotonically with . Consider now the case when ␣ D lim Ͻ ␣ Ͻ ␣ S lim , more precisely, when ␣ ϭ 0.18. Fig-

FIG. 6 .

 6 FIG. 6. ␣ ϭ 0.25: ͑a͒ ⌬p vs for static and dynamic inflations, ͑b͒ stable phase-plane trajectories corresponding to init ϭ 1 and v init ϭ 0, and ͑c͒ temporal behavior of the circumferential stretch ratio , for three ⌬ p ranges: ͑-͒ 0 Ͻ ⌬p Ͻ ⌬ p S min , ͑--͒ ⌬p S min Ͻ ⌬p Ͻ ⌬p S C , and ͑¯͒ ⌬p у ⌬p S C .

FIG. 7 .

 7 FIG. 7.Influence of the pressure step on the period of motion for different values of ␣.