Razvan Barbulescu

Jérémie Detrey

Nicolas Estibals

Paul Zimmermann

Finding Optimal Formulae for Bilinear Maps

Keywords: optimal algorithms, polynomial multiplication and squaring, finite field arithmetic, tensor rank, bilinear map, bilinear rank

We describe a unified framework to search for optimal formulae evaluating bilinear -or quadratic -maps. This framework applies to polynomial multiplication and squaring, finite field arithmetic, matrix multiplication, etc. We then propose a new algorithm to solve problems in this unified framework. With an implementation of this algorithm, we prove the optimality of various published upper bounds, and find improved upper bounds.

Introduction

This article studies the bilinear rank problem. Given a field K, this problem naturally arises when considering the computational complexity of formulae (or algorithms) for evaluating bilinear maps over K [START_REF] Bürgisser | Algebraic complexity theory[END_REF]Ch. 14]. For instance, typical bilinear maps include -but are not limited to -the multiplication of two n-term polynomials of K [X], or of two r × r square matrices of K r×r .

Note that, given an algorithm for computing a bilinear map over K, this algorithm can be naturally extended to compute the same bilinear map over any larger K-algebra K. For instance, given a formula for computing the product of two n-term polynomials over K = F 2 , the exact same formula can also be used to compute the product of two n-term polynomials over any field extension K = F 2 m , or even the product of two n-term polynomials of F r×r 2 [X], i.e., polynomials having binary r × r matrices as coefficients. A nice illustration of the latter is given by Albrecht in [START_REF] Albrecht | The M4RIE library for dense linear algebra over small fields with even characteristic[END_REF], where the multiplication of large r × r matrices over F 2 n is implemented by representing F r×r 2 n as F r×r 2 [X]/ f (X) • F r×r 2 [X] , with f an irreducible polynomial over F 2 of degree n, and then using Montgomery's n-term polynomial multiplication formulae [START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF].

Therefore, when dealing with such formulae, a crucial distinction is made between a full multiplication a•b of two elements a and b derived from the inputs of the bilinear map, and thus possibly living in the larger algebra K; and an addition a + b or a scalar multiplication λa, where λ is given by the bilinear map only, and thus belongs to the smaller coefficient field K.

Intuitively, the first one is expected to have a much higher computational cost than the last two. Table 3 in [START_REF] Albrecht | The M4RIE library for dense linear algebra over small fields with even characteristic[END_REF] shows that the time needed to multiply two matrices over fields of characteristic 2 is proportional to the number of full multiplications in the formulae used for this purpose. A sensible way to optimize the overall complexity of the formula at hand is therefore to minimize the number of these full multiplications, which corresponds to solving a particular instance of the bilinear rank problem. Most of the fast algorithms for polynomial or integer multiplication can be expressed in terms of the bilinear rank problem, for example Karatsuba's algorithm [START_REF] Karatsuba | Multiplication of multi-digit numbers on automata (in Russian)[END_REF] and Toom's algorithm [START_REF] Toom | The complexity of a scheme of functional elements realizing the multiplication of integers[END_REF]. But this applies to other problems as well, like the middle product [START_REF] Hanrot | The middle product algorithm, I. Speeding up the division and square root of power series[END_REF], or matrix multiplication [START_REF] Bläser | On the complexity of the multiplication of matrices of small formats[END_REF][START_REF] Courtois | A new general-purpose method to multiply 3 × 3 matrices using only 23 multiplications[END_REF][START_REF] Strassen | Gaussian elimination is not optimal[END_REF].

Definition of the problem. Let K be a field. Given an n × m bilinear map φ : K n × K m → K , the bilinear rank problem consists in finding formulae for evaluating φ involving a minimal number k of full multiplications. This optimal k is called the bilinear rank of the map φ.

More formally, consider B(K n , K m ; K) the set of n×m bilinear forms from K n ×K m to K. Each such bilinear form γ can be written as γ : a, b → i,j γ i,j a i b j , with the two input vectors a = (a 0 , a 1 , . . . , a n-1 ) ∈ K n , b = (b 0 , b 1 , . . . , b m-1 ) ∈ K m , and the coefficients γ i,j ∈ K. As γ is uniquely determined by the γ i,j 's, we can represent B(K n , K m ; K) as the vector space of dimension nm over K, which we will denote by V in the following. The representation of the bilinear form γ in V is therefore the nm-dimensional vector γ = (γ 0,0 , γ 0,1 , . . . , γ n-1,m-1 ), and a bilinear map φ from K n × K m to K can then be represented as a tuple of vectors of V : φ = (φ 0 , φ 1 , . . . , φ -1 ) ∈ V .

Consider now an n × m bilinear form γ. We say that γ is of rank 1 if there exist two linear forms α :

K n → K, a → i α i a i and β : K m → K, b → j β j b j such that γ(a, b) = α(a) • β(b) for all a ∈ K n and b ∈ K m .
Rank-1 bilinear forms are particularly relevant in the context of the bilinear rank problem: indeed, each such form can be evaluated using only one full multiplication, as evaluating the two linear forms α and β at a and b, respectively, only requires scalar multiplications and additions. Furthermore, since the maps e i,j : a, b → a i b j are all rank-1 bilinear forms, any bilinear form can then be written as a linear combination of rank-1 bilinear forms: γ = i,j γ i,j e i,j . Therefore, one can ask the question: What is the minimal number of rank-1 bilinear forms necessary to evaluate a given bilinear form? A given bilinear map?

In the case of a (single) bilinear form γ, the answer is easy: the bilinear rank is given by the actual rank of the n × m matrix (γ i,j ) 0≤i<n,0≤j<m formed by the coefficients of γ. However, the minimal number of full multiplications is much more difficult to compute when evaluating ≥ 2 bilinear forms simultaneously, such as is the case with bilinear maps: this is the aforementioned bilinear rank problem, formalized below.

Definition 1 (Bilinear rank problem [START_REF] Bürgisser | Algebraic complexity theory[END_REF]Ch. 14]). Using the above notations, and given a finite generator set G ⊂ V composed of rank-1 n × m bilinear forms, along with a finite target set of bilinear forms T = {t 0 , t 1 , . . . , t -1 } ⊆ Span(G), the bilinear rank problem consists in generating all the elements of T by K-linear combinations of a minimal number k of elements of G or, alternatively, to find all solutions for that optimal k.

Without loss of generality, the target vectors t i can further be assumed to be linearly independent, as reconstructing an extra target vector by linear combination would only require scalar multiplications and additions.

Note that this problem can be adapted to also encompass the case of sets of quadratic forms. Indeed, any quadratic form σ from K n to K, σ : a → 0≤i≤j<n σ i,j a i a j , is uniquely represented by its coefficients σ i,j ∈ K. The set of n-ary quadratic forms Q(K n ; K) can thus be seen as a vector space of dimension n(n+1)/2 over K. It suffices then to take V to be this vector space, and to define the rank-1 quadratic forms to be the quadratic forms σ for which there exist two linear forms over K n , α : a → i α i a i and α :

a → i α i a i , such that σ(a) = α(a) • α (a) for all a ∈ K n .
In fact, any n × m bilinear form can be seen as an (n + m)-ary quadratic form σ such that no product a i a j occurs between the first n or between the last m input variables: σ i,j = 0 for all 0 ≤ i ≤ j < n and for all n ≤ i ≤ j < n + m. Hence, computing the rank of a set of quadratic forms is more general than computing the bilinear rank. Also, the bilinear rank problem is NP-hard. Indeed, Bürgisser et al. show in [START_REF] Bürgisser | Algebraic complexity theory[END_REF]Sec. 14.2] its equivalence to the decomposition of an order-3 tensor, known to be NPhard [START_REF] Hastad | Tensor rank is NP-complete[END_REF].

Notation. Note that in the rest of this document, by an abuse of notation, we will omit writing the input vectors of the considered bilinear or quadratic forms, as they will always be a, b or a, respectively. Hence, we will simply write γ = i,j γ i,j a i b j or σ = i≤j σ i,j a i a j to refer to the corresponding forms in the following. This amounts to implicitly considering the a i 's and b j 's as formal variables over K.

Related results. Several authors have considered special instances of the bilinear rank problem. We do not claim an exhaustive state-of-the-art here, but we mention the main results related to our work.

For polynomial multiplication, evaluation-interpolation algorithms like Karatsuba's and Toom's algorithms [START_REF] Karatsuba | Multiplication of multi-digit numbers on automata (in Russian)[END_REF][START_REF] Toom | The complexity of a scheme of functional elements realizing the multiplication of integers[END_REF] are special cases of the problem, where the only full multiplications involved correspond to evaluations of the product of the two polynomials at different points.

Even though the bilinear rank problem was already well known in the algebraic complexity community [START_REF] Bürgisser | Algebraic complexity theory[END_REF], until only recently, all the formulae used in the computer arithmetic community were based on this evaluation-interpolation scheme. In [START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF], instead of using only the evaluation products, i.e., ( i κ i a i )( j κ j b j ) for some κ ∈ K, Montgomery considered other rank-1 bilinear forms. This allowed him to find formulae with less products for the 5-, 6-, and 7-term polynomial multiplications, which do not fit in any evaluation-interpolation scheme. Montgomery however restricted the exploration to the set of "symmetric" generators of the form ( i α i a i )( j α j b j ). Our work shows how to improve Montgomery's exploration and proves that the formulae in [START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF] are optimal for the 5-term polynomial multiplication.

In [START_REF] Chung | Asymmetric squaring formulae[END_REF], Chung and Hasan propose asymmetric squaring formulae over Q, found using an exhaustive search method similar to that of Montgomery, but starting from an adhoc set of rank-1 generators. In [START_REF] Fan | Comments on five, six, and seven-term Karatsuba-like formulae[END_REF], Fan and Hasan use Montgomery's formulae and the Chinese Remainder Theorem together with a short product for the high degree terms to improve the bounds over F 2 ; those results were further improved by Cenk and Özbudak [5]. In [START_REF] Oseledets | Optimal Karatsuba-like formulae for certain bilinear forms in GF(2)[END_REF], using both exhaustive and heuristic search algorithms, Oseledets considers the multiplication of n-term polynomials, and also their product modulo X n -which corresponds to Mulders' short product [START_REF] Mulders | On short multiplications and divisions[END_REF] -and modulo X n + 1. Some of the linear algebra routines he proposed are a building block in our method.

Multiplication in finite field extensions F q n for n ≥ 2 reduces to polynomial multiplication modulo a degree-n irreducible polynomial f over F q . Thus, one could first compute a full product of two n-term polynomials, then reduce it modulo f using only scalar multiplications; however it is sometimes faster to directly compute the n terms of the product modulo f [START_REF] Cenk | On multiplication in finite fields[END_REF].

Fast matrix multiplication is another application of the bilinear rank problem, with the smallest unsolved problem being the multiplication of two 3 × 3 matrices (for two square matrices): we know that at least 19 products are needed, and that 23 are enough, but the minimal rank is still unknown [START_REF] Bläser | On the complexity of the multiplication of matrices of small formats[END_REF][START_REF] Courtois | A new general-purpose method to multiply 3 × 3 matrices using only 23 multiplications[END_REF].

Contributions. The contributions of the article are the following. First, we present the bilinear rank problem, as it is known to the algebraic complexity community [START_REF] Bürgisser | Algebraic complexity theory[END_REF], showing why it is pertinent to the computer algebra and computer arithmetic communities. Following Montgomery [START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF], we see this problem in the light of linear algebra, which enables one to search for optimal formulae for a large set of applications such as polynomial multiplication or squaring, matrix multiplication, etc. Second, we propose a new algorithm to solve the bilinear rank problem; this algorithm is faster than Montgomery's exploration algorithm [START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF] and is an improvement on Oseledets' heuristic Minbas [START_REF] Oseledets | Optimal Karatsuba-like formulae for certain bilinear forms in GF(2)[END_REF], which might not find the minimal rank. Last, using this new algorithm we prove the optimality of several known formulae, and exhibit new bounds for some bilinear maps.

Roadmap. The article is organized as follows. After detailing a few instances of the bilinear rank problem in Section 2, we show in Section 3 how this problem can be translated into a linear algebra problem. Section 4 gives an efficient algorithm solving this linear algebra problem. Finally Section 5 presents experimental results proving that some bounds from the literature are optimal, along with an improved bound for the multiplication in

F 3 5 .

Some Instances of the Bilinear Rank Problem

Throughout the rest of this paper, we will use the following running example:

Example 1 (2 × 3-term polynomial product). We want to multiply A = a 1 X + a 0 by B = b 2 X 2 + b 1 X + b 0 in K[X]
. The naive algorithm requires 6 products, while only 5 products are necessary when using Karatsuba's trick: To illustrate the generality of the proposed framework, we give a few more examples:

A • B = g 3 X 3 + (g 1 + g 2 )X 2 + (g 4 -g 2 -g 0 )X + g 0 , where g 0 = a 0 b 0 , g 1 = a 0 b 2 , g 2 = a 1 b 1 , g 3 =
Example 2 (3-term polynomial squaring [START_REF] Chung | Asymmetric squaring formulae[END_REF]). We want to square the polynomial A = a 2 X 2 + a 1 X + a 0 ∈ K [X]. For this quadratic map, we have n = 3, the generator set G consists of the 28 products

a 2 0 , a 0 a 1 , a 0 a 2 , a 0 (a 0 + a 1 ), a 0 (a 0 + a 2 ), a 0 (a 1 + a 2 ), a 0 (a 0 + a 1 + a 2 ), a 2 1 , a 1 a 2 , a 1 (a 0 + a 1 ), a 1 (a 0 + a 2 ), a 1 (a 1 + a 2 ), a 1 (a 0 + a 1 + a 2 ), a 2 2 , a 2 (a 0 + a 1 ), a 2 (a 0 + a 2 ), a 2 (a 1 + a 2 ), a 2 (a 0 + a 1 + a 2 ), (a 0 + a 1 ) 2 , (a 0 + a 1 )(a 0 + a 2 ), (a 0 + a 1 )(a 1 + a 2 ), (a 0 + a 1 )(a 0 + a 1 + a 2 ), (a 0 + a 2 ) 2 , (a 0 + a 2 )(a 1 + a 2 ), (a 0 + a 2 )(a 0 + a 1 + a 2 ), (a 1 + a 2 ) 2 , (a 1 + a 2 )(a 0 + a 1 + a 2 ), (a 0 + a 1 + a 2 ) 2 ,
and the target set is

T = {t 0 := a 2 0 , t 1 := 2a 0 a 1 , t 2 := a 2 1 + 2a 0 a 2 , t 3 := 2a 1 a 2 , t 4 := a 2 2 },
which are the 5 quadratic forms corresponding to the coefficients of A 2 .

Example 3 (Middle product [START_REF] Hanrot | The middle product algorithm, I. Speeding up the division and square root of power series[END_REF]). Assume we only want to compute the degree-1 and -2 coefficients of the product A • B from Example 1. We have n = 2, m = 3, and G is the set of all rank-1 bilinear forms of the form (α

0 a 0 + α 1 a 1 )(β 0 b 0 + β 1 b 1 + β 2 b 2 ) for α i , β j ∈ K. The target set is T = {t 0 := a 0 b 1 + a 1 b 0 , t 1 := a 0 b 2 + a 1 b 1 }.
A solution with k = 3 considers the subset

{g 0 := a 0 (b 1 + b 2 ), g 1 := a 1 (b 0 + b 1 ), g 2 := (a 1 -a 0 )b 1 } ⊆ G,
which gives the formulae t 0 = g 1 -g 2 and t 1 = g 0 + g 2 .

Example 4 (3 × 3 matrix multiplication). Here, n = m = 9, and we consider the product A • B over K 3×3 of the two matrices

A =   a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8   and B =   b 0 b 1 b 2 b 3 b 4 b 5 b 6 b 7 b 8   .
The target set T consists of the 9 bilinear forms {a

0 b 0 + a 1 b 3 + a 2 b 6 , ..., a 6 b 2 + a 7 b 5 + a 8 b 8 },
and the generator set G consists of the non-zero rank-1 bilinear forms (α

0 a 0 + • • • + α 8 a 8 )(β 0 b 0 + • • • + β 8 b 8 ) for α i , β j ∈ K, i.e., (2 9 -1) 2 = 261 121 forms for K = F 2 .
3 From Bilinear Applications to Linear Algebra

We focus in the rest of the paper on bilinear maps, the case of quadratic maps being similar. Recall that we denote by V the vector space of dimension nm over K isomorphic to the space of n × m bilinear forms B(K n , K m ; K). Thus, each element of T and G (see Definition 1) being such a bilinear form, it can be written as a row vector of dimension nm, where the (im + j)-th column corresponds to the coefficient of the a i b j term.

Example 1 (Cont'd). For our running example of the 2×3-term polynomial multiplication over K = F 2 , the canonical basis of the vector space V is (a

0 b 0 , a 0 b 1 , a 0 b 2 , a 1 b 0 , a 1 b 1 , a 1 b 2 ). The target set is T = {a 0 b 0 , a 0 b 1 + a 1 b 0 , a 0 b 2 + a 1 b 1 , a 1 b 2 }
, and the set of generators G consists of the 21 products:

G = { a 0 b 0 , a 1 b 0 , (a 0 + a 1 )b 0 , a 0 b 1 , a 1 b 1 , (a 0 + a 1 )b 1 , a 0 (b 0 + b 1 ), a 1 (b 0 + b 1 ), (a 0 + a 1 )(b 0 + b 1 ), a 0 b 2 , a 1 b 2 , (a 0 + a 1 )b 2 , a 0 (b 0 + b 2 ), a 1 (b 0 + b 2 ), (a 0 + a 1 )(b 0 + b 2 ), a 0 (b 1 + b 2 ), a 1 (b 1 + b 2 ), (a 0 + a 1 )(b 1 + b 2 ), a 0 (b 0 + b 1 + b 2 ), a 1 (b 0 + b 1 + b 2 ), (a 0 + a 1 )(b 0 + b 1 + b 2 ) }.
We can then rewrite T as the following matrix of 4 row vectors:

T =     a 0 b 0 a 0 b 1 a 0 b 2 a 1 b 0 a 1 b 1 a 1 b 2 1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1     a 0 b 0 a 0 b 1 + a 1 b 0 a 0 b 2 + a 1 b 1 a 1 b 2
The same applies to the set of generators G, which gives the 21 × 6 matrix:

G =        a 0 b 0 a 0 b 1 a 0 b 2 a 1 b 0 a 1 b 1 a 1 b 2 1 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 . . . . . . . . . . . . . . . . . . 1 1 1 1 1 1        a 0 b 0 a 1 b 0 (a 0 + a 1 )b 0 . . . (a 0 + a 1 )(b 0 + b 1 + b 2 )
The bilinear rank problem is then stated as follows: For a given integer k, we want to find the subsets W of k elements of G such that the subspace W spanned by W contains T .

In linear algebra terms, we want to find the subsets of k rows of G whose K-linear span contains all the row vectors of T , i.e., contains the subspace T = Span(T ).

Solving the Linear Algebra Problem

We recall the notations from Definition 1: V is the vector space of dimension nm spanned by the a i b j products, G ⊂ V is a finite set of generators, T ⊆ Span(G) is the set of target vectors, and T = Span(T ) is the corresponding target space spanned by the target vectors, which has dimension .

In the previous section, we reduced our problem of finding all formulae with k full multiplications for evaluating a given bilinear map to a linear algebra problem: given a finite subset G of a K-vector space V and a target subspace T = Span(T ) of Span(G), we want to find all the rank-k subspaces W of V containing T and which can be generated by elements of G only (i.e., Span(W ∩ G) = W ). One should note that this linear algebra problem is more general than the original problem of optimizing the computation of bilinear maps.

Naive Algorithm

A trivial approach to the linear algebra problem is to enumerate all the #G k subsets W = {g 1 , . . . , g k } of G and, for each of them, test if its span covers T . The most efficient way to test this inclusion of spaces is to construct W = Span(W), represented by its basis as a matrix in row echelon form, and then test if each vector t ∈ T reduces to 0 against this matrix. For simplicity's sake, in this paper we focus on the combinatorial complexity of the algorithms, considering that all matrix operations (e.g., computing the dimension of a vector space, checking if a vector is in a vector space, computing the row-echelon form of a matrix) have constant complexity. The complexity of the naive algorithm is thus

#G k .
Example 1 (Cont'd). For the 2 × 3-term polynomial multiplication, #G = 21, k = 5. We therefore have to test 20 349 subsets W.

Improved Algorithm

The main drawback of the naive algorithm is that different subsets of generators W i may be linearly dependent and span the same vector space W = Span(W i ). For instance, for the multiplication of two 3-term polynomials over Q, Montgomery [START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF] gives the following family of solutions:

(a 0 + a 1 X + a 2 X 2 )(b 0 + b 1 X + b 2 X 2 ) = a 0 b 0 (C + 1 -X -X 2 ) (1) +a 1 b 1 (C -X + X 2 -X 3 ) + a 2 b 2 (C -X 2 -X 3 + X 4 ) +(a 0 + a 1 )(b 0 + b 1 )(-C + X) + (a 0 + a 2 )(b 0 + b 2 )(-C + X 2 ) +(a 1 + a 2 )(b 1 + b 2 )(-C + X 3 ) + (a 0 + a 1 + a 2 )(b 0 + b 1 + b 2 )C.
Taking different values for the arbitrary polynomial C ∈ Q[X], we see that any subset of 6 out of the 7 generators

G = {a 0 b 0 , a 1 b 1 , a 2 b 2 , (a 0 + a 1 )(b 0 + b 1 ), (a 0 + a 2 )(b 0 + b 2 ), (a 1 + a 2 )(b 1 + b 2 ), (a 0 + a 1 + a 2 )(b 0 + b 1 + b 2 )
} yields a solution. However, these 7 solutions correspond to the same vector space W , since the 7 generators are linearly dependent:

(a 0 + a 1 + a 2 )(b 0 + b 1 + b 2 ) = (a 0 + a 1 )(b 0 + b 1 ) + (a 0 + a 2 )(b 0 + b 2 ) + (a 1 + a 2 )(b 1 + b 2 ) -a 0 b 0 -a 1 b 1 -a 2 b 2 .
We propose an algorithm that takes advantage of this redundancy by looking directly for the subspaces W -instead of the subsets of generators -that cover our target space T . More formally, we search all the subspaces W of V such that:

(i) T ⊂ W , i.e., the target space T is covered by W ; (ii) Span(W ∩ G) = W , i.e., W is spanned by elements of G; and (iii) dim W = k, i.e., only k generators are needed.

A first remark is that from (i), the target space T is contained in each solution space W . Thus we should look for all the W 's by extending this vector space. Unfortunately, there are a lot of spaces above T , possibly more than there are subsets W of G. Instead, we might consider only those spaces W which are obtained by adding generators to T . This technique has already been used by Oseledets [START_REF] Oseledets | Optimal Karatsuba-like formulae for certain bilinear forms in GF(2)[END_REF] in some heuristic algorithms. In order to use this idea without losing the exhaustiveness of our algorithm, we introduce a new condition:

(ii ) ∃W ⊂ G such that W = T ⊕ Span(W ).

Lemma 1. All subspaces W that verify (i) and (ii) also verify (ii ).

Proof. Let W a subspace verifying (i) and (ii), the result follows directly from Proposition 1 with H being W ∩ G, and F being the free family T ⊂ W .

Proposition 1 ([2, Prop. 3.15]). Let W be a vector space over a field K spanned by a finite set of generators H. Let F be a free family of W . Then there exists a subset J of H such that F ∪ J is a basis of W .

Proof. If F generates W then we take J = ∅.

Let now assume that F does not span W . If we had H ⊂ Span(F), then it would follow that W = Span(H) ⊂ Span(F), which contradicts the assumption. Therefore, there exists h ∈ H such that h / ∈ Span(F), and thus the family F := F ∪ {h} is also free. We then iterate the process with F , until it eventually spans the whole vector space W . Since one cannot choose the same element h ∈ H twice and since H is finite, the process terminates. The set J is then composed of all the generators h added to F when constructing F .

Algorithm 1 Find all the subspaces of dimension k generated by elements of G only and that contain the target space T .

Input: A vector space V , a finite set G of elements of V , a subspace T of V included in Span(G), and an integer k such that dim T ≤ k ≤ rk(G). Output: The set S of all subspaces W of V such that T ⊂ W , Span(W ∩ G) = W , and dim W = k.

1. S ← ∅ 2. procedure expand subspace(W ) 3.

if dim W = k and rk(W ∩ G) = k then 4.
S ← S ∪ {W } 5.

else if dim W < k then 6.

for each g ∈ G \ W do 7.

expand subspace(W ⊕ Span(g)) 8. end procedure 9. expand subspace(T ) 10. return S Thanks to Lemma 1, finding all subspaces W of V verifying conditions (i), (ii), and (iii) can be achieved by first enumerating all subsets W of G such that W = T ⊕Span(W ) verifies conditions (i), (ii ), and (iii), and then keeping only those for which W also verifies condition (ii). This method is implemented in Algorithm 1, which is proven correct in Theorem 1.

Theorem 1 (Correctness of Algorithm 1). For any given k, Algorithm 1 returns all the subspaces W of V verifying conditions (i), (ii) and (iii). In particular, when Algorithm 1 fails to find any solutions for all k under a bound k 0 , then no solutions exist below this bound, and the bilinear rank is greater or equal to k 0 .

Proof. Let W be a subspace of V verifying (i), (ii) and (iii). We first prove that W is included in the output of Algorithm 1. By Lemma 1, W also satisfies (ii ). Thus there is a subset W = {g 1 , ..., g k } of G such that W = T ⊕ Span(g 1 , ..., g k ), and we can choose g 1 , ..., g k to be linearly independent. Therefore, Algorithm 1 will, at some point in the enumeration, consider the candidate subspace W . Since dim W = k and, by condition (ii) we have rk(W ∩ G) = dim W = k, Algorithm 1 will include W in S.

Conversely, let W in the output of Algorithm 1; let us show that conditions (i), (ii), and (iii) are fulfilled. Condition (i) is trivially verified, since Algorithm 1 constructs W by expanding the initial subspace T . Furthermore, since W is in S, we have dim W = k, thus (iii) is also verified. Finally, we have rk(W ∩ G) = k, which implies Span(W ∩ G) = W , and thus W also fulfills condition (ii).

Example 1 (Cont'd). For our running example, Algorithm 1 immediately shows that no solutions exist with k = 4 over F 2 . Indeed, we have dim T = 4, thus dim W = 4 at the very first call of expand subspace, but the rank of W ∩ G is only 3.

Example 5 (

Step-by-step illustration of Algorithm 1). Take K = Q, V = Q 3 , T = {t 0 , t 1 } with t 0 = (1, 1, 1) and t 1 = (1, 1, 0); and G = {g 0 , g 1 , g 2 , g 3 } with g 0 = (1, 0, 0), g 1 = (0, 1, 0), g 2 = (0, 0, 1), and g 3 = (1, 0, 1), as depicted in Figure 1. In this example, for clarity's sake, we denote the recursive depth of the current call by a parenthesized superscript: e.g., W (i) represents the subspace W at recursive depth i, starting at i = 0.

Let us first consider the case k = 2. The algorithm initializes W (0) to T = Span(t 0 , t 1 ), which is already of dimension 2. It then computes W ∩ G = {g 2 }, which is of rank 1. Therefore, there are no solutions for k = 2.

Consider now the case k = 3. Starting again with W (0) = T , since dim T = 2, the algorithm will then have to expand W (0) by adding a linear independent vector g from

g 0 g 3 g 1 t 1 g 2 t 0 0 G \ W (0) = {g 0 , g 1 , g 3 }.
Taking g = g 0 and adding it to W (0) , the algorithm recurses down with W (1) = T ⊕ Span(g 0 ). We have W (1) ∩ G = {g 0 , g 1 , g 2 , g 3 }, which has rank 3, meaning that this W (1) is a solution subspace, which is then added to the set S. The algorithm can now reiterate the process for the subspaces W (0) ⊕ Span(g 1 ) and W (0) ⊕ Span(g 3 ), finding both of them to also be solutions.

Note that there is no point to try the algorithm for k > 3: since dim V = 3, the algorithm will never find subspaces W of V of dimension k.

Remark also that this simple example illustrates rank leaps for W ∩G: if rk(W (0) ∩G) = 1, when taking W (1) = W (0) ⊕ Span(g 0 ), the rank of W (1) ∩ G jumps directly to 3. In fact, since T cannot usually be generated by elements of G only, we have dim W (0) > rk(W (0) ∩G) at the beginning of the algorithm. Therefore, since dim W is only incremented by 1 at each recursive call, rank leaps between W (i) ∩ G and W (i+1) ∩ G gradually close the gap between dim W and rk(W ∩ G).

Example 6 (Evaluation-interpolation schemes). In the case of the n-term polynomial multiplication over a field K, evaluating the resulting polynomial at a point κ ∈ K only involves a linear combination of the target coefficients t i . Therefore, the rank-1 bilinear forms ( i κ i a i )( j κ j b j ) corresponding to the evaluation of the two input polynomials at κ are also in T . This is also the case when evaluating "at infinity", i.e., considering the leading coefficient a n-1 b n-1 . Therefore, since these rank-1 bilinear forms are linearly independent (by Vandermonde determinant), if #K ≥ 2n -2, then we have rk(T ∩ G) = 2n -1 = , and T itself is a solution subspace.

Implementation Issues

Algorithm 1 relies extensively on operations on vector spaces, such as computing the rank of the intersection W ∩ G (line 3), or excluding elements of W from G (line 6). In order to perform these operations efficiently, each subspace W is represented in the algorithm by the r × nm matrix of its basis in row echelon form, where r = dim W . Note that we also denote by W the matrix in row echelon form corresponding to the subspace W .

Computing W G = Span(W ∩ G) is done as follows: first reduce every vector g ∈ G by the matrix W ; if g reduces to the null vector, then g ∈ W . For each such g, reduce it by the current matrix W G (initially empty), and if the reduced vector is not null, add it to W G . Once W G is constructed, the rank of W ∩ G is then given by the number of rows of the matrix.

Enumerating all the generators g ∈ G \ W (line 6) can be achieved using a similar technique. However, it is more efficient to maintain the set H of all generators g ∈ G reduced by W , and pick only from these reduced generators for the recursive calls. Indeed, it might very well be the case that two generators g and g ∈ G are in fact the same vector after reduction by W . In that case, the two subspaces W ⊕ Span(g) and W ⊕ Span(g ) are identical, and recursing down the algorithm for both of them would just be a waste of time. Maintaining a set H of reduced generators prevents this to happen, as g and g will then correspond to a single reduced generator in H. Furthermore, at recursion depth i, since all elements g of H (i) are reduced by W (i) , when recursing down the algorithm:

constructing W (i+1) = W (i) ⊕ Span(g) only requires inserting the row vector g into the matrix W (i) ; and reducing each element h ∈ H (i) by W (i+1) to construct H (i+1) only requires reducing h by Span(g), as it is already reduced by W (i) .

Example 1 (Cont'd). On our running example, T = Span(T ) is the 4 × 6 matrix -which is already in row echelon form -and G consists of the 21 generators from Section 3. Let us follow the execution of the algorithm in the case k = 5.

We first construct H (0) by reducing all generators of G by W (0) = T . Only 3 unique reduced generators remain:

H (0) =   a 0 b 0 a 0 b 1 a 0 b 2 a 1 b 0 a 1 b 1 a 1 b 2 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 1 0   a 1 b 0 a 1 b 1 a 1 (b 0 + b 1 )
Since dim W (0) = 4 < k, these generators in H (0) are then enumerated, and the procedure expand subspace is called recursively for each of them. For instance, when taking g = a 1 b 0 , the subspace W (1) is constructed by inserting g into the matrix W (0) :

W (1) =       a 0 b 0 a 0 b 1 a 0 b 2 a 1 b 0 a 1 b 1 a 1 b 2 1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 1       a 0 b 0 a 0 b 1 + a 1 b 0 a 0 b 2 + a 1 b 1 g = a 1 b 0 a 1 b 2
Computing W (1) ∩ G yields a set of 9 generators and of rank 5; W (1) is thus a solution subspace. Finally, continuing the enumeration for g = a 1 b 1 then g = a 1 (b 0 + b 1 ), the algorithm will find two other solution subspaces of rank 5.

From Solution Subspaces back to Formulae

Given a solution subspace W verifying conditions (i), (ii), and (iii), we still need to retrieve the corresponding formulae. Note that since T ⊂ W , each formula corresponds to a basis of W using only generators from G. A simple solution is thus to first compute the set W ∩ G, then enumerate all subsets of k linearly independent vectors in W ∩ G.

For every basis of W , obtaining the corresponding formula is then just a matter of finding the coordinates of the vectors of T in this basis, which will give us the corresponding linear combinations of generators required to compute them.

Example 1 (Cont'd). In our running example, each of the 3 solution subspaces yields 54 different 5-subsets of G that span T . We thus find a total of 162 formulae for evaluating 2 × 3-term polynomial multiplications over binary fields.

Complexity Analysis

Algorithm 1 allows us to drastically reduce the number of operations required to find all the formulae for evaluating a given bilinear map thanks to the two following key ideas:

1. searching for subspaces W of Span(G) instead of subsets W of G; and 2. constructing these subspaces starting from T instead of {0}.

The algorithm consists in expanding T to a vector space of dimension k by adding k -vectors from G, where = dim T . The combinatorial complexity of the algorithm is thus bounded by: #G k -.

Note however that this complexity bound does not reflect the gain due to the first idea -searching for subspaces instead of subsets -since that gain depends on the particular problem to be solved, and cannot be expressed simply. In Algorithm 1, once say g 1 , g 2 , . . . have been added to the initial W (0) = T , this idea will prune the remaining set of independent generators G \W (without this idea, instead of considering all g ∈ G \W at line 6 of Algorithm 1, we would consider all g ∈ G, which would be very similar to the naive algorithm). Note however that Algorithm 1 does not guarantee that all subspaces W explored and found to be solutions are different. Duplicates have to be detected and removed upon adding solution subspaces to the set S.

In practice, this combinatorial complexity is usually not attained, as can be seen in the experimental results reported in Section 5.

Example 1 (Cont'd).

For the 2×3-term polynomial multiplication over F 2 , with #G = 21, k = 5, and = 4, the complexity estimate predicts 21 subspaces W to consider, whereas in practice, only 3 of these subspaces need be explored by the algorithm. This has to be compared with the 20 349 subsets explored by the naive algorithm.

Special Case of Symmetric Bilinear Maps

An n × n bilinear form γ : a, b → i,j γ i,j a i b j is said to be symmetric if γ i,j = γ j,i for all 0 ≤ i, j < n. Similarly, an n × n bilinear map φ is said to be symmetric if the bilinear forms corresponding to its coefficients are all symmetric. For instance, it is the case for the n-term polynomial multiplication.

Example 7 (3-term polynomial multiplication). Given

A = a 2 X 2 + a 1 X + a 0 and B = b 2 X 2 + b 1 X + b 0 ∈ K[X], the target set for the product A • B is T = {a 0 b 0 , a 0 b 1 + a 1 b 0 , a 0 b 2 + a 1 b 1 + a 2 b 0 , a 1 b 2 + a 2 b 1 , a 2 b 2 },
all of whose bilinear forms are symmetric.

In the case where the target set T is composed only of symmetric forms, an idea to accelerate the search for formulae is to restrict G to be the set of rank-1 symmetric bilinear forms, which are of the form ( i α i a i )( j α j b j ). Indeed, while there are (#K) n 2 -1 nonzero bilinear forms of rank 1 in B(K n , K n ; K), only (#K) n -1 of them are symmetric, which yields a huge speedup.

However, the symmetric rank-1 forms alone may fail to produce optimal formulae for symmetric maps: for instance, a 0 b 1 +a 1 b 0 cannot be computed with 2 symmetric products over F 2 . (For an introduction on the topic of decomposing symmetric maps as sums of symmetric tensors see for instance [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF].) Therefore, restricting to symmetric forms has to be considered as a heuristic, and cannot be used to prove lower bounds on the bilinear rank of symmetric maps.

Experimental Results

We present here some results obtained for various instances of the bilinear rank problem using a multi-threaded C implementation of Algorithm 1. The reported timings correspond to the total execution time on a single core of a 2.2 GHz Xeon L5640. We reproduce or improve some known results [START_REF] Cenk | On multiplication in finite fields[END_REF][START_REF] Chung | Asymmetric squaring formulae[END_REF][START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF] but also find the complete list of solutions for each problem or prove lower bounds thanks to Theorem 1.

In all the following tables, only one value for k is reported for each problem. It is the smallest possible value for which solutions were found or, if no solutions were found, it is the largest value which was tried with our algorithm. In both cases, it means that there are no solutions for smaller values of k. Rows without a value for k were deemed too computationally expensive and were not tried at all with our implementation.

Under the heading "# of tests" are reported the number of subspaces W for which the algorithm had to check whether rk(W ∩ G) = k or not. It corresponds to the number of leaves in the tree of recursive calls. The column "# of solutions" reports the number #S of solution subspaces returned by Algorithm 1. Finally, we indicate in the column "# of formulae" the corresponding number of formulae, or ∞ when this number is really large.

Finally, where applicable and when an exhaustive search using all rank-1 bilinear forms for G was too expensive, we enumerated subspaces using only symmetric generators. This is indicated by a "(Sym.)" mention next to the cardinality of G.

Thus, Table 1 gives experimental results for n × m-term polynomial multiplication over F 2 . For instance, the 6 × 6 row indicates that there are no formulae with only k = 14 full multiplications over F 2 ; however, symmetric formulae exist for k = 17. Similarly, the 5 × 5 row shows that Montgomery's bound of k = 13 is optimal over F 2 [START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF]. Table 2 gives similar results over F 3 . For example, the 5 × 5 row with no solution for k = 11 proves that the bound M 3 (5) = 12 from [6, Table 2] is optimal.

Table 3 corresponds to small extensions of F 2 or F 3 , where we consider the multiplication of two elements in polynomial basis. This is relevant for example for multiplication of matrices over F 2 e [START_REF] Albrecht | The M4RIE library for dense linear algebra over small fields with even characteristic[END_REF]. In particular the rows F 2 4 with k = 9 and F 3 4 with k = 6 confirm the values from [7, Table 1].

Furthermore, according to the literature, the best known formula for computing the product of two elements over the finite field F 3 5 uses 12 full multiplications [START_REF] Cenk | On multiplication in finite fields[END_REF]. As indicated in Table 3, our algorithm found 121 symmetric formulae using only 11 full multiplications. One such formula is given in Algorithm 2.

Finally, Table 4 considers the product of two n-term polynomials modulo X n and X n -1 over F 2 -as in [START_REF] Oseledets | Optimal Karatsuba-like formulae for certain bilinear forms in GF(2)[END_REF] -and also over F 3 . This proves the tensor rank from [START_REF] Oseledets | Optimal Karatsuba-like formulae for certain bilinear forms in GF(2)[END_REF] is the optimal one for n = 2, 3, 4 over F 2 modulo both X n and X n -1. One should note that over F 3 , computing a product modulo X 4 -1 can be done with the same number of products (5) than modulo X 3 -1.

Conclusion

Long considered to be "magical", Karatsuba's and Strassen's formulae can be understood in a unified framework: the bilinear rank problem. Small instances of this NP-hard problem can be tackled by exhaustive search methods, as proposed by Montgomery [START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF]. Our improved algorithm proved that many known formulae from the literature, some of which discovered by ad-hoc methods, are optimal. Not relying on heuristic methods, the algorithm we presented here finds all the formulae using the optimal number of full multiplications, and gives lower bounds. For example, for the product of 3-term polynomials over F 2 , in addition to the 7 possible formulae from Eq. ( 1) already found by Montgomery in [START_REF] Montgomery | Five, six, and seven-term Karatsuba-like formulae[END_REF], we found two new asymmetric formulae, the first one using the generators It would be interesting to give a simple mathematical explanation for these new formulae.

We were also able to improve the bound from [START_REF] Cenk | Multiplication of polynomials modulo x n[END_REF] for the multiplication in F 3 5 from 12 to 11 multiplications, using again a completely generic method.

If one wants to go further, one has to either speed up the computations or to use more heuristic or proven techniques borrowed to the algebraic complexity community: restricting the set G of products, e.g., to symmetric products for polynomial multiplication, enlarging the target set T by imposing that some products have to occur in the formula or using the symmetries of the problem.

written as W = T ⊕ g 1 , ..., g k-dim T . One recognizes the idea of replacing condition (ii) by (ii ) in Section 4.2. Oseledets did not realize this change still performs an exhaustive search, as demonstrated by Theorem 1, since he says on page 2061 "Partial exhaustive search -search only for the complement space C of R".

a 1 b 2 ,

 2 and g 4 = (a 0 + a 1 )(b 0 + b 1 ). Since A • B has 4 terms, if the base field K contains at least 3 elements, evaluating A • B at those elements and at infinity -i.e., the a 1 b 2 product -enables to recover A • B by Lagrange interpolation. However, if K = F 2 , are 5 products optimal?

Fig. 1 .

 1 Fig. 1. The sets G and T from Example 5.

g 1 :

 1 = a 0 b 0 , g 2 := a 2 b 2 , g 3 := (a 0 + a 1 )(b 0 + b 2 ), g 4 := (a 0 + a 2 )(b 1 + b 2 ), g 5 := (a 1 + a 2 )(b 0 + b 1 ), g 6 := (a 0 + a 1 + a 2 )(b 0 + b 1 + b 2 ), with (a 0 + a 1 X + a 2 X 2 )(b 0 + b 1 X + b 2 X 2) being equal to:g 1 + (g 2 + g 3 + g 5 + g 6 )X + (g 3 + g 4 + g 6 )X 2 + (g 1 + g 4 + g 5 + g 6 )X 3 + g 2 X 4, and the second one using the generators a 0 b 0 , a 2 b 2 , (a 0 + a 1 )(b 1 + b 2 ), (a 0 + a 2 )(b 0 + b 1 ), (a 1 + a 2 )(b 0 + b 2 ), (a 0 + a 1 + a 2 )(b 0 + b 1 + b 2 ).

Table 1 .

 1 Experimental results for n × m polynomial multiplication over F2.

	n × m	#G	k	# of tests	# of solutions formulae time [s] # of Calculation
	2 × 2		9 3	1		1		1	0.00
	3 × 2		21 5	2		3	162	0.00
	3 × 3		49 6	9		3		9	0.00
	4 × 2		45 6	5		4	108	0.00
	4 × 3		105 8 700	33	423	0.00
	4 × 4		225 9 6.60 • 10 3		4		4	0.03
	5 × 2		93 8	56	28	790 272	0.00
	5 × 3		217 10 1.46 • 10 5	366	48 195	0.51
	5 × 4		465 12 3.13 • 10 8	4 113	66 153	2.86 • 10 3
	5 × 5		961 13 9.65 • 10 9	27		27	2.28 • 10 5
	6 × 2		189 9 250	64	1 404 928	0.00
	6 × 3		441 11 2.05 • 10 6		3	243	11.5
	6 × 4		945 13 7.69 • 10 9		9		15	1.62 • 10 5
	6 × 5	1 953 14 2.01 • 10 11	-		-	9.97 • 10 6
	6 × 6	3 969 14 4.37 • 10 9 (Sym.) 63 17 8.08 • 10 6	-	6	-	54	6.03 • 10 5 17.7
	7 × 2		381 11 9.14 • 10 3	960		∞	0.07
	7 × 3		889 13 2.52 • 10 9	87	63 423	3.66 • 10 4
	7 × 4	1 905 14 1.47 • 10 11	-		-	6.34 • 10 6
	7 × 7 (Sym.) 127 22 3.38 • 10 12	2 618	19 550	1.59 • 10 7
	8 × 2		765 12 7.80 • 10 4	4 096		∞	0.75
	8 × 3	1 785 14 5.27 • 10 10	-		-	-

Table 2 .

 2 Experimental results for n × m polynomial multiplication over F3.

	n × m	#G	k	# of tests	# of solutions formulae time [s] # of Calculation
	2 × 2		16 3	1	1	4	0.00
	3 × 2		52 4	1	1	1	0.00
	3 × 3		169 6	24	22	1 493	0.00
	4 × 2		160 6	9	13	38 880	0.00
	4 × 3		520 7 164	12	48	0.00
	4 × 4	1 600 9 4.11 • 10 5	726	50 640	14.9
	5 × 2		484 7	24	36	93 312	0.00
	5 × 3	1 573 9 2.81 • 10 5	1 116	94 629	9.33
	5 × 4	4 840 10 4.75 • 10 6	48	768	1.01 • 10 3
	5 × 5	14 641 11 4.89 • 10 7 (Sym.) 121 12 3.93 • 10 4	-31	-6 460	4.02 • 10 4 0.14
	6 × 2	1 456 8	69	81	104 976	0.01
	6 × 3	4 732 10 3.24 • 10 6	240	4 272	566
	6 × 4	14 560 11 4.55 • 10 7	-	-	3.31 • 10 4
	6 × 5	44 044 12 4.58 • 10 8	-	-	1.31 • 10 6
	6 × 6 (Sym.) 364 15 2.37 • 10 8	4	1 024	3.79 • 10 3
	7 × 2	4 372 10 2.27 • 10 4	10 530	∞	2.84
	7 × 3	14 209 11 3.15 • 10 7	-	-	1.84 • 10 4
	7 × 4	43 720 12 4.16 • 10 8	-	-	1.03 • 10 6
	7 × 7 (Sym.) 1 093 17 2.69 • 10 10	-	-	1.50 • 10 6
	8 × 2	13 120 11 2.01 • 10 5	85 293	∞	53.6
	8 × 3	42 640 12 2.90 • 10 8	-	-	5.46 • 10 5

Table 3 .

 3 Experimental results for multiplication over small extensions of F2 and F3.

	Finite field	#G	k	# of tests	# of solutions formulae time [s] # of Calculation
	F 2 2		9 3	3	3	3	0.00
	F 2 3		49 6 7.03 • 10 3	105	147	0.01
	F 2 4		225 9 2.57 • 10 9	2 025	2 025	1.13 • 10 4
	F 2 5	961 9 3.10 • 10 10 (Sym.) 31 13 3.49 • 10 6	-2 015	-2 015	8.11 • 10 5 6.24
	F 2 6	(Sym.) 63 15 2.21 • 10 10	21	21	6.63 • 10 4
	F 2 7	(Sym.) 127 15 1.34 • 10 12	-	-	6.17 • 10 6
	F 3 2		16 3	3	4	16	0.00
	F 3 3		169 6 2.42 • 10 5	11 843	105 963	1.08
	F 3 4	1 600 8 2.27 • 10 11 (Sym.) 40 9 1.10 • 10 5	-234	-615 240	1.08 • 10 7 0.45
	F 3 5	(Sym.) 121 11 2.66 • 10 9	121	121	1.45 • 10 4
	F 3 6	(Sym.) 364 12 3.01 • 10 12	-	-	4.50 • 10 7

Table 4 .

 4 Experimental results for the multiplication of two n-term polynomials in Fp[X]/(X n ) and Fp[X]/(X n -1), with p = 2 and 3.

	Ring	n	#G	k	# of tests	# of solutions formulae time [s] # of Calculation
		2		9 3	3	3		10	0.00
		3		49 5 590	12		40	0.00
		4		225 8 5.17 • 10 7	1 440	9 248	230
	F2[X]/(X n )	5	961 9 2.66 • 10 10 (Sym.) 31 11 3.64 • 10 5	-112	-736	6.70 • 10 5 0.48
		6 (Sym.) 63 14 2.63 • 10 9	384	2 816	7.66 • 10 3
		7 (Sym.) 127 15 1.16 • 10 12	-	-	5.46 • 10 6
		2		9 3	3	3		10	0.00
		3		49 4	21	3		3	0.00
		4		225 8 2.69 • 10 7	1 440	9 248	124
	F2[X]/(X n -1)	5	961 9 1.39 • 10 10 (Sym.) 31 10 7.46 • 10 4	-25	-	25	3.65 • 10 5 0.09
		6 (Sym.) 63 12 2.33 • 10 7	31	148	50.0
		7 (Sym.) 127 13 2.55 • 10 9	1		49	1.24 • 10 4
		2		16 3	4	4		39	0.00
		3		169 5 7.94 • 10 3	90	1 539	0.07
	F3[X]/(X n )	4	1 600 7 5.54 • 10 8 (Sym.) 40 8 3.17 • 10 5	-252	-40 095	3.22 • 10 4 0.14
		5 (Sym.) 121 10 1.45 • 10 8	243	13 122	2.28 • 10 3
		6 (Sym.) 364 11 4.79 • 10 10	-	-	8.22 • 10 5
		2		16 2	1	1		1	0.00
		3		169 5 4.45 • 10 3	90	1 539	0.04
	F3[X]/(X n -1)	4	1 600 5 767	4		16	0.07
		5 (Sym.) 121 10 8.74 • 10 7	234	615 240	1.39 • 10 3
		6 (Sym.) 364 10 3.37 • 10 8	9	2 025	1.68 • 10 4

Acknowledgements. The authors would like to thank Marc Mezzarobba for the interesting and fruitful discussions we had on the subject, especially for his knowledge of the many results on this topic.

Appendix

We compare here Algorithm 1 to Oseledets' work. The main similarity is that both methods use a linear algebra setting and search for solution spaces W rather than solution sets W of generators. In terms of differences, we note that Algorithm 1 is proven, making it a tool for proving lower bounds (cf. Section 5). Oseledets' work presents a main algorithm called Minbas, a subroutine called Rank-1 and some more tricks.

Minbas starts by computing the rank of all #K dim T -1 elements of T (denoted C in [START_REF] Oseledets | Optimal Karatsuba-like formulae for certain bilinear forms in GF(2)[END_REF]), where the rank of a bilinear form is the minimal k such that t = g 1 + g 2 + • • • + g k with g i ∈ G. Minbas continues by constructing a basis of T by adding in order linearly independent elements of rank 1, 2 and so on. In our setting, the rank-1 elements from [START_REF] Oseledets | Optimal Karatsuba-like formulae for certain bilinear forms in GF(2)[END_REF] are elements of the target space which are also in the set G of generators. Since Algorithm 1 starts by W = T , those forms are automatically included in our algorithm. From this point on, the two algorithms diverge since Algorithm 1 explores spaces W of higher dimension whereas Minbas adds forms of rank 2 or more.

The second algorithm of Oseledets' work is called Rank-1 and can be used as a subroutine in algorithms solving the bilinear rank problem. Given a candidate space W of dimension k lying above T , the Rank-1 procedure tests if W contains a basis made out of elements in G by writing the lines of W and G in row echelon form. In our setting this is equivalent to the test rk(W ∩ G) = k.

Without giving a general description of the method, nor analyzing its correctness, Oseledets implements the idea which we presented in Section 4.2. For example, in the problem of 5-term polynomial multiplication we have dim T = 9 and we search solution spaces of dimension k = 13; Oseledets notes that "we have to add 3 [there is a typo here, read 4] matrices to the current basis. There are 31 -3 = 28 possible symmetric rank-1 matrices to add". The implicit idea is that we restrict our search to spaces which can be


Heidi Ruiz

Pedro Cobo

Thomas Dupont

Bruno Martin

Philippe Leclaire

Acoustic properties of plates with unevenly distributed macro perforations backed by woven meshes

Keywords: numbers: 43.55.Ev, 43.50.Gf, 43.20.Hq

A hybrid model describing the acoustic properties of plates with macroperforations that can be unevenly distributed on the plate surface and backed by woven or precision woven meshes with microscopic perforations is proposed. The plate perforations may be of circular or rectangular shapes.

Since the perforated plate may not necessarily be considered as an equivalent fluid, its surface impedance is calculated by the Maa model [Noise Control Eng. J. 29, 77-84 (1987)], whereas the Johnson-Champoux-Allard model [J. Appl. Phys. 70, 1975-1979 (1991)] is used for the mesh, considered as an equivalent fluid. A simple model for the elementary cell of the mesh structure is proposed in order to calculate a parameter that can be considered as the thermal length. It is shown that, this thermal length can be smaller than the length considered as the viscous length. An effective airflow resistivity is introduced to account for the increase of particle velocity through the mesh placed behind the carrying macroperforated plate and is used in the transfer matrix approach to obtain the impedance of the whole multilayer system. The hybrid model seems to represent a good approach of this multilayer system. The theoretical predictions are compared with experimental measurements.

I. INTRODUCTION

The acoustical behavior of perforated panels has been an important research topic over the last decades. The analytical model proposed by Maa [START_REF] Maa | Microperforated-panel wideband absorbers[END_REF] is frequently used to predict the acoustic impedance of a microperforated panel (MPP) with circular perforations. The equation for the impedance combines the model of sound propagation in narrow tubes studied by Rayleigh and Crandall [START_REF] Rayleigh | The Theory of sound[END_REF][START_REF] Crandall | Theory of vibrating systems and sound[END_REF] , and includes the end corrections done by Ingard [START_REF] Ingard | On the theory and design on acoustic resonators[END_REF] . In the useful audible frequency range, the best performance for perforated panels is obtained for submillimetric holes. As recommended by Maa 1 , the diameter of perforation -to -thickness ratio should be close to one to obtain optimal absorption properties. Manufacturing panels with submillimetric holes and submillimetric thicknesses can be a difficult task in practical acoustic engineering applications. Furthermore, drilling a large amount of submillimetric holes in a panel (usually done by laser technology) can be expensive. An alternative method to reduce the amount of perforations while keeping the porosity equal to that of an MPP is to use slits [START_REF] Maa | Theory of microslit absorbers(in Chinese)[END_REF][START_REF] Randeberg | Perforated panel absorbers with viscous energy dissipation enhanced by orifice design[END_REF] . Manufacturing a panel bearing slits (MSP) can be less expensive than an MPP.

A perforated plate can in some cases also be modeled as a rigid frame porous medium described by the well established Johnson-Champoux-Allard (JCA) theory [START_REF] Allard | Propagation of Sound in Porous Media[END_REF][START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF][START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF] that depends on five physical parameters: the thermal Λ and viscous Λ characteristic lengths, the porosity φ, the tortuosity α ∞ , and the static air flow resistivity σ. Atalla and Sgard 10 presented an equivalent fluid approach that considers these parameters. The analytical models of Maa and the equivalent fluid approach agree in what concerns the impedance inside the perforations. Nevertheless, the end corrections differ. While Maa [START_REF] Maa | Microperforated-panel wideband absorbers[END_REF][START_REF] Maa | Theory of microslit absorbers(in Chinese)[END_REF] introduces the end correction as an additional complex impedance with resistive and reactive components, Atalla and Sgard 10 introduce the end correction by means of the geometric tortuosity α ∞ as a multiplying factor of the effective density ρe .

To increase the sound attenuation, perforated panels can be backed by porous materials or textiles [START_REF] Kang | Predicting the absorption of open weave textiles and micro-perforated membranes backed by an air space[END_REF][START_REF] Pfretzschner | Microperforated Insertion Units: An alternative strategy to design microperforated panels[END_REF][START_REF] Chevillote | Controlling sound absorption by an upstream resistive layer[END_REF] . Combining a carrying plate with a commercial woven mesh makes it possible to design an absorber with a performance similar to that of a single-layer MPP. Indeed, Pfretzschner et a) Electronic address: heidi@caend.upm-csic.es al. [START_REF] Pfretzschner | Microperforated Insertion Units: An alternative strategy to design microperforated panels[END_REF] proposed a design procedure that consisted on combining a perforated panel and a micrometric mesh with appropriate constitutive parameters, which they named a microperforated insertion unit (MIU). The impedance of the MIU combination resulted from the addition of the acoustic impedances of the perforated plate and of the wire mesh divided by the open area fraction of the plate φ. The impedance of the system was then obtained by the impedance transfer method [START_REF] Cobo | Comparison of models describing double layer microperforated absorbers[END_REF] .

Recently, Chevillote [START_REF] Chevillote | Controlling sound absorption by an upstream resistive layer[END_REF] studied the role of a resistive layer upstream on a porous medium.

The resistive layer, was then modeled using the JCA equivalent fluid model [START_REF] Allard | Propagation of Sound in Porous Media[END_REF][START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF][START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF] . The acoustic impedance of such resistive layer was obtained using the transfer matrix method [START_REF] Allard | Propagation of Sound in Porous Media[END_REF] . Chevillote used a 2D microscopic characterization with help of a picture processing to determine the hydraulic radius and deduce the thermal characteristic length of the textile.

In this article, we study a multilayer system made of macroperforated metal plates backed by microperforated woven textile meshes. This solution involving fairly high perforation radii for the carrying plate and widely available textiles can significantly decrease the costs of manufacturing of the acoustic system. The prefixes "micro" or "macro" are related to the relative importance of the radius with respect to the viscous skin depth inside the perforation.

In the present study, a new approach involving a simple model for the elementary cell of the woven textile structure is proposed to estimate more precisely the thermal characteristic length. Also, a hybrid model capable of describing non-homogenously distributed macroperforated panels backed by woven textiles is proposed. The carrying plate is modeled using the analytical approach 1,5 of Maa whereas the JCA model is used in the rigid frame approximation for the mesh.

The hybrid model is valid for carrying panels with circular or slit-like macroperforations regardless of their distribution over the surface. Another original aspect developed in this article is the introduction of an effective resistivity. The flow resistivity of the mesh is calculated by means of the shape factor described by Stinson and Champoux [START_REF] Stinson | Propagation of sound and the assignment of shape factores in model porous materials having simple pore geometries[END_REF] . In addition, whereas in Ref. 12 the fact that bonding the membrane to the carrying plate modified the perforation rate of the screen, it is shown in the present study that the modification can be accounted for through the use of an effective resistivity.

The description of the theory used to obtain the impedance of the carrying plate is presented in section II. The calculation of the thermal characteristic length of microperforated woven meshes is described in section III. Section IV presents the hybrid model using a combination of the models of Maa and JCA in the context of the transfer matrix method and also the derivation of the effective flow-resistivity. Experimental validation is presented in the last section.

II. PERFORATED PLATE

The surface impedance of a perforated panel can be written as the sum of two contributions [START_REF] Maa | Microperforated-panel wideband absorbers[END_REF][START_REF] Maa | Theory of microslit absorbers(in Chinese)[END_REF] : the impedance at the interior of the perforation and the impedance at the borders which accounts for the end corrections:

Z 1 = √ 2µy φd + jωρ 0 φ Rψ(ξ)r + t [1 + S ] , (1)
where µ is the air kinematic viscosity coefficient, ρ 0 the air density, ω the angular frequency, t the panel thickness, r the perforation radius, φ the perforation ratio, and y = d/ 4µ/ρω with d the perforation diameter. R and S are parameters that depend on the geometry of the perforations, in particular on the length of the slit l. In the case of periodic hole arrangement, ψ(ξ) is a correction of the reactance term of the end correction that takes into account the interaction between holes when the perforation ratio is high [START_REF] Melling | The acoustic impedance of perforates at medium and high sound pressure lev-els[END_REF][START_REF] Tayong | Experimental investigation of holes interaction effect on the sound absoprtion coefficient of micro-perforated panels under high and medium sound levels[END_REF] .

The end corrections studied by Rayleigh and Ingard [START_REF] Rayleigh | The Theory of sound[END_REF][START_REF] Ingard | On the theory and design on acoustic resonators[END_REF] are included in the first and second term of Eq. ( 1). The influence of a mass of air contained in a cylindrical volume for which the apparent length is greater than the plate thickness is included in the numerator of the second term. The third term represents the impedance of the perforations which has been studied by Crandall 3 for an array of tubes and by Ingard [START_REF] Ingard | On the theory and design on acoustic resonators[END_REF] for an elliptic aperture. The velocity profile of the flow through the holes is determined by the term y. The values of R, S and Ψ(ξ) are shown in Table I where J 0 , J 1 are Bessel functions of the first kind and order 0 and 1, respectively and F (e) is an elliptic integral with e = 1 -(r/l) [START_REF] Rayleigh | The Theory of sound[END_REF] .

F (e) = π 2       1 + 1 2 2 e 2 + 1 × 3 2 × 4 2 e 4 + 1 × 3 × 5 2 × 4 × 6 2 e 6 + ...       , (2)
The equation of impedance described by Randeberg in Ref. 6, includes the interaction between perforations as well as the resistive and reactance end corrections deduced for the modified depth of the slit [START_REF] Ruiz | Proposal of alternative designs for microperforated panels[END_REF] . For the present study, the equations proposed by Maa will be used to model a macroperforated carrying plate.

III. MICROMETRIC MESH

The absorption of textiles and microperforated membranes has been calculated by Kang and Fuchs using the equation of impedance given by Maa. In a recent study, Atalla and Sgard considered a plate or screen as an equivalent fluid using the rigid frame porous model developed by

Johnson-Champoux-Allard [START_REF] Allard | Propagation of Sound in Porous Media[END_REF][START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF][START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF] . The model by Maa and that of Atalla and Sgard are similar in what concerns the impedance inside the perforations. However, in the equivalent fluid approach, the end corrections are introduced by means of the tortuosity allowing to define the boundary conditions in multilayer configurations. Latest works apply this model for woven or non-woven fabrics [START_REF] Chevillote | Controlling sound absorption by an upstream resistive layer[END_REF][START_REF] Jaouen | Acoustical characterization of perforated facings[END_REF] .

In many cases, a homogeneous cross section of perforations is considered (for example cylindrical or square cross section). Nevertheless, as the woven meshes studied here exhibit more complex squared pore geometries, a more specific approach is needed (three dimensions approach is required for some parameter characterizations).

In the JCA model, the tortuosity, resistivity, porosity and characteristic lengths of the pores can be obtained [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF] or measured for each specific equivalent fluid. The visco-thermal effects associated with a wave impinging onto a porous medium can be modeled introducing a complex dynamic density i.e, the effective density ρ e (ω) (see ref . 7).

ρe (ω) = α ∞ ρ 0          1 + σφ jωρ 0 α ∞ 1 + jω4α 2 ∞ µρ 0 σ 2 Λ 2 φ 2          , (3)
where α ∞ is the tortuosity, ρ 0 the air density, σ the flow resistivity, φ the porosity, which is considered equivalent to the perforation ratio, ω the angular frequency, µ is the air kinematic viscosity coefficient and Λ is the viscous length given by [START_REF] Johnson | Theory of dynamic permeability and tortuosity in fluid-saturated porous media[END_REF]

2 Λ = s v 2 i (r w )dS v v 2 i (r)dV . ( 4
)
The integrations of the numerator and the denominator of Eq.( 4) are over the walls of the pore perforation and over the pore volume respectively.

Equation (3) indicates that the parameters α ∞ and σ are important factors in the calculation of the effective density. While the flow resistivity depends on the pore shape, the tortuosity depends on the pore shape and the layer arrangement. For a screen with air at both sides,

α ∞ = 1 + (2ε/t)
where ε represents a correction length.

The classical MPP's formulation [START_REF] Maa | Theory and design of microperforated panel sound absorbing constructions[END_REF] and the equivalent fluid MPP's approach 10 take into account the increase of kinetic energy caused by the sound radiation effect at the end openings of perforations and the distortion effect of the acoustic flow at the panel surface. As a result, the MPP's reactance must be supplemented by an additional mass by way of the length correction ε on the plate thickness, in the imaginary part of the MPP's impedance [START_REF] Ingard | On the theory and design on acoustic resonators[END_REF][START_REF] Allard | Propagation of Sound in Porous Media[END_REF] . In the case of a single perforation, the classical approach employs an analogy with the radiation impedance of a rigid piston with a circular cross section flanged in an infinite rigid wall. For this approach the end termination of perforation is considered to have a sharp edge. In the present study, however, the end termination of the screen perforation has rounded edges. This kind of edge tends to diminish the length correction [START_REF] Dalmont | Radiation of tubes with different flanges: numerical and experimental investigations[END_REF] . Moreover, when there are several perforations in close proximity, one must consider that the interaction between perforations tends to diminish the length correction [START_REF] Ingard | On the theory and design on acoustic resonators[END_REF][START_REF] Randeberg | Perforated panel absorbers with viscous energy dissipation enhanced by orifice design[END_REF] . For the meshes treated in this study, the flange between two perforations is small compared to their size, which causes a strong interaction between the rounded edges of the pore. The length correction is thought to be very small and is therefore disregarded. As a result, the effective tortuosity is equal to 1.

The flow resistivity σ has been described by Stinson and Champoux [START_REF] Stinson | Propagation of sound and the assignment of shape factores in model porous materials having simple pore geometries[END_REF] for simple pore shapes as

α = 4k 0 µα 2 ∞ h 2 r φ , (5)
where, h r is the hydraulic radius and k 0 is a constant parameter for a given pore shape. When the perforations are assumed of cylindrical shape with constant cross section, the viscous characteristic length and the thermal characteristic length are equal to the hydraulic radius of the perforation [START_REF] Allard | Propagation of Sound in Porous Media[END_REF][START_REF] Atalla | Modeling of perforated plates and screens using rigid frame porous models[END_REF][START_REF] Chevillote | Controlling sound absorption by an upstream resistive layer[END_REF][START_REF] Jaouen | Acoustical characterization of perforated facings[END_REF] Λ=Λ =h r =d/2. However as will be seen in the next section, this assumption cannot be made for other pore geometries.

Considering woven membranes as materials with identical uniform pores, the effective density can be linked to the dynamic bulk modulus given by 7,9

K(ω) = γP 0 γ -(γ -1)(1 + H 1 + jω/H/2 jω) , ( 6
)
where γ is the fluid specific heat ratio, P 0 is the fluid equilibrium pressure, ω the angular frequency and H = 16µ/(N Pr Λ 2 ρ 0 ). N Pr is the Prandtl number. The thermal length Λ characterizes the high frequency behavior of the bulk modulus [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF] . A simple model to calculate Λ is presented in section III.A.

Equations ( 3) and ( 6) allow to calculate the characteristic impedance Zc of the fluid,

Z c = K(ω)/ ρe (ω), (7)
and the associated wave number

k = ω ρe (ω)/K(ω). (8)
A. Simple model for calculating the thermal length of woven meshes

The thermal characteristic length has been defined as the surface to pore volume ration of the pore solid interface [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF] 2

Λ = S dS V dV = S V . ( 9
)
The volume of air passing through the pore depends on the end termination of the yarn geometry.

Assuming isotropic meshes with cylindrical yarns, an original model to obtain the thermal length is proposed. The elementary cell of the screen perforation is shown in Figure 1.

A simple calculation can be made by observing a single perforation as the one in Figure 2a.

The volume of air through the cell is observed in Figure 2b.

The volume of each truncated cylinder in the elementary cell of Figure 2a is

V h = πr 2 y b -4 ry -r y √ r 2 y -x 2 0 ydxdy, (10)
where r y is the yarn radius and b = d + 2r y . Solving the integral of the second term in the right hand side of Eq. (10) gives

V h = πr 2 y b -4 2r 3 y 3 . ( 11
)
The volume V p of the pore (Figure 2b) is then

V p = b 2 Y d - 4V h 2 . ( 12
)
Figure 2b represents Eq.( 12). The surface of each truncated cylinder is

S h = 2πr y b -4 π 0 r 2 y sinθdθ. ( 13
)
Solving the integral in Eq.( 13) we obtain

S h = 2πr y b -8r 2 y . (14)
The interior surface S p of the pore is then The thermal length can be obtained replacing S p in S and V p in V in Eq.( 9).

S p = 4S h 2 . ( 15
Λ = b 2 -bπr y + (8/3)r 2 y πb -4r y . (16)
In the examples studied, it is found that the thermal characteristic length calculated via Eq.( 16) can be smaller than the viscous characteristic length. This result tends to show that for the microstructure studied, the thermal effect can have a large influence on the bulk modulus at high frequencies.

Assuming identical pores for the particular geometry of a square pore of woven meshes, the physical parameters are: Λ = r, α ∞ = 1, σ is obtained with Eq.( 5) where [START_REF] Stinson | Propagation of sound and the assignment of shape factores in model porous materials having simple pore geometries[END_REF] k 0 = 1.78, Λ is obtained with Eq.( 16) and

φ = (d/(Y d + d)) 2 .

IV. THE HYBRID MODEL AND THE TRANSFER MATRIX APPROACH.

In the work of Pfretzschner et al., the impedance of a system composed of a perforated plate carrying a micrometric mesh was obtained using the impedance transfer method [START_REF] Pfretzschner | Microperforated Insertion Units: An alternative strategy to design microperforated panels[END_REF][START_REF] Cobo | Comparison of models describing double layer microperforated absorbers[END_REF] . Kang and Fuchs 11 used the electrical circuit approach for the same purpose. In the present work, the transfer matrix approach [START_REF] Allard | Propagation of Sound in Porous Media[END_REF] is used and the acoustic impedance of the whole system is obtained considering the impedance matrix of each element.A hybrid model combining the Maa and the JCA model to obtain the input impedance of a system composed of a macro-perforated plate and a micrometric commercial woven mesh is proposed.

The absorption curve of the hybrid model Maa-JCA depends on 7 parameters

(d 1 , t 1 , φ 1 , d 2 , t 2 , φ 2 , D)
, where d 1 , t 1 , φ 1 , d 2 , t 2 , φ 2 are the pore diameter, thickness and porosity of the carrying plate and the mesh respectively, and D is the air cavity thickness. Using Eq.( 1), the impedance matrix for a carrying plate with circular or slits perforations can be obtained as

T plate =            1 j(imag(Z 1 )) 0 1            . ( 17
)
In this study, carrying plates with perforations sizes greater than 1 mm and perforations rates close to 10% are considered, the viscous characteristic frequency is very low, so the thermoviscous effect can be disregarded 7 ; however, the mass effect is very high. As a consequence the impedance expression could be reduced in such a way that only the imaginary part is taken into account. Indeed, the experimental sound absorption coefficient of this plate (MPP or MSP without screen) coupled with an air cavity and rigid wall is weak.

The matrix for a mesh in the equivalent fluid context is given by

T Mesh =            cos(k M t M ) jZ cM sin(k M t M )/φ M φ M Z cM sin(k M t M ) cos(k M t M )            , ( 18
)
where t M ,φ M ,Z cM ,k M are the thickness, the porosity, the characteristic impedance and the wave number of the fluid in the mesh, respectively Z cM ,k M can be obtained with Eq.( 7) and ( 8) respectively.

For the air cavity, t M ,φ M ,Z cM ,k M are replaced by the properties of air i.e. D, 1, Z cA , k A which represent the thickness, the porosity, the characteristic impedance and wave number of the air in the cavity. The corresponding matrix is

T Air =            cos(k A t A ) jZ cA sin(k A t A )/φ A φ A Z cA sin(k A t A ) cos(k A t A )            . ( 19
)
A transfer matrix can be adapted for each new layer of plate, mesh or air. The multilayer configuration is obtained with Eq.( 24).

The schematic represented in Figure 3 shows the effect of the flow inside the perforations. The particle velocity is the same at both sides of the hybrid system. As φ 1 and φ 2 are always smaller than 1, the particle velocity inside the pores is always greater than outside them. The relationship between the pressure gradient ∆P and the particle velocity is described by Darcy's law as: FIG. 3. (color online)Viscous flow of the fluid through the perforations. Downstream the perforation, the "apparent" fluid velocity in the mesh is increased because of the squeezing due to the presence of the perforations upstream. This results in an "effective flow resistivity" that accounts for this diaphragm effect.

∆ P = -φϑσ. ( 20
)
This equation can be used to describe the flow in the woven mesh. However, when a macroperforated plate is coupled with a mesh without airspace between the two elements, a correction should be applied to this law to incorporate a "diaphragm effect" (Fig. 3). As the flow is squeezed inside the plate macroperforations, the flow velocity at the entrance of the pores of the mesh is increased and the flow resistivity of the mesh should be modified to account for this phenomenon.

The correction factor on the flow resistivity is precisely the porosity of the carrying plate. As a result, the resistivity of the mesh should be replaced by an "effective resistivity" of the combination plate-mesh

σ e f f = σ mesh σ plate . ( 21
)
Therefore the modified Darcy law for the mesh is

∆P = -φ mesh ϑ mesh σ mesh σ plate , (22)
and

∆P = -φ mesh ϑ mesh σ e f f . ( 23
)
The multiplication of the individual matrices results in the global transfer matrix of the hybrid model:

T h = [T Plate ][T Mesh ][T ]. ( 24
)
The input impedance of the system and absorption coefficient are then obtained by

Z h = T h 11 T h 21 , (25)
α = 1 - Z h -Z 0 Z h + Z 0 2 . ( 26
)
The effect of modifying the mesh resistivity by the effective resistivity is shown in Figure 4.

V. EXPERIMENTAL RESULTS

The experiments were carried out in a laboratory where the atmospheric conditions were:

T=23 • C, P o = 103070 Pa, c = 344.9 m/s, ρ 0 = 1.21 kg/m 3 , µ = 1.83 × 10 -5 Kg m -1 s -1 , γ = 1.4
and N pr = 0.71. Figure 5 shows the absorption curve of a commercial woven mesh placed in front of a rigid wall (d, t, φ, Y d , D)=(36 µm,50 µm,28 %,33 µm,5 cm), where d is the diameter of perforations, t the thickness, φ the perforation ratio, Y d the yarn diameter and D the air cavity thickness. A good agreement between the JCA model and the measurement is achieved. For this example, the manufacturer's nominal parameters where used and good results were obtained.

The real and imaginary parts of the dynamic tortuosity and bulk modulus of air in the mesh are observed in Figure 6. The experimental validation of the hybrid model is now proposed. Figure with a digital caliper. The porosity was also evaluated by the expression φ = (d/(Y d + d)) [START_REF] Rayleigh | The Theory of sound[END_REF] . The precision on the evaluation of the porosity when compared to the manufacturer's value was 5%.

Results are shown in Table 2. It is seen that when Y d ≥ 3.3d, the thermal length calculated by Eq. ( 16) is smaller than the viscous length taken to be equal to half the interfiber distance d. This is the case for all the meshes studied except for mesh a. woven mesh used for filtering, the second one is a simple woven curtain. The difference in costs between the two meshes is about 10 times and the results are very similar.

Figure 10 seems to validate the model also for a carrying plate with circular perforations. Mesh e is then combined with plates 2 and 4 (P2-HF and P4-HF, see Figure 11). It appears that using the same mesh with two different plates with the same perforation ratio but different perforations distribution provides very similar absorption curve. Deviations from the two curves in Figure 11 are due to the fact that the attached mass in plate 2 is slightly bigger than in plate 4. A simple equation for the calculation of the thermal characteristic length with square pores was used. This equation resulted from a simple modeling of the woven mesh elementary cell. It is found that in the woven mesh geometries studied, the thermal characteristic length can be smaller than the viscous characteristic length. This result tends to show that the thermal effects have a great importance in woven meshes at high frequencies.

A correction of the air flow resistivity of the mesh when placed behind a carrying plate was proposed. This correction accounts for a "diaphragm effect" due to the squeezing of the air flow inside the macroperforation, which results in an effective flow resistivity for the mesh. The transfer matrix method allows to handle multilayer configurations. The macroperforated plate can be studied in the framework of the classical Maa model while the porous layer can be modeled as an equivalent fluid with the JCA approach.

The costs of an absorbing system can be considerably reduced by using macro perforated panels and commercial meshes. Panels with perforations of the millimetric order are more manageable and can be drilled easily or even purchased. Micrometric meshes used for filtering are commercially manufactured and even woven curtains can be used for this purpose.

  FIG. 2. (color online)Simple model for a single mesh perforation geometry. a) fiber structure (the dotted rectangle represents the elementary cell), b) inside volume.

FIG. 4 .

 4 FIG.4. Absorption curve of a carrying plate with (d 1 , t 1 , φ 1 ) = (3 mm, 1 mm, 10 %) with a micrometric mesh (d 2 , t 2 , φ 2 )=(50 µm, 50 µm, 37 %) and D=3.5 cm. (-) with correction for resistivity, (--) without correction.

Figure 8 FIG. 5 .

 85 Figure8shows the five meshes used. The pores and yarn diameter were obtained using the software Optika vision lite and an Euromex Holland microscope. The microscopic analysis allowed to calculate the thermal and viscous characteristic lengths. The thickness was measured

Figure 9

 9 Figure 9 illustrates the results of combining plate 1 with meshes a, b and d. Two interesting results are highlighted. First, the hybrid model presents a very good agreement with the experimental measurements for precision woven and non-precision woven meshes. Second, a very high absorption is achieved by using mesh b and d and while the first one is a commercial precision

FIG. 6 .

 6 FIG. 6. (color online)Dynamic tortuosity and bulk modulus of air for a precision woven micrometric mesh Sefar Nitex 36/28. a) and b) real parts, c) and d) imaginary parts.

FIG. 7 .FIG. 9 .

 79 FIG. 7. (color online)MPPs and MSPs with different perforations distribution.

FIG. 10 .

 10 FIG. 10. (color online)Absorption curve of plate 3 combined with mesh c. Air cavity D=2.4 cm.

FIG. 11 .

 11 FIG. 11. (color online)Absorption curve of plate 2(+) and plate 4(o) combined with mesh e. Air cavity D=3.13 cm.

FIG. 1 ( 9 FIG. 2 ( 10 FIG. 3 ( 12 FIG. 4 14 FIG. 5 ( 15 FIG. 6 ( 16 FIG. 7 (

 192103124145156167 FIG. 1 (color online)Microscopic image of an elementary cell of a woven mesh. . . . . 9 FIG. 2 (color online)Simple model for a single mesh perforation geometry. a) fiber structure (the dotted rectangle represents the elementary cell), b) inside volume. . . . 10 FIG. 3 (color online)Viscous flow of the fluid through the perforations. Downstream the perforation, the "apparent" fluid velocity in the mesh is increased because of the squeezing due to the presence of the perforations upstream. This results in an "effective flow resistivity" that accounts for this diaphragm effect. . . . . . . . . 12 FIG. 4 Absorption curve of a carrying plate with (d 1 , t 1 , φ 1 ) = (3 mm, 1 mm, 10 %) with a micrometric mesh (d 2 , t 2 , φ 2 )=(50 µm, 50 µm, 37 %) and D=3.5 cm. (-) with correction for resistivity, (--) without correction. . . . . . . . . . . . . . . . . . 14 FIG. 5 (color online)Absorption curve of a mesh with parameters (d, t, φ, Y d , D)=(36 µm,50 µm,28 %,33 µm,5 cm) . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 FIG. 6 (color online)Dynamic tortuosity and bulk modulus of air for a precision woven micrometric mesh Sefar Nitex 36/28. a) and b) real parts, c) and d) imaginary parts. 16 FIG. 7 (color online)MPPs and MSPs with different perforations distribution. . . . . . . 17

FIG. 8 (

 8 FIG. 8 (color online)Microscopic image of the meshes tested a)Woven screen, b) sefar R acoustic BHY, c) sefar R Nitex 41/31, d) curtain e) sefar acoustic HF 34-39. 17

FIG. 9 ( 18 FIG. 10 ( 19 FIG. 11 (

 918101911 FIG. 9 (color online)Absorption curve of plate 1 combined with mesh a (+), b (o), and d (x). Air cavity D=2 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 FIG. 10 (color online)Absorption curve of plate 3 combined with mesh c. Air cavity D=2.4 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 FIG. 11 (color online)Absorption curve of plate 2(+) and plate 4(o) combined with mesh e. Air cavity D=3.13 cm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

  (µm) (µm) (%) (µm) (µm) (µm) (N sm -4 ) (N sm -4 ) a 59.18 110 5.29 198.21 29.59 29.84 2.82 × 10 6 2.83 × 10 7 b 39 49 34.9 33 19.50 13.59 9.84 × 10 5 9.88 × 10 6 c 45.7 45.87 36.35 30.1 22.85 15.59 6.85 × 10 5 6.92 × 10 6 d 94.25 120 19.7 117.8 47.12 34.78 2.97 × 10 5 2.99 × 10 6 e 39 45.37 34.90 30.29 19.50 13.46 9.84 × 10 5 4 9.88 × 10 6

		TABLE II. Characteristics of the meshes
		d		t	φ	Y d	Λ	Λ	σ	σ e f f
		Mesh
	Slit	F (e)	y	2 √ -j tanh(y J 1 J 0 √ √ -j √ -j 1-1.4092(ξ)+0.33818(ξ) 3 +0.06793(ξ) 5 -0.02287(ξ) 6 +0.03015(ξ) 7 -0.01641(ξ) 8 j y √ j

ACKNOWLEDGMENTS

This work has been supported by the Spanish Ministry of Science and Innovation (MICINN)

through Projects TRA2008-05654-C03-03 and TRA2011-26261-C04-01.