Guilhem Castagnos
email: guilhem.castagnos@math.u-bordeaux1.fr

Antoine Joux
email: antoine.joux@m4x.org

Fabien Laguillaumie

Phong Q Nguyen

Factoring pq 2 with Quadratic Forms: Nice Cryptanalyses

Keywords: Public-key Cryptanalysis, Factorisation, Binary Quadratic Forms, Homogeneous Coppersmith's Root Finding, Lattices

We present a new algorithm based on binary quadratic forms to factor integers of the form N = pq 2 . Its heuristic running time is exponential in the general case, but becomes polynomial when special (arithmetic) hints are available, which is exactly the case for the so-called NICE family of public-key cryptosystems based on quadratic fields introduced in the late 90s. Such cryptosystems come in two flavours, depending on whether the quadratic field is imaginary or real. Our factoring algorithm yields a general key-recovery polynomial-time attack on NICE, which works for both versions: Castagnos and Laguillaumie recently obtained a total break of imaginary-NICE, but their attack could not apply to real-NICE. Our algorithm is rather different from classical factoring algorithms: it combines Lagrange's reduction of quadratic forms with a provable variant of Coppersmith's lattice-based root finding algorithm for homogeneous polynomials. It is very efficient given either of the following arithmetic hints: the public key of imaginary-NICE, which provides an alternative to the CL attack; or the knowledge that the regulator of the quadratic field Q( √ p) is unusually small, just like in real-NICE.

Introduction

Many public-key cryptosystems require the hardness of factoring large integers of the special form N = pq 2 , such as Okamoto's Esign [START_REF] Okamoto | A fast signature scheme based on congruential polynomial operations[END_REF], Okamoto and Uchiyama's encryption [START_REF] Okamoto | A New Public-Key Cryptosystem as Secure as Factoring[END_REF], Takagi's fast RSA variants [START_REF] Takagi | Fast RSA-Type Cryptosystem Modulo p k q[END_REF], and the large family (surveyed in [START_REF] Buchmann | Number Field Cryptography[END_REF]) of cryptosystems based on quadratic fields, which was initiated by Buchmann and Williams' key exchange [START_REF] Buchmann | A Key-Exchange System based on Imaginary Quadratic Fields[END_REF], and which includes NICE 1 cryptosystems [HPT99,PT99,PT00,JSW08] (whose main feature is a quadratic decryption). These moduli are popular because they can lead to special functionalities (like homomorphic encryption) or improved efficiency (compared to RSA). And no significant weakness has been found compared to standard RSA moduli of the form N = pq: to the best of our knowledge, the only results on pq 2 factorisation are [START_REF] Peralta | Faster Factoring of Integers of a Special Form[END_REF][START_REF] Peralta | Elliptic curve factorization using a "partially oblivious" function. Cryptography and computational number theory[END_REF][START_REF] Boneh | Factoring N = p r q for large r[END_REF]. More precisely, [START_REF] Peralta | Faster Factoring of Integers of a Special Form[END_REF][START_REF] Peralta | Elliptic curve factorization using a "partially oblivious" function. Cryptography and computational number theory[END_REF] obtained a linear speed-up of Lenstra's ECM, and [BDH99, Sect. 6] can factor in time Õ(N 1/9 ) when p and q are balanced. Furthermore, computing the "squarefree part" of an integer (that is, given N ∈ N as input, compute (r, s) ∈ N 2 such that N = r 2 s with s squarefree) is a classical problem in algorithmic number theory (cf. [START_REF] Adleman | Open problems in number theoretic complexity[END_REF]), because it is polynomial-time equivalent to determining the ring of integers of a number field [START_REF] Chistov | The complexity of constructing the ring of integers of a global field[END_REF].

However, some of these cryptosystems actually provide additional information (other than N ) in the public key, which may render factorisation easy. For instance, Howgrave-Graham [START_REF] Howgrave-Graham | Approximate Integer Common Divisors Proc[END_REF] showed that the public key of [START_REF] Okamoto | Fast public-key cryptosystem using congruent polynomial equations[END_REF] disclosed the secret factorisation in polynomial time, using the gcd extension of Coppersmith's root finding method [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF]. Very recently, Castagnos and Laguillaumie [START_REF] Castagnos | On the Security of Cryptosystems with Quadratic Decryption: The Nicest Cryptanalysis[END_REF] showed that the public key in the imaginary version [HPT99, PT99, PT00] of NICE allowed to retrieve the secret factorisation in polynomial time. And this additional information in the public key was crucial to make the complexity of decryption quadratic in imaginary-NICE, which was the main claimed benefit of NICE. But surprisingly, the attack of [START_REF] Castagnos | On the Security of Cryptosystems with Quadratic Decryption: The Nicest Cryptanalysis[END_REF] does not work against REAL-NICE [JSW08], which is the version of NICE with real (rather than imaginary) quadratic fields, and which also offers quadratic decryption. In particular, the public key of REAL-NICE only consists of N = pq 2 , but the prime p has special arithmetic properties. Our Results. We present a new algorithm to factor integers of the form N = pq 2 , based on binary quadratic forms (or equivalently, ideals of orders of quadratic number fields). In the worst case, its heuristic running time is exponential, namely Õ(p 1/2 ). But in the presence of special hints, it becomes heuristically polynomial. These hints are different from the usual ones of lattice-based factoring methods [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF][START_REF] Boneh | Factoring N = p r q for large r[END_REF][START_REF] Howgrave-Graham | Approximate Integer Common Divisors Proc[END_REF] where they are a fraction of the bits of the secret prime factors. Instead, our hints are arithmetic, and correspond exactly to the situation of NICE, including both the imaginary [HPT99, PT99, PT00] and real versions [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]. This gives rise to the first general key-recovery polynomial-time attack on NICE, using only the public key.

More precisely, our arithmetic hints can be either of the following two:

i. The hint is an ideal equivalent to a secret ideal of norm q 2 in an imaginary quadratic field of discriminant -pq 2 : in NICE, such an ideal is disclosed by the public key. This gives an alternative attack of NICE, different from [START_REF] Castagnos | On the Security of Cryptosystems with Quadratic Decryption: The Nicest Cryptanalysis[END_REF]. ii. The hint is the knowledge that the regulator of the quadratic field Q( √ p) is unusually small, just like in REAL-NICE. Roughly speaking, the regulator is a

1 for New Ideal Coset Encryption real number which determines how "dense" the units of the ring of integers of the number field Q( √ p) are. This number is known to lie in the large interval log 1 2 ( √ p -4 + √ p) , 1 2 p 1 2 log p + 1 . But for infinitely many p (including square-free numbers of the form p = k 2 + r, where p > 5, r|4k and -k < r ≤ k, see [START_REF] Degert | Uber die Bestimmung der Grundeinheit gewisser reell-quadratischer Zhalkörper[END_REF]), the regulator is at most polynomial in log p. For these unusually small regulators, our algorithm heuristically runs in time polynomial in the bit-length of N = pq 2 , which gives the first total break of REAL-NICE [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]. We stress that although such p's are easy to construct, their density is believed to be arbitrary small. Interestingly, our algorithm is rather different from classical factoring algorithms. It is a combination of Lagrange's reduction of quadratic forms with a provable variant of Coppersmith's lattice-based root finding algorithm [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF] for homogeneous polynomials. In a nutshell, our factoring method first looks for a reduced binary quadratic form f (x, y) = ax 2 + bxy + cy 2 representing properly q 2 with small coefficients, i.e. there exist small coprime integers x 0 and y 0 such that q 2 = f (x 0 , y 0 ). In case i., such a quadratic form is already given. In case ii., such a quadratic form is found by a walk along the principal cycle of the class group of discriminant pq 2 , using Lagrange's reduction of (indefinite) quadratic forms. Finally, the algorithm finds such small coprime integers x 0 and y 0 such that q 2 = f (x 0 , y 0 ), by using the fact that gcd(f (x 0 , y 0 ), pq 2 ) is large. This discloses q 2 and therefore the factorisation of N . In both cases, the search for x 0 and y 0 is done with a new rigorous homogeneous bivariate variant of Coppersmith's method, which might be of independent interest: by the way, it was pointed out to us that Bernstein [START_REF] Bernstein | List decoding for binary Goppa codes[END_REF] independently used a similar method in the different context of Goppa codes decoding.

Our algorithm requires "natural" bounds on the roots of reduced quadratic forms of a special shape. We are unable to prove rigorously these bounds, which makes our algorithm heuristic (like many factoring algorithms). But we have performed many experiments supporting such bounds, and the algorithm works very well in practice.

Factorisation and Quadratic Forms. Our algorithm is based on quadratic forms, which share a long history with factoring (see [START_REF] Crandall | Prime Numbers: A Computational Perspective[END_REF]). Fermat's factoring method represents N in two intrinsically different ways by the quadratic form x 2 + y 2 . It has been improved by Shanks with SQUFOF, whose complexity is Õ(N 1/4 ) (see [START_REF] Gower | Square form factorization[END_REF] for a detailed analysis). Like ours, this method works with the infrastructure of a class group of positive discriminant, but is different in spirit since it searches for an ambiguous form (after having found a square form), and does not focus on discriminants of a special shape. Schoof's factoring algorithms [START_REF] Schoof | Quadratic fields and factorization[END_REF] are also essentially looking for ambiguous forms. One is based on computation in class groups of complex quadratic orders and the other is close to SQUFOF since it works with real quadratic orders by computing a good approximation of the regulator to find an ambiguous form. Like SQUFOF, this algorithm does not takes advantage of working in a non-maximal order and is rather different from our algorithm. Both algorithms of [START_REF] Schoof | Quadratic fields and factorization[END_REF] runs in Õ(N 1/5 ) under the generalised Riemann hypothesis. McKee's method [START_REF] Mckee | Speeding Fermat's factoring method[END_REF] is a speedup of Fermat's algorithm (and was presented as an alternative to SQUFOF) with a heuristic complexity of Õ(N 1/4 ) instead of Õ(N 1/2 ).

SQUFOF and other exponential methods are often used to factor small numbers (say 50 to 100 bits), for instance in the post-sieving phase of the Number Field Sieve algorithm. Some interesting experimental comparisons can be found in [START_REF] Milan | Factoring Small Integers: An Experimental Comparison[END_REF]. Note that the currently fastest rigorous deterministic algorithm actually has exponential complexity: it is based on a polynomial evaluation method (for a polynomial of the form x(x -1) • • • (x -B + 1) for some bound B) and its best variant is described in [START_REF] Bostan | Linear Recurrences with Polynomial Coefficients and Application to Integer Factorization and Cartier-Manin Operator[END_REF]. Finally, all sieve factoring algorithms are somewhat related to quadratic forms, since their goal is to find random pairs (x, y) of integers such that x 2 ≡ y 2 mod N . However, these algorithms factor generic numbers and have a subexponential complexity. Road Map. The rest of the paper is organised as follows. The first section recalls facts on quadratic fields and quadratic forms, and present our heuristic supported by experiments. The next section describes the homogeneous Coppersmith method and the following exhibits our main result: the factoring algorithm. The last section consists of the two cryptanalyses of cryptosystems based on real quadratic fields (REAL-NICE) and on imaginary quadratic fields (NICE).

2 Background on Quadratic Fields and Quadratic Forms

Quadratic Fields

Let D = 0, 1 be a squarefree integer and consider the quadratic number field K = Q( √ D). If D < 0 (resp. D > 0), K is called an imaginary (resp. a real ) quadratic field. The fundamental discriminant ∆ K of K is defined as ∆ K = D if D ≡ 1 (mod 4) and ∆ K = 4D otherwise. An order O in K is a subset of K such that O is a subring of K containing 1 and O is a free Z-module of rank 2. The ring O ∆ K of algebraic integers in K is the maximal order of K. It can be written as Z + ω K Z, where

ω K = 1 2 (∆ K + √ ∆ K ). If we set f = [O ∆ K : O] the finite index of any order O in O ∆ K , then O = Z + f ω K Z. The integer f is called the conductor of O. The discriminant of O is then ∆ f = f 2 ∆ K . Now, let O ∆ be an order of discriminant ∆ and a be a nonzero ideal of O ∆ , its norm is N(a) = |O ∆ /a|. A fractional ideal is a subset a ⊂ K such that da is an ideal of O ∆ for d ∈ N. A fractional ideal a is said to be invertible if there exists an another fractional ideal b such that ab = O ∆ . The ideal class group of O ∆ is C(O ∆ ) = I(O ∆ )/P (O ∆ )
, where I(O ∆ ) is the group of invertible fractional ideals of O ∆ and P (O ∆ ) the subgroup consisting of principal ideals. Its cardinality is the class number of O ∆ denoted by h(O ∆ ). A nonzero ideal a of O ∆ , a is said to be prime to f if a + f O ∆ = O ∆ . We denote by I(O ∆ , f ) the subgroup of I(O ∆ ) of ideals prime to f . The group O ∆ of units in O ∆ is equal to {±1} for all ∆ < 0, except when ∆ is equal to -3 and -4 (O -3 and O -4 are respectively the group of sixth and fourth roots of unity). When ∆ > 0, then

O ∆ = -1, ε ∆ where ε ∆ > 0 is called the fundamental unit. The real number R ∆ = log(ε ∆ ) is the regulator of O ∆ .
The following important bounds on the regulator of a real quadratic field can be found in [START_REF] Jacobson | An investigation of bounds for the regulator of quadratic fields[END_REF]:

log 1 2 ( √ ∆ -4 + √ ∆) ≤ R ∆ < 1 2 ∆ 1 2 log ∆ + 1 . (1)
The lower bound is reached infinitely often, for instance with ∆ = x 2 + 4 with 2 x. Finally, this last proposition is the heart of both NICE and REAL-NICE.

Proposition 1 ([Cox99, Proposition 7.20] [Wei04, Theorem 2.16]). Let O ∆ f be an order of conductor f in a quadratic field K. i. If A is an O ∆ K -ideal prime to f , then A ∩ O ∆ f is an O ∆ f -ideal prime to f of the same norm. ii. If a is an O ∆ f -ideal prime to f , then aO ∆ K is an O ∆ K -ideal prime to f of the same norm. iii. The map ϕ f : I(O ∆ f , f ) → I(O ∆ K , f ), a → aO ∆ K is an isomorphism. The map ϕ f from Proposition 1 induces a surjection φf : C(O ∆ f ) C(O ∆ K )
which can be efficiently computed (see [START_REF] Paulus | A New Public-Key Cryptosystem over a Quadratic Order with Quadratic Decryption Time[END_REF]). In our settings, we will use a prime conductor f = q and consider ∆ q = q 2 ∆ K , for a fundamental discriminant ∆ K . In that case, the order of the kernel of φq is given by the classical analytic class number formula (see for instance [START_REF] Buchmann | Binary Quadratic Forms[END_REF])

h(O ∆q ) h(O ∆ K ) = q -(∆ K /q) if ∆ k < -4, (q -(∆ K /q))R ∆ K /R ∆q if ∆ k > 0.
(2)

Note that in the case of real quadratic fields, ∆q = t ∆ K for a positive integer t, hence R ∆q /R ∆ K = t and t | (q -(∆ K /q)).

Representation of the Classes

Working with ideals modulo the equivalence relation of the class group is essentially equivalent to work with binary quadratic forms modulo SL 2 (Z) (cf. Section 5.2 of [START_REF] Cohen | A Course in Computational Algebraic Number Theory[END_REF]). Moreover, quadratic forms are more suited to an algorithmic point of view. Every ideal a of O ∆ can be written as a = m aZ

+ -b+ √ ∆ 2 Z with m ∈ Z, a ∈ N and b ∈ Z such that b 2 ≡ ∆ (mod 4a
). In the remainder, we will only consider primitive integral ideals, which are those with m = 1. This notation also represents the binary quadratic form ax 2 + bxy + cy 2 with b 2 -4ac = ∆. This representation of the ideal is unique if the form is normal (see below). We recall here some facts about binary quadratic forms. Definition 1. A binary quadratic form f is a degree 2 homogeneous polynomial f (x, y) = ax 2 +bxy +cy 2 where a, b and c are integers, and is denoted by

[a, b, c]. The discriminant of the form is ∆ = b 2 -4ac. If a > 0 and ∆ < 0, the form is called definite positive and indefinite if ∆ > 0. Let M ∈ SL 2 (Z) with M = α β γ δ
, and f = [a, b, c], a binary quadratic form, then f.M is the equivalent binary quadratic form f (αx + βy, γx + δy).

Definite Positive Forms. Let us first define the crucial notion of reduction.

Definition 2. The form f = [a, b, c] is called normal if -a < b ≤ a. It is called reduced if it is normal, a ≤ c, and if b ≥ 0 for a = c.
The procedure which transforms a form f = [a, b, c] into a normal one consists in setting s such that b + 2sa belongs to the right interval (see [START_REF] Buchmann | Binary Quadratic Forms[END_REF](5.4)]) and producing the form [a, b+2sa, as

2 +bs+c]. Once a form f = [a, b, c] is normalised, a reduction step consists in normalising the form [c, -b, a].
We denote this form by ρ(f ) and by Rho a corresponding algorithm. The reduction then consists in normalising f , and then iteratively replacing f by ρ(f ) until f is reduced. The time complexity of this (Lagrange-Gauß) algorithm is quadratic (see [START_REF] Vallée | Lattice Reduction in Two Dimensions: Analyses under Realistic Probabilistic Models[END_REF]). It returns a reduced form g which is equivalent to f modulo SL 2 (Z). We will call matrix of the reduction, the matrix M such that g = f.M . The reduction procedure yields a uniquely determined reduced form in the class modulo SL 2 (Z). Indefinite Forms. Our main result will deal with forms of positive discriminant. Here is the definition of a reduced indefinite form.

Definition 3. The form f = [a, b, c] of positive discriminant ∆ is reduced if √ ∆ -2|a| < b < √ ∆ and normal if -|a| < b ≤ |a| for |a| ≥ √ ∆, and √ ∆ -2|a| < b < √ ∆ for |a| < √ ∆.
The reduction process is similar to the definite positive case. The time complexity of the algorithm is still quadratic (see [BV07, Theorem 6.6.4]). It returns a reduced form g which is equivalent to f modulo SL 2 (Z). The main difference with forms of negative discriminant is that there will in general not exist a unique reduced form per class, but several organised in a cycle structure i. e., when f has been reduced then subsequent applications of ρ give other reduced forms. Definition 4. Let f be an indefinite binary quadratic form, the cycle of f is the sequence (ρ i (g)) i∈Z where g is a reduced form which is equivalent to f . From Theorem 6.10.3 from [START_REF] Buchmann | Binary Quadratic Forms[END_REF], the cycle of f consists of all reduced forms in the equivalence class of f . Actually, the complete cycle is obtained by a finite number of application of ρ as the process is periodic. It has been shown in [START_REF] Buchmann | Short Representation of Quadratic Integers[END_REF] that the period length of the sequence of reduced forms in each class of a class group of discriminant ∆ satisfies R ∆ log ∆ ≤ ≤ 2R ∆ log 2 + 1. Our factoring algorithm will actually take place in the principal equivalence class. The following definition exhibits the principal form of discriminant ∆.

Definition 5. The reduced form [1, √ ∆ , ( √ ∆ 2 -∆)/4
] of discriminant ∆ is called the principal form of discriminant ∆, and will be denoted 1 ∆ .

2.3 Reduction of the Forms [q 2 , kq, (k 2 ± p)/4] and Heuristics

In this subsection, p and q are two distinct primes of the same bit-size λ and p ≡ 1 mod 4 (resp. p ≡ 3 mod 4) when we deal with positive (resp. negative) discriminant. Our goal is to factor the numbers pq 2 with the special normalised quadratic forms [q 2 , kq, (k 2 + p)/4] or [q 2 , kq, (k 2 -p)/4], depending whether we work with a negative discriminant ∆ q = -pq 2 or with a positive one ∆ q = pq 2 . If p and q have the same size, these forms are clearly not reduced neither in the imaginary setting nor in real one. But as we shall see, we can find the reduced forms which correspond to the output of the reduction algorithm applied on these forms.

Suppose that we know a form fk , either definite positive or indefinite, which is the reduction of a form f k = [q 2 , kq, (k 2 ± p)/4] where k is an integer. Then fk represents the number q 2 . More precisely, if

M k = α β γ δ ∈ SL 2 (Z) is the matrix such that fk = f k .M k , then fk .M -1 k = f k and q 2 = f k (1, 0) = fk (δ, -γ).
In Section 3, we will see that provided they are relatively small compared to ∆ q , the values δ and -γ can be found in polynomial time with a new variant of Coppersmith method. Our factoring algorithm can be sketched as follows: find such a form fk and if the coefficients of M k are sufficiently small, retrieve δ and -γ and the non-trivial factor q 2 of ∆ q . In this paragraph, we give some heuristics on the size of such a matrix M k and discuss their relevance. If M is a matrix we denote by |M | the max norm, i. e., the maximal coefficient of M in absolute value.

In the imaginary case, it is showed in the proof of [CL09, Theorem 2] that the forms f k belong to different classes of the kernel of the map φq , depending on k, so the reduced equivalent forms fk are the unique reduced elements of the classes of the kernel. To prove the correctness of our attack on NICE, we need the following heuristic (indeed, the root finding algorithm of Section 3 recovers roots up to |∆ q | 1/9 ): Heuristic 1 (Imaginary case) Given a reduced element fk of a nontrivial class of ker φq , the matrix of reduction M k is such that |M k | < |∆ q | 1/9 with probability asymptotically close to 1.

From Lemma 5.6.1 of [BV07], |M k | < 2 max{q 2 , (k 2 + p)/4}/ pq 2 . As f k is normalised, |k| ≤ q and |M k | < 2q/ √ p ≈ |∆ q | 1/6
. Note that we cannot reach such a bound with our root finding algorithm. Experimentally, for random k, |M k | can be much smaller. For example, if the bit-size λ of p and q equals 100, the mean value of |M k | is around |∆ q | 1/11.7 . Our heuristic can be explained as follows.

A well-known heuristic in the reduction of positive definite quadratic forms (or equivalently, two-dimensional lattices) is that if [a, b, c] is a reduced quadratic form of discriminant ∆, then a and c should be close to √ ∆. This cannot hold for all reduced forms, but it can be proved to hold for an overwhelming majority of reduced forms. Applied to fk = [a, b, c], this means that we expect a and c to be close to |∆ q | 1/2 . Now, recall that q 2 = fk (δ, -γ) = aδ 2 -bδγ + cγ 2 , which leads to δ and γ close to q 2 /a = q/ √ a ≈ q/|∆ q | 1/4 ≈ |∆ q | 1/12 . Thus, we expect that |M k | ≤ |∆ q | 1/12 . And this explains why we obtained experimentally the bound |∆ q | 1/11.7 . Figure 1(a) shows a curve obtained by experimentation, which gives the probability that |M k | < |∆ q | 1/9 for random k, in function of λ. This curve also supports our heuristic.

In the real case, we prove in the following theorem that R ∆q /R ∆ K forms f k are principal and we exhibit the generators of the corresponding primitive ideals.

Theorem 1. Let ∆ K be a fundamental positive discriminant, ∆ q = ∆ K q 2 where q is an odd prime conductor. Let ε ∆ K (resp. ε ∆q ) be the fundamental unit of O ∆ K (resp. O ∆q ) and t such that ε t ∆ K = ε ∆q . Then the principal ideals of O ∆q generated by qε i ∆ K correspond to quadratic forms f k(i) = [q 2 , k(i)q, (k(i) 2 -p)/4] with i ∈ {1, . . . , t -1} and k(i) is an integer defined modulo 2q computable from ε i ∆ K mod q.

Proof. Let α i = qε i ∆ K with i ∈ {1, . . . , t -1}. Following the proof of [BTW95, Proposition 2.9], we detail here the computation of a i = α i O ∆q . Let x i and y i be two integers such that

ε i ∆ K = x i + y i ω K . Then α i = qx i + y i q∆ K (1 -q)/2 + y i 1 2 (∆ q + ∆ q ),
and α i is an element of O ∆q . Let m i , a i and b i be three integers such that

a i = m i a i Z + -bi+ √ ∆q 2
. As mentioned in the proof of [BTW95, Proposition 2.9], m i is the smallest positive coefficient of ∆ q /2 in a i . As O ∆q is equal to Z + (∆ q + ∆ q )/2Z, α i O ∆q is generated by α i and α i (∆ q + ∆ q )/2 as a Z-module. So a simple calculation gives that m i = gcd(y i , q(x i + y i ∆ K /2)). As ε i ∆ K is not an element of O ∆q , we have gcd(y i , q) = 1 so m i = gcd(y i , x i + y i ∆ K /2). The same calculation to find m i for the ideal

ε i ∆ K O ∆ K reveals that m i = m i . As ε i ∆ K O ∆ K = O ∆ K we must have m i = 1. Now, N(a i ) = |N(α i )| = q 2 and N(a i ) = m 2
i a i = a i and therefore a i = q 2 . Let us now find b i . Note that b i is defined modulo 2a i . Since α i ∈ α i O ∆q , there exist µ i and ν i such that α i = a i µ i + (-b i + ∆ q )/2ν i . By identification in the basis (1, ∆ q ), ν k = 1 and by a multiplication by 2, we obtain 2qx i + qy i ∆ K ≡ -b i y i (mod 2a i ). As b i ≡ ∆ q (mod 2), we only have to determine b i modulo q 2 . As y i is prime to q, we have b i ≡ k(i)q (mod q 2 ) with k(i) ≡ -2x i /y i -∆ K (mod q). Finally, as we must have -a i < b ≤ a i if a i > ∆ q and else ∆ q -2a i < b < ∆ q , k(i) is the unique integer with k(i) ≡ ∆ q (mod 2) and k(i) ≡ -2x i /y i -∆ K (mod q), such that b = k(i)q satisfies that inequalities. Eventually, the principal ideal of O ∆q generated by

qε i ∆ K corresponds to the form [q 2 , k(i)q, c i ] with c i = (b 2 i -∆ q )/(4a i ) = (k(i) 2 -∆ K )/4.
From this theorem, we see that if we go across the cycle of principal forms, then we will find reduced forms fk . To analyse the complexity of our factoring algorithm, we have to know the distribution of these forms on the cycle. An appropriate tool is the Shanks distance d (see [BV07, Definition 10.1.4]) which is close to the number of iterations of Rho between two forms. One has all i = 1, 2, . . . , t -1. Let j be the smallest integer such that 0 < jR ∆ K -2 log q, then as jR

d(1 ∆q , f k(i) ) = iR ∆ K . From Lemma 10.1.8 of [BV07], |d( fk(i) , f k(i) )| < log q, for
∆ K = d(f k(i) , f k(i+j) ) = d(f k(i) , fk(i) ) + d( fk(i) , fk(i+j) ) + d( fk(i+j) , f k(i+j)
), from the triangle inequality, one has jR ∆ K < 2 log(q) + |d( fk(i) , fk(i+j) )|. So, |d( fk(i) , fk(i+j) )| > jR ∆ K -2 log q > 0. This inequality proves that f k(i) and f k(i+j) do not reduce to the same form. Experiments actually show that asymptotically, |d( fk(i) , f k(i) )| is very small on average (smaller than 1). As a consequence, as pictured in figure 2,d

(1 ∆q , fk(i) ) ≈ iR ∆ K . R∆ K f k(1) fk(1) f k(2) fk(2) f k(3) fk(3) 1∆ q

Fig. 2. Repartition of the forms fk(i) along the principal cycle

Moreover, as in the imaginary case, experiments show that asymptotically the probability that the norm of the matrices of reduction, |M k | is smaller than ∆ 1/9 q is close to 1 (see figure 1(b)). This leads to the following heuristic.

Heuristic 2 (Real case) From the principal form 1 ∆q , a reduced form fk such that the matrix of the reduction, M k , satisfy |M k | < ∆ 1/9 q , can be found in O(R ∆ K ) successive applications of Rho.

We did also some experiments to investigate the case where the bit-sizes of p and q are unbalanced. In particular when the size of q grows, the norm of the matrix of reduction becomes larger. For example, for a 100-bit p and a 200-bit q (resp. a 300-bit q), more than 95% (resp. 90%) of the fk have a matrix M k with

|M k | < ∆ 1/6.25 q (resp. |M k | < ∆ 1/5.44 q ).

A Rigorous Homogeneous Variant of Coppersmith's Root Finding Method

Our factoring algorithm searches many times for small modular roots of degree two homogeneous polynomials and the most popular technique to find them is based on Coppersmith's method (see [START_REF] Coppersmith | Small Solutions to Polynomial Equations, and Low Exponent RSA Vulnerabilities[END_REF] or May's survey [START_REF] May | Using LLL-Reduction for Solving RSA and Factorization Problems: A Survey[END_REF]).

Our problem is the following: Given f (x, y) = x 2 + bxy + cy 2 a (monic) binary quadratic form and N = pq 2 an integer of unknown factorisation, find (x 0 , y 0 ) ∈ Z 2 such that f (x 0 , y 0 ) ≡ 0 (mod q 2 ), while |x 0 |, |y 0 | ≤ M , where M ∈ N. The usual technique for this kind of problems is only heuristic, since it is the gcd extension of bivariate congruences. Moreover, precise bounds cannot be found in the litterature. Fortunately, because our polynomial is homogeneous, we will actually be able to prove the method. This homogenous variant is quite similar to the one-variable standard Coppersmith method, but is indeed even simpler to describe and more efficient since there is no need to balance coefficients. We denote as • the usual Euclidean norm for polynomials. The main tool to solve this problem is given by the following variant of the widespread elementary Howgrave-Graham's lemma [START_REF] Howgrave-Graham | Finding small roots of univariate modular equations revisited[END_REF].

Lemma 1. Let B ∈ N and g(x, y) ∈ Z[x, y] be a homogeneous polynomial of total degree δ. Let M > 0 be a real number and suppose that ||g(x, y)|| <

B √ δ+1M δ
then for all x 0 , y 0 ∈ Z such that g(x 0 , y 0 ) ≡ 0 (mod B) and |x 0 |, |y 0 | ≤ M , g(x 0 , y 0 ) = 0.

Proof. Let g(x, y) = δ i=0 g i x i y δ-i where some g i s might be zero. We have

|g(x 0 , y 0 )| ≤ δ i=0 |g i ||x i 0 y δ-i 0 | ≤ M δ δ i=0 |g i | ≤ M δ √ δ + 1 g(x, y) < B
and therefore g(x 0 , y 0 ) = 0.

The trick is then to find only one small enough bivariate homogeneous polynomial satisfying the conditions of this lemma and to extract the rational root of the corresponding univariate polynomial with standard techniques. On the contrary, the original Howgrave-Graham's lemma suggests to look for two polynomials of small norm having (x 0 , y 0 ) as integral root, and to recover it via elimination theory. The usual way to obtain these polynomials is to form a lattice spanned by a special family of polynomials, and to use the LLL algorithm (cf. [START_REF] Lenstra | Factoring Polynomials with Rational Coefficients[END_REF]) to obtain the two "small" polynomials. Unfortunately, this reduction does not guarantee that these polynomials will be algebraically independent, and the elimination can then lead to a trivial relation. Consequently, this bivariate approach is heuristic. Fortunately, for homogeneous polynomials, we can take another approach by using Lemma 1 and then considering a univariate polynomial with a rational root. This makes the method rigorous and slightly simpler since we need a bound on g(x, y) and not on g(xX, yY ) if X and Y are bounds on the roots and therefore the resulting lattice has smaller determinant than in the classical bivariate approach.

To evaluate the maximum of the bound we can obtain, we need the size of the first vector provided by LLL which is given by: Lemma 2 (LLL). Let L be a full-rank lattice in Z d spanned by an integer basis B = {b 1 , . . . , b d }. The LLL algorithm, given B as input, will output in time O(d 6 log 3 (max b i )) a non-zero vector u satisfying u ≤ 2 (d-1)/4 det(L) 1/d .

We will now prove the following general result regarding the modular roots of bivariate homogeneous polynomials which can be of independent interest. Theorem 2. Let f (x, y) ∈ Z[x, y] be a homogeneous polynomial of degree δ with f (x, 0) = x δ , N be a non-zero integer and α be a rational number in [0, 1], then one can retrieve in polynomial time in log N , δ and the bit-size of α, all the rationals x 0 /y 0 , where x 0 and y 0 are integers such that gcd(f (x 0 , y 0 ),

N ) ≥ N α and |x 0 |, |y 0 | ≤ N α 2 /2δ .
Proof. Let b be a divisor of N for which their exists (x 0 , y 0 ) ∈ Z 2 such that b = gcd(f (x 0 , y 0 ), N ) ≥ N α . We define some integral parameters (to be specified later) m, t and t with t = m + t and construct a family of δt + 1 homogeneous polynomials g and h of degree δt such that (x 0 , y 0 ) is a common root modulo b m . More precisely, we consider the following polynomials g i,j (x, y) = x j y δ(t-i)-j f i N m-i for i = 0, . . . , m -1, j = 0, . . . , δ -1 h i (x, y) = x i y δt -i f m for i = 0, . . . , δt .

We build the triangular matrix L of dimension δt + 1, containing the coefficients of the polynomials g i,j and h i . We will apply LLL to the lattice spanned by the rows of L. The columns correspond to the coefficients of the monomials y δt , xy δt-1 , . . . , x δt-1 y, x δt . Let β ∈ [0, 1] such that M = N β . The product of the diagonal elements gives det(L) = N δm(m+1)/2 . If we omit the quantities that do not depend on N , to satisfy the inequality of Lemma 1 with the root bound M , the LLL bound from Lemma 2 implies that we must have

δm(m + 1)/2 ≤ (δt + 1)(αm -δtβ) (3)
and if we set λ such that t = λm, this gives asymptotically β ≤ α δλ -1 2δλ 2 , which is maximal when λ = 1 α , and in this case, β max = α 2 /2δ. The vector output by LLL gives a homogeneous polynomial f (x, y) such that f (x 0 , y 0 ) = 0 thanks to Lemma 1. Let r = x/y, any rational root of the form x 0 /y 0 can be found by extracting the rational roots of f (r) = 1/y δt f (x, y) with classical methods.

For the case we are most interested in, δ = 2, N = pq 2 with p and q of the same size, i. e., α = 2/3 then λ = 3/2 and we can asymptotically get roots up to N β with β = 1 9 . If we take m = 4 and t = 6, i. e., we work with a lattice of dimension 13, we get from (3) that β ≈ 1 10.63 and with a 31-dimensional lattice (m = 10 and t = 15), β ≈ 1 9.62 . If the size of q grows compared to p, i. e., α increases towards 1, then β increases towards 1/4. For example, if q is two times larger than p, i. e., α = 4/5 then β = 1/6.25. For α = 6/7, we get β ≈ 1/5.44.

We will call HomogeneousCoppersmith the algorithm which implements this method. It takes as input an integer N = pq 2 and a binary quadratic form [a, b, c], from which we deduce the unitary polynomial x 2 +b xy +c y 2 , by dividing both b and c by a modulo N , and the parameters m and t. In fact, this method will only disclose proper representations of q 2 , those for which x and y are coprime, but we note that f k properly represents q 2 , and therefore so does our form [a, b, c].

The case α = 1 of Theorem 2 can already be found in Joux's book [START_REF] Joux | Algorithmic Cryptanalysis[END_REF] and we mention that a similar technique has already been independently investigated by Bernstein in [START_REF] Bernstein | List decoding for binary Goppa codes[END_REF].

4 A Õ(p 1/2 )-Deterministic Factoring Algorithm for pq 2
We detail our new quadratic form-based factoring algorithm for numbers of the form pq 2 . In this section, p and q will be of same bit-size, and p ≡ 1 (mod 4).

The Algorithm

Roughly speaking, if ∆ q = N = pq 2 , our factoring algorithm, depicted in Fig. 3, exploits the fact that the non-reduced forms f k = [q 2 , kq, -] reduce to forms fk for which there exists a small pair (x 0 , y 0 ) such that q 2 | fk (x 0 , y 0 ) while q 2 | N . From Theorem 1, we know that these reduced forms appear on the principal cycle of the class group of discriminant ∆ q . To detect them, we start a walk in the principal cycle from the principal form 1 N , and apply Rho until the Coppersmith-like method finds these small solutions.

Input: N = pq 2 , m, t Output: p, q 1. h ← 1N 2. while (x0, y0) not found do 2.1. h ← Rho(h) 2.2. x0/y0 ← HomogeneousCoppersmith(h, N, m, t) 3. q ← Sqrt(Gcd(h(x0, y0), N )) 4. return (N/q 2 , q) Fig. 3. Factoring N = pq 2

Heuristic Correctness and Analysis of Our Algorithm

Assuming Heuristic 2, starting from 1 N , after O(R p ) iterations, the algorithm will stop on a reduced form whose roots will be found with our Coppersmithlike method (for suitable values of m and t) since they will satisfy the expected N 1/9 bound. The computation of gcd(h(x 0 , y 0 ), N ) will therefore expose q 2 and factor N . The time complexity of our algorithm is then heuristically O(R p Poly(log N )), whereas the space complexity is O(log N ). The worst-case complexity is O(p 1/2 log p Poly(log N )). For small regulators, such as in REAL-NICE cryptosystem (see. Subsection 5.1), the time complexity is polynomial.

This algorithm can be generalised with a few modifications to primes p such that p ≡ 3 (mod 4), by considering ∆ q = 4pq 2 . Moreover if the bit-sizes of p and q are unbalanced, our experiments suggest that the size of the roots will be small enough (see end of Subsection 2.3 and Section 3), so the factoring algorithm will also work in this case, with the same complexity.

Comparison with other Deterministic Factorisation Methods. Boneh, Durfee and Howgrave-Graham presented in [START_REF] Boneh | Factoring N = p r q for large r[END_REF] an algorithm for factoring integers N = p r q. Their main result is the following: Lemma 3 ([BDH99]). Let N = p r q be given, and assume q < p c for some c. Furthermore, assume that P is an integer satisfying |P -p| < p 1-c r+c -2 r d . Then the factor p may be computed from N , r, c and P by an algorithm whose running time is dominated by the time it takes to run LLL on a lattice of dimension d.

For r = 2 and c = 1, this leads to a deterministic factoring algorithm which consists in exhaustively search for an approximation P of p and to solve the polynomial equation (P + X) 2 ≡ 0 (mod p 2 ) with a method à la Coppersmith. The approximation will be found after O(p 1/3 ) = O(N 1/9 ) iterations.

The fastest deterministic generic integer factorisation algorithm is actually a version of Strassen's algorithm [START_REF] Strassen | Einige Resultate über Berechnungskomplexität[END_REF] from Bostan, Gaudry and Schost [START_REF] Bostan | Linear Recurrences with Polynomial Coefficients and Application to Integer Factorization and Cartier-Manin Operator[END_REF], who ameliorates a work of Chudnovsky and Chudnovsky [START_REF] Chudnovsky | Approximations and Complex Multiplication According to Ramanujan[END_REF] and proves a complexity of O(M int ( 4 √ N log N )) where M int is a function such that integers of bit-size d can me multiplied in M int (d) bit operations. More precisely, for numbers of our interest, Lemma 13 from [START_REF] Bostan | Linear Recurrences with Polynomial Coefficients and Application to Integer Factorization and Cartier-Manin Operator[END_REF] gives the precise complexity: In particular, for b = N 1/3 , the complexity is Õ(N 1/6 ), with a very large space complexity compared to our algorithm. Moreover, none of these two last of algorithms can actually factor an integer of cryptographic size. The fact that a prime divisor has a small regulator does not help in these algorithms, whereas it makes the factorisation polynomial in our method.

Cryptanalysis of the NICE Cryptosystems

Hartmann, Paulus and Takagi proposed the elegant NICE encryption scheme (see [HPT99,PT99,PT00]), based on imaginary quadratic fields and whose main feature was a quadratic decryption time. Later on, several other schemes, including (special) signature schemes relying on this framework have been proposed. The public key of these NICE cryptosystems contains a discriminant ∆ q = -pq 2 together with a reduced ideal h whose class belongs to the kernel of φq . The idea underlying the NICE cryptosystem is to hide the message behind a random element [h] r of the kernel. Applying φq will make this random element disappear, and the message will then be recovered.

In [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF], Jacobson, Scheidler and Weimer embedded the original NICE cryptosystem in real quadratic fields. Whereas the idea remains essentially the same as the original, the implementation is very different. The discriminant is now ∆ q = pq 2 , but because of the differences between imaginary and real setting, these discriminant will have to be chosen carefully. Among these differences, the class numbers are expected to be small with very high probability (see the Cohen-Lenstra heuristics [START_REF] Cohen | Heuristics on class groups[END_REF]). Moreover, an equivalence class does not contain a unique reduced element anymore, but a multitude of them, whose number is governed by the size of the fundamental unit. The rough ideas to understand these systems and our new attacks are given in the following. The full description of the systems is omitted for lack of space but can be found in [HPT99, JSW08].

Polynomial-Time Key Recovery in the Real Setting

The core of the design of the REAL-NICE encryption scheme is the very particular choice of the secret prime numbers p and q such that ∆ K = p and ∆ q = pq 2 . They are chosen such that the ratio R ∆q /R ∆ K is of order of magnitude of q and that R ∆ K is bounded by a polynomial in log(∆ K ). To ensure the first property, it is sufficient to choose q such that q -∆ K q is a small multiple of a large prime. If the second property is very unlikely to naturally happen since the regulator of p is generally of the order of magnitude of √ p, it is indeed quite easy to construct fundamental primes with small regulator. The authors of [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF] suggest to produce a prime p as a so-called Schinzel sleeper, which is a positive squarefree integer of the form p = a 2 x 2 + 2bx + c with a, b, c, x in Z, a = 0 and b 2 -4ac dividing 4 gcd(a 2 , b) 2 . Schinzel sleepers are known to have a regulator of the order log(p) (see [START_REF] Cheng | Some Results Concerning Certain Periodic Continued Fractions[END_REF]). Some care must be taken when setting the (secret) a, b, c, x values, otherwise the resulting ∆ q = pq 2 is subject to factorisation attacks described in [START_REF] Weimer | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]. We do not provide here more details on these choices since the crucial property for our attack is the fact that the regulator is actually of the order log(p). The public key consists of the sole discriminant ∆ q . The message is carefully embedded (and padded) into a primitive O ∆q -ideal so that it will be recognised during decryption. Instead of moving the message ideal m to a different equivalence class (like in the imaginary case), the encryption actually hides the message in the cycle of reduced ideal of its own equivalent class by multiplication of a random principal O ∆q -ideal h (computed during encryption). The decryption process consists then in applying the (secret) map φq and perform an exhaustive search for the padded message in the small cycle of φq ([mh]). This exhaustive search is actually possible thanks to the choice of p which has a very small regulator. Like in the imaginary case, the decryption procedure has a quadratic complexity and significantly outperforms an RSA decryption for any given security level (see Table 3 from [START_REF] Jacobson | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]).

Unfortunately, due to the particular but necessary choice of the secret prime p, the following result states the total insecurity of the REAL-NICE system.

Result 1 Algorithm 3 recovers the secret key of REAL-NICE in polynomial time in the security parameter under Heuristic 2 since the secret fundamental discriminant p is chosen to have a regulator bounded by a polynomial in log p.

We apply the cryptanalysis on the following example. The Schinzel polynomial S(X) = 2725 2 X 2 + 2 • 3815X + 2 produces a suitable 256-bit prime p for the value X 0 = 103042745825387139695432123167592199. This prime has a regulator R ∆ K 90.83. The second 256-bit prime q is chosen following the recommendations from [START_REF] Weimer | An Adaptation of the NICE Cryptosystem to Real Quadratic Orders[END_REF]. This leads to a the discriminant ∆q = 28736938823310044873380716142282073396186843906757463274792638734144060602830510 80738669163489273592599054529442271053869832485363682341892124500678400322719842 63278692833860326257638544601057379571931906787755152745236263303465093

Our algorithm recovers the prime q = 60372105471499634417192859173853663456123015267207769653235558092781188395563 from ∆ q after 45 iterations in 42.42 seconds on a standard laptop. The rational root is x0 y0 equal to -2155511611710996445623 3544874277134778658948 , where x 0 and y 0 satisfy

Polynomial-Time Key Recovery of the Original NICE

As mentioned above, the public key of the original NICE cryptosystem contains the representation of a reduced ideal h whose class belongs to the kernel of the surjection φq . The total-break of the NICE cryptosystem is equivalent to solving the following kernel problem.

Definition 6 (Kernel Problem [START_REF] Biehl | Efficient Undeniable Signature Schemes based on Ideal Arithmetic in Quadratic Orders[END_REF]). Let λ be an integer, p and q be two λ-bit primes with p ≡ 3 (mod 4). Fix a non-fundamental discriminant ∆ q = -pq 2 . Given an element [h] of ker φq , factor the discriminant ∆ q .

Castagnos and Laguillaumie proposed in [START_REF] Castagnos | On the Security of Cryptosystems with Quadratic Decryption: The Nicest Cryptanalysis[END_REF] a polynomial-time algorithm to solve this problem. We propose here a completely different solution within the spirit of our factorisation method and whose complexity is also polynomial-time. As discuss in Subsection 2.3, the idea is to benefit from the fact that the public ideal h corresponds to a reduced quadratic form, fk , which represents q 2 . We thus find these x 0 and y 0 such that gcd( fk (x 0 , y 0 ), ∆ q ) = q 2 with the Coppersmith method of Section 3.

Fig. 1 .

 1 Fig. 1. Probability that |M k | < |∆q| 1/9 in function of the bit-size λ of p and q

  Lemma 4 ([BGS07]). Let b, N be two integers with 2 ≤ b < N . One can compute a prime divisor of N bounded by b, or prove that no such divisor exists in O M int ( √ b log N ) + log bM int (log N ) log log N bit operations and space O( √ b log N ) bits.

Acknowledgements. We warmly thank Denis Simon and Brigitte Vallée for helpful discussions and the reviewers for their useful comments. Part of this work was supported by the Commission of the European Communities through the ICT program under contract ICT-2007-216676 ECRYPT II.

This work was done while

Result 2 The Homogeneous Coppersmith method from Section 3 solves the Kernel Problem in polynomial time in the security parameter under Heuristic 1.

We apply our key recovery on the example of NICE proposed in [JJ00,CL09]: The public key consists in ∆ q and h = (a, b). Our Coppersmith method finds in less that half a second the root u 0 = -103023911 349555951 = x0 y0 and

All our experiments have been run on a standard laptop under Linux with software Sage. The lattice reduction have been performed with Stehlé's fplll [Ste].


Carole Le Guyader

Dominique Apprato

Christian Gout

Spline approximation of gradient fields: applications to wind velocity fields

Keywords: Approximation, smoothing problems, vectorial splines, convergence, Sobolev space, conservative vector fields AMS subject classification: 41A15, 41A25, 68U10

We study a spline-based approximation of vector fields in the conservative case. This problem appears for instance when approximating current or wind velocity fields, the data deriving in those cases from a potential (pressure for the wind, etc..). In the modelling, we introduce a minimization problem on an Hilbert space for which the existence and uniqueness of the solution are provided. A convergence result in the introduced Sobolev space is established using norm equivalence and compactness arguments, as well as an approximation error estimate of the involved smoothing D m -splines.

Introduction

The problem of vector field approximation emerges in a wide range of fields such as: motion control, computer vision, geometrical analysis, geometrical design, analysis of acoustic and/or electromagnetic waves, as well as in geophysics, medical imaging, fluid mechanics and so on... Many different approaches have been introduced to solve each specific problem occuring in the above fields of investigation. One can mention for instance: finite element approximation (see Dzhabrailov et al. [START_REF] Dzhabrailov | The finite element approximation of vector fields in curvilinear coordinates[END_REF]), PDE-based methods (see Amodei and Benbourhim [START_REF] Amodei | A Vector Spline Approximation[END_REF]), spline and Rational Basis Function approximations (see Awanou and Lai [START_REF] Awanou | Trivariate spline approximations of 3D Navier-Stokes equations[END_REF], Dodu and Rabut [START_REF] Dodu | Vectorial interpolation using radial-basis-like functions[END_REF], Benbourhim and Bouhamidi [START_REF] Benbourhim | Approximation of vector fields by thin plate splines with tension[END_REF]- [START_REF] Benbourhim | Pseudo-polyharmonic vectorial approximation for divcurl and elastic semi-norms[END_REF], Ettl-Lowitsch et al. [START_REF] Ettl-Lowitzsch | Shape reconstruction from gradient data[END_REF][START_REF] Lowitzsch | Approximation and Interpolation employing divergence-free radial basis functions with applications[END_REF][START_REF] Lowitzsch | Error estimate for matrix-valued radial basis functions interpolation[END_REF], Zhou et al. [START_REF] Zhou | Energy Minimization Method for Scattered Data Hermite Interpolation[END_REF]...).

The originality of this work consists:

• in considering that the vector field derives from a potential (conservative vector field): it occurs for instance in meteorology (winds derive from temperature potentials), oceanography (currents derive from pression potentials), image processing...

• of a rigorous study of existence-uniqueness of the solution of the problem phrased as an energy minimization,

• in establishing a convergence result (while many approaches only give algorithms without mathematical study) and providing an approximation error estimate.

From the best of our knowledge, our approach is the first one of this kind. Moreover, we propose a method easy to implement and parallelize, using classical tools like the least-square method and the finite element method (for the discretization). More precisely, in this work, we do not want to compute a potential that could generate the vector field data. We only want to get a global approximation of the vector field dataset on a bounded domain, taking into account in the modelling that this approximation derives from a potential. Furthermore, contrary to interpolation methods, we prefer to fit the vector field dataset in the case of realistic data (when the number of vectors is large or when the data are corrupted by noise). To achieve this, we introduce a minimization problem defined as a regularized leastsquare problem formulated on a Sobolev space of potentials (see subsection 2.1). Obviously, this problem has an infinite number of solutions, but we derive from it a problem expressed in terms of the gradient vectors. We prove that the associated problem in terms of vectors has a unique solution which is the corresponding approximation of the vector field dataset (subsection 2.2). In subsection 2.3, we give a convergence result when the number of vectors increases to infinity. In Section 3, we give the discretization complemented by an approximation error estimate of the involved smoothing D m -splines. We then give numerical results in Section 4.

Let Ω be a non-empty connected bounded open subset of R n with Lipschitz boundary. Let A = a i i=1,••• ,N be a set of N points of Ω which contains a P m-1 -unisolvent subset. Let also w i i=1,••• ,N be a set of N vectors of R n that corresponds to the vector field dataset. We aim at finding a function Φ : Ω → R smooth enough such that for i

= 1, • • • , N , ∇Φ(a i ) ≃ w i ∈ R n .
Let us denote by ρ the operator defined by:

ρ : H m (Ω, R n ) → (R n ) N v → ρ(v) = (v(a 1 ), • • • , v(a N )) T . (1)
Let us assume that m > n 2 so that the following Sobolev 's embeddings hold:

H m (Ω, R n ) C 0 ( Ω, R n ) (2) H m+1 (Ω, R) C 1 ( Ω, R) (3)
One

sets ∀ξ ∈ (R n ) N and ∀η ∈ (R n ) N , ξ, η N = N i=1 ξ i , η i n , where •, • n denotes the Euclidean scalar product in R n and ∀ξ ∈ (R n ) N , ξ N = ξ, ξ 1 2 N . 2
2. Modelling

The fitting model

We first introduce a regularized least-square problem defined on a space of potentials (realvalued functions) to fit the vector field dataset. For any ε > 0, we introduce the functional J ε defined as follows:

J ε : H m+1 (Ω, R) → R v → ρ(∇v) -w 2 N + ε|v| 2 m+1,Ω,R , (4)
where

w = (w 1 , • • • , w N ) T ∈ (R n ) N is the vector field dataset and | • | m+1,Ω,R , the semi-norm on H m+1 (Ω, R).
This kind of approach is related to the smoothing D m -spline for surface approximation introduced by Arcangéli [START_REF] Arcangéli | Multidimensional minimizing splines. Theory and applications[END_REF] (see also Duchon [12] for a general introduction or Gout [START_REF] Gout | C k surface reconstruction from surface patches[END_REF], and López de Silanes and Arcangéli [START_REF] López De Silanes | Estimations de l'erreur d'approximation par splines d'interpolation et d'ajustement d'ordre (m, s), (French) [Approximation error estimates for interpolating and smoothing (m, s)-splines[END_REF][START_REF] López De Silanes | Sur la convergence des D m -splines d'ajustement pour des données exactes ou bruitées[END_REF] for convergence results). We consider the following fitting problem using potentials:

Search for σ ε ∈ H m+1 (Ω, R) such that ∀v ∈ H m+1 (Ω, R), J ε (σ ε ) ≤ J ε (v), (5)
where J ε is defined in [START_REF] Arcangéli | Multidimensional minimizing splines. Theory and applications[END_REF]. Obviously, the minimization problem [START_REF] Arcangeli | Sur l'évaluation de l'erreur d'interpolation de Lagrange dans un ouvert de R n[END_REF] does not have a unique solution, the potential function σ ε being defined up to a constant.

Remark 2.1. An alternative model could be introduced to circumvent the lack of uniqueness of the potential function by including a few interpolation constraints in the modelling. Indeed, it would guarantee uniqueness of the solution expressed in terms of the potential function. In this case, the problem is phrased as a functional minimization problem on a convex subset K (or a vector space if we retain only some positions for which the potential function is zero) of the Hilbert space H m+1 (Ω, R) as follows:

Search for σ ǫ ∈ K such that ∀v ∈ K, J ǫ (σ ǫ ) ≤ J ǫ (v), with K = v ∈ H m+1 (Ω, R), ∀i ∈ {1, • • • , M } , v(b i ) = β i and β i given at b i for all i ∈ {1, • • • , M }.
Existence and uniqueness of the solution are obtained by applying Stampacchia's theorem (after rewriting the functional to be minimized and using norm equivalence arguments resulting from Necas theorem application). The introduction of Lagrange's multipliers allows to obtain the variational formulation of the problem on the whole space H m+1 (Ω, R) instead of a variational inequality.

The main hindrance of this modelling comes from the fact that σ ǫ ∈ H m+1 (Ω, R). In practice m = 2 in order to have

H 3 (Ω, R) C 1 ( Ω, R), yielding σ ǫ ∈ H 3 (Ω, R).
For instance, if we consider as generic finite element the C 2 Bogner-Fox-Schmit element, there are 9 basis functions per node, leading to block matrices of size 81 × 81 to compute the stiffness matrix, which can turn out to be tricky from a numerical viewpoint. The idea of reformulating the problem in terms of the gradient vector field itself avoids this issue.

So, in order to overcome the lack of uniqueness, we introduce the following problem expressed in terms of the gradients of the potentials: as ε is a positive arbitrary parameter and as |.| m+1,Ω,R and |∇.| m,Ω,R n are equivalent in H m+1 (Ω, R), we derive from problem [START_REF] Arcangeli | Sur l'évaluation de l'erreur d'interpolation de Lagrange dans un ouvert de R n[END_REF], the following problem of minimization:

Search for σ ε ∈ H m+1 (Ω, R) such that ∀v ∈ H m+1 (Ω, R), ρ(∇σ ε ) -w 2 N + ε|∇σ ε | 2 m,Ω,R n ≤ ρ(∇v) -w 2 N + ε|∇v| 2 m,Ω,R n . (6)
This formulation leads to the following problem associated to problem [START_REF] Awanou | Trivariate spline approximations of 3D Navier-Stokes equations[END_REF]. To state this new problem, we first define the following functional denoted by F ε :

F ε : H m (Ω, R n ) → R f → ρ(f ) -w 2 N + ε|f | 2 m,Ω,R n , (7)
and we consider the following problem stated by:

Search for u ε ∈ H m (Ω, R n ) such that ∀v ∈ H m (Ω, R n ), F ε (u ε ) ≤ F ε (v). ( 8
)
This formulation is motivated by the following theorem (in the proof, we use the connectedness of Ω and the unisolvence property):

Theorem 2.2. If σ 1 ε and σ 2 ε are two distinct solutions of problem (6) then:

∇σ 1 ε = ∇σ 2 ε = u ε , (9)
where u ε is the unique solution of problem (8).

Existence and uniqueness results

We first give the following variational problem (which is equivalent to (8)):

Search for u ε ∈ H m (Ω, R n ) such that ∀v ∈ H m (Ω, R n ), ρ(u ε ), ρ(v) N + ε(u ε , v) m,Ω,R n = w, ρ(v) N . (10)
We now give a lemma that will be of interest in the following.

Lemma 2.3. The mapping defined by:

|| • || A,m,Ω,R n :        H m (Ω, R n ) → R f → ||f || A,m,Ω,R n = ρ(f ) 2 N + |f | 2 m,Ω,R n 1 2 = N i=1 f (a i ) 2 n + |f | 2 m,Ω,R n 1 2 (11) is a Hilbert norm equivalent to the norm ||f || m,Ω,R n in H m (Ω, R n ) (H m (Ω, R n ) is
∀f ∈ H m (Ω, R n ), ∀a i ∈ A, i = 1, • • • , N , f (a i ) n ≤ ||f || C 0 ( Ω,R n ) ≤ c||f || H m (Ω,R n ) (Sobolev ′ s embedding), so ||f || A,m,Ω,R n ≤ (1 + c 2 N ) 1 2 ||f || m,Ω,R n .
Let us take k = m and p = 2 in Necas' theorem (Chapter 2, section 7.1 from [START_REF] Necas | Les méthodes directes en théorie des équations elliptiques[END_REF]) and let us take ρ as functional f i . Then, using the property of unisolvence of the set A, we obtain from Necas' theorem that there exists a positive constant c 1 such that:

||f || m,Ω,R n ≤ 1 c 1 |f | 2 m,Ω,R n + N i=1 f (a i ) 2 n 1 2 = 1 c 1 ||f || A,m,Ω,R n ,
which concludes the proof.

Using Lax-Milgram lemma and norm equivalence results, we can now establish the following theorem:

Theorem 2.4. Variational problem (10) admits a unique solution.

Convergence result

In the following, we consider classical convergence hypotheses (see [START_REF] Arcangéli | Multidimensional minimizing splines. Theory and applications[END_REF] for more details).

Let D be a subset of R + * for which 0 is an accumulation point.

For any d ∈ D, let A d be a set of N = N (d) distinct points from Ω that contains a P m-1 -unisolvent subset. We assume that

sup x∈Ω δ(x, A d ) = d, (12)
where δ is the Euclidean distance in R n . Thus d is the radius of the biggest sphere included in Ω that contains no point from A d (Hausdorff distance). Also d is bounded and

lim d→0 sup x∈Ω δ(x, A d ) = 0. ( 13
)
For any d ∈ D, let us denote by ρ d the mapping defined by:

ρ d :    H m (Ω, R n ) → (R n ) N v → ρ d (v) = (v(a)) a∈A d T .
Introducing the norm || • || A d ,m,Ω,R n defined by:

||f || A d ,m,Ω,R n = ρ d (f ) 2 N + |f | 2 m,Ω,R n 1 2 , as shown in Lemma 2.3, the norm || • || A d ,m,Ω,R n is equivalent to the norm || • || m,Ω,R n in H m (Ω, R n ).
To prove the convergence result, we need the following lemma.

Lemma 2.5.

Let A 0 = {b 01 , b 02 , • • • , b 0ℵ } be a fixed P m-1 -unisolvent subset of Ω (in this case, ℵ = dim P m-1 ).
By hypothesis, 0 ∈ D and (13) holds so:

∀j = 1, • • • , ℵ, ∃ a d 0j d∈D , ∀d ∈ D, a d 0j ∈ A d and b 0j = lim d→0 a d 0j . ( 14
)
For any d ∈ D, let A d 0 be the set a d 01 , • • • , a d 0ℵ and let || • || A d 0 ,m,Ω,R n be the mapping defined by: ∀f ∈ H m (Ω, R n ), ||f || A d 0 ,m,Ω,R n =   ℵ j=1 f (a d 0j ) 2 n + |f | 2 m,Ω,R n   1 2
.

Then, there exists η > 0 such that for any d ≤ η, the set

A d 0 is P m-1 -unisolvent and || • || A d 0 ,m,Ω,R n is a norm on H m (Ω, R n ) uniformly equivalent over D∩]0, η] to the norm || • || m,Ω,R n .

Proof

The Sobolev's embedding gives that:

∃C 1 > 0, ∀d ∈ D, ∀f ∈ H m (Ω, R n ), ||f || A d 0 ,m,Ω,R n ≤ C 1 ||f || m,Ω,R n . The constant C 1 is independent of d and depends on ℵ which is the dimension of P m-1 . Let us now find a constant C 2 independent of d such that the inequality ||f || m,Ω,R n ≤ C 2 ||f || A d 0 ,m,Ω,R n holds. For f ∈ H m (Ω, R n ), 1 2 ℵ j=1 f (b 0j ) 2 n = 1 2 ℵ j=1 f (b 0j ) -f a d 0j + f a d 0j 2 n , ≤ ℵ j=1 f (b 0j ) -f a d 0j 2 n + ℵ j=1 f a d 0j 2 n . (15)
The open subset Ω having a Lipschitz-continuous boundary and m being such that m > n 2 , the space H m (Ω, R n ) satisfies the Sobolev's Hölder embedding theorem (continuous injection), namely:

∃λ ∈]0, 1], H m (Ω, R n ) C 0,λ ( Ω, R n ). Thus f ∈ C 0,λ ( Ω, R n ) and ∃C > 0, ∀j = 1, • • • , ℵ, ∀d ∈ D, f (b 0j ) -f a d 0j 2 n ≤ ||f || 2 C 0,λ ( Ω,R n ) b 0j -a d 0j 2λ n ≤ C 2 ||f || 2 m,Ω,R n b 0j -a d 0j 2λ n . (16)
Besides, from ( 14) it comes:

∀j = 1, • • • , ℵ, ∀β j > 0, ∃η βj > 0, ∀d ∈ D, d ≤ η βj ⇒ a d 0j -b 0j n ≤ β j . Then ∀j = 1, • • • , ℵ, ∀β j > 0, ∃η βj > 0, ∀d ∈ D, d ≤ η βj ⇒ f (b 0j ) -f a d 0j 2 n ≤ C 2 ||f || 2 m,Ω,R n β 2λ j .
Let β > 0 and let us take

β j = β, ∀j = 1, • • • , ℵ and η = min (η β1 , • • • , η β ℵ ) then ∀d ∈ D,   d ≤ η ⇒ ℵ j=1 f (b 0j ) -f a d 0j 2 n ≤ C 2 β 2λ ℵ||f || 2 m,Ω,R n   .
This implies that:

∀β > 0, ∃η > 0, ∀d ∈ D, ∀f ∈ H m (Ω, R n ),   d ≤ η ⇒ 1 2 ℵ j=1 f (b 0j ) 2 n + |f | 2 m,Ω,R n -C 2 β 2λ ℵ||f || 2 m,Ω,R n ≤ ||f || 2 A d 0 ,m,Ω,R n   .
As previously shown, the mapping

f ∈ H m (Ω, R n ) →   1 2 ℵ j=1 f (b 0j ) 2 n + |f | 2 m,Ω,R n   1 2
is a norm on

H m (Ω, R n ) equivalent to the norm || • || m,Ω,R n so: ∀β > 0, ∃η > 0, ∀d ∈ D, ∀f ∈ H m (Ω, R n ), d ≤ η ⇒ C ′ 2 -C 2 β 2λ ℵ ||f || 2 m,Ω,R n ≤ ||f || 2 A d 0 ,m,Ω,R n .
By choosing β adequately, the norm equivalence is obtained.

Theorem 2.6. Assume that there exists a function f ∈ H m (Ω, R n ) such that for any d ∈ D: ρ d ( f ) = w = w d (we recall that w = w d is the given data set -a gradient vector set, that is, the input of our problem), and ε = ε(d) ∈]0, ε 0 ], ε 0 > 0. For any d ∈ D, we denote by u d ε the unique solution of problem [START_REF] Benbourhim | Pseudo-polyharmonic vectorial approximation for divcurl and elastic semi-norms[END_REF]. Then under the above assumptions we have:

lim d→0 ||u d ε -f || m,Ω,R n = 0. ( 17
)
Proof 2.5.

The proof is divided into four steps that we detail hereafter.

First step: We start by proving that the sequence

u d ε d∈D∩]0,η] ε∈]0,ε0]
is bounded in H m (Ω, R n ) so that one can extract a subsequence u d l ε l l∈N with lim l→+∞ d l = 0 (since 0 is an accumulation point of D)

and ε l ∈]0, ε 0 ], ∀l ∈ N (-we assume that ε = ε(d), so ε l = ε(d l )-) that weakly converges to an element of H m (Ω, R n ) denoted by f * .
Second step: In the second step of the proof, we prove that f * = f .

Third step: Then we aim at proving that the sequence u d l ε l l∈N strongly converges to f in H m (Ω, R n ). The Rellich-Kondrašov compact embedding theorem gives that:

∀r, r ′ ∈ R, r > r ′ , H r (Ω, R n )⊂H r ′ (Ω, R n ) (compact injection).
In our case, it means that the sequence u d l ε l l∈N that weakly converges to

f in H m (Ω, R n ), strongly converges to f in H m-1 (Ω, R n ). Therefore, as ||u d l ε l -f || 2 m,Ω,R n = ||u d l ε l -f || 2 m-1,Ω,R n + |u d l ε l -f | 2 m,Ω,R n , we just need to prove that lim l→∞ |u d l ε l -f | m,Ω,R n = 0.
One has:

|u d l ε l -f | 2 m,Ω,R n = |u d l ε l | 2 m,Ω,R n + | f | 2 m,Ω,R n -2(u d l ε l , f ) m,Ω,R n , |u d l ε l -f | 2 m,Ω,R n ≤ 2| f | 2 m,Ω,R n -2(u d l ε l , f ) m,Ω,R n , ∀α ∈ N n , |α| = m, ∂ α u d l ε l weakly converges to ∂ α f in L 2 (Ω, R n ). So taking as test function f , it follows that lim l→+∞ (u d l ε l , f ) m,Ω,R n = | f | 2 m,Ω,R n . Consequently, lim l→∞ |u d l ε l -f | m,Ω,R n = 0 (18)
and finally

lim l→∞ ||u d l ε l -f || m,Ω,R n = 0. ( 19
)
Fourth step: Let us assume that ||u d εf || m,Ω,R n does not tend to 0 when d tends to 0. It means that there exists a real number α > 0 and two sequences (d k ) k∈N and (ε k ) k∈N such that

d k →0 when k → +∞ and ε k = ε(d k ) and ∀k ∈ N, ||u d k ε k -f || m,Ω,R n > α. ( 20
)
Following the same steps as previously done, there exists a subsequence of (u d k ε k ) k∈N that strongly converges to f in H m (Ω, R n ). This last point is in contradiction with (20).

Discretization

Introduction of the notations

In this part, we give a discretization of the variational problem [START_REF] Clément | Approximation by Finite Element Functions Using Local Regularisation[END_REF]. We choose here a finite element space and we need the following notations.

• ∆ is the dimension of P m-1 , space of n-variable polynomials with total degree less than m -1.

• Let H be a bounded subset of R * + for which 0 is an accumulation point.

• Let Ω be a bounded polyhedral open subset of R n containing Ω.

• Let m ′ be an integer such that m ′ > m. We recall that m > n 2 .

• For any h ∈ H, let Th be a triangulation of Ω by n-simplices or n-rectangles K with diameter h K ≤ h and let Ṽh n be a finite-element-type space constructed on Th , finite-dimensional subspace of H m ( Ω, R n ). We assume that the family Th h∈H is regular.

• ∀h ∈ H, let Ω h be the open set defined as the interior of the union of elements K ∈ Th such that K ∩ Ω = ∅. It can be seen that ∀h ∈ H, Ω ⊂ Ω h ⊂ Ω and lim h→0 mes (Ω h \Ω) = 0.

• For any h ∈ H, let (V h ) n be the space of restrictions to Ω h of functions in Ṽh

n . It is a finite dimensional subspace of H m (Ω h , R n ) with dimension M = M (h).
We consider the variational formulation:

   Find u d ε h ∈ (V h ) n such that ∀v h ∈ (V h ) n , ρ d u d ε h , ρ d (v h ) N =N (d) + ε u d ε h , v h m,Ω h ,R n = w d , ρ d (v h ) N , (21)
which is a discretisation of problem [START_REF] Clément | Approximation by Finite Element Functions Using Local Regularisation[END_REF] and which admits a unique solution u d ε h . We aim to establish an error estimate for the solution u d ε h of problem [START_REF] Clément | Approximation by Finite Element Functions Using Local Regularisation[END_REF]. We mainly use results from Duchon ( [START_REF] Duchon | Sur l'erreur d'interpolation des fonctions de plusieurs variables par les D m splines[END_REF]) and Apprato and Lopez de Silanes ( [START_REF] Lopez De Silanes | Estimations de l'erreur d'approximation sur un domaine borné de R n par D m splines d'interpolation et d'ajustement discrètes[END_REF]). We first recall two results from Duchon ([13]). For any t ∈ R n , and for R > 0, we denote by B (t, R) the closed ball centered at t and with radius R. We have the following propositions (see Duchon [START_REF] Duchon | Sur l'erreur d'interpolation des fonctions de plusieurs variables par les D m splines[END_REF] for the proofs): Proposition 3.1. Let Ω be an open set of R n with Lipschitz boundary. Then there exist constants M > 1, M 1 > 1 and λ 0 > 0 such that, for any λ ∈ ]0, λ 0 ], there exists T λ ⊂ Ω such that:

         i) ∀t ∈ T λ , B (t, λ) ⊂ Ω ii) Ω ⊂ ∪ t∈T λ B (t, M λ) iii) t∈T λ 1 B(t,M λ) ≤ M 1
where 1 E is the characteristic function of the set E.

Proposition 3.2. There exist R > 0 (depending on n and m) and (for any M ≥ 1 and any integer l such that l ≤ m) C > 0 (depending on M , m, n and l) satisfying: for any d > 0 and any t ∈ R n , the ball B (t, Rd) contains ∆ balls β 1 , β 2 , ..., β ∆ of radius d such that:

|u| l,B(t,M Rd) ≤ Cd m-l |u| m,B(t,M Rd)
for all u ∈ H m (B (t, M Rd)) that vanishes at at least one point of each ball

β i , i ∈ {1, • • • , ∆}.
The proof shows that the ball B(t, Rd) results from the mapping of a ball B( a, R) by an affine transformation.

Let

γ 0 1 , γ 0 2 , • • • , γ 0 ∆ be a P m-1 -unisolvent set. The set of P m-1 -unisolvent ∆-tuples being an open set in (R n ) ∆ , there exists (see Duchon ([13])) δ > 0 such that {γ 1 , γ 2 , ..., γ ∆ } is still P m-1 - unisolvent if: ∀i = 1, • • • , ∆, γ i -γ 0 i ≤ δ.
By applying an homothety with ratio 1 δ and setting

α i = 1 δ γ 0 i , we obtain ∆ balls B ( α i , 1): the product ∆ i=1 B ( α i , 1) is a compact set of (R n ) ∆ made up of P m-1 -unisolvent ∆-tuples. The set ∆ i=1 B ( α i , 1
) is bounded and thus included in a ball B ( a, R).

Estimation of the approximation error fu d ε h|Ω

Let f ∈ H m ′ (Ω, R n ), m ′ > m and let u d ε h|Ω (u d ε h in the sequel) be the solution of problem (21) with w d = ρ d f . For any t ∈ T Rd ∈ Ω, we define the mapping F t : x → t + d(xa), which is affine bijective and which maps the ball B( a, M R) into the ball B(t, M Rd).

Let { a 1 , a 2 , • • • , a ∆ } be a subset of B ( a, R) such that a i ∈ B ( α i , 1) and F t ( a i ) = a d i ∈ A d , for i = 1, • • • , ∆.
This is made possible since the mapping F t maps the ball B ( α i , 1) into a ball of radius d and since d is defined by [START_REF] Duchon | Splines minimizing rotation-invariant semi-norms in Sobolev spaces[END_REF]. Note that the set a d 1 , • • • , a d ∆ is P m-1 -unisolvent. This property results from the closing remarks of subsection 3.1.

Let us set v

ε = f -u d ε h and v ε = P v ε , where P is the m-extension operator of H m (Ω, R n ) into H m (R n , R n ). For any t ∈ T Rd , let Π t be the P m-1 -Lagrange interpolation operator such that for v ε ∈ H m (B(t, M Rd), R n ), Π t v ε ∈ P m-1 and Π t v ε (a d i ) = v ε (a d i ), i = 1, • • • , ∆.
We have the following propositions (see Lopez de Silanes and Apprato ( [START_REF] Lopez De Silanes | Estimations de l'erreur d'approximation sur un domaine borné de R n par D m splines d'interpolation et d'ajustement discrètes[END_REF]) for the adaptation of the proofs to our problem):

Proposition 3.3. There exist positive constants R, λ 0 , θ and C such that for all d ≤ min λ 0 R , θ , for all ε and for all h, one has:

∀l = 0, • • • , m, t∈T Rd |Π t v ε | 2 l,B(t,M Rd),R n ≤ C 2 d n-2l ε.
Proposition 3.4. There exist positive constants R, λ 0 , η and C such that for all d ≤ min λ 0 R , η , for all ε and for all h, one has:

∀l = 0, • • • , m, t∈T Rd | v ε -Π t v ε | 2 l,B(t,M Rd),R n ≤ C 2 |v ε | 2 m,Ω,R n d 2(m-l) . Theorem 3.5. Let f ∈ H m ′ (Ω, R n ), m ′ >
m and let u d ε h be the solution of ( 21) with w d = ρ d ( f ). Then there exist positive constants R, λ 0 , θ, η and C such that for all d ≤ min λ 0 R , θ, η , for all ε and for all h, one has:

∀l = 0, • • • , m, | f -u d ε h | l,Ω,R n ≤ C | f -u d ε h | m,Ω,R n d m-l + d n 2 -l ε 1 2
.

As a consequence, when u d ε h converges to f in H m (Ω, R n ) (which is the case under some additional assumptions), one has:

∀l = 0, • • • , m -1, | f -u d ε h | l,Ω,R n = o(d m-l ) + O(d n 2 -l ε 1 2 ), | f -u d ε h | l,Ω,R n = o(1) + O(d n 2 -l ε 1 2 ) since | f -u d ε h | l,Ω,R n ≤ C || f -u d ε h || m,Ω,R n + d n 2 -l ε 1 2 .
This error estimate is only significant when the term d n/2-l ε 1/2 tends to zero. Let us note that Theorem 3.5 only gives a convergence result based on the density (d) of the given data samples. Of course, it is also possible to link the finite element parameter h with a convergence result, and therefore to link both d and h in a theoretical convergence result.

In effect, from Theorem 3.5, it is possible to obtain error estimates (depending on h) on fu d ǫh l,Ω,R n : l = 0, ..., m. To do that, we need the following hypotheses:

∃C > 0, ∃Θ > 0, ∀d ∈ D, d ≤ Θ ⇒ N (d) ≤ C d n . ( 22
)
This hypothesis implies an asymptotic distribution of the data on Ω (see Arcangeli et al. [START_REF] Arcangéli | Multidimensional minimizing splines. Theory and applications[END_REF]).

The generic finite elements (K, P K , Σ K ) of V h h∈H must satisfy

P K ⊃ P m ′ -1 (Ω) (23)
and

∀v ∈ H m ′ ( Ω, R n ), ∃ (v h ) h∈H ⊂ V h n , ∀h ∈ H, ∀l = 0, ..., m, |v -v h | l,Ω,R n ≤ Ch m ′ -l |v| m ′ , e Ω,R n . ( 24
)
Let us note that inequation ( 24) is needed to establish the convergence of the approximation but it does not require the usual regularity hypothesis of the finite element method:

H m ′ ( Ω, R n ) ֒→ C s ( Ω, R n ).
Of course, hypothesis [START_REF] Strang | Approximation in the Finite Element Method[END_REF] requires that the triangulation is regular (see Ciarlet [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]). These hypotheses [START_REF] Necas | Les méthodes directes en théorie des équations elliptiques[END_REF][START_REF] Strang | Approximation in the Finite Element Method[END_REF] are satisfied in the usual finite element framework (see Clément [START_REF] Clément | Approximation by Finite Element Functions Using Local Regularisation[END_REF], Strang [START_REF] Strang | Approximation in the Finite Element Method[END_REF]).

We also consider f ∈ H m ′ ( Ω, R n ) a m ′extension of f on Ω (existing since Ω has a Lipschitz boundary):

f |Ω = f in H m ′ (Ω, R n ), ( 25
)
with

H m ′ (Ω, R n ) C 0 (Ω, R n ).
Then we have the following theorem:

Theorem 3.6. We keep the hypotheses of Theorem 3.5, and we suppose satisfied the hypotheses [START_REF] Monier | Dynamic of the currents in the Mozambique canal[END_REF][START_REF] Necas | Les méthodes directes en théorie des équations elliptiques[END_REF][START_REF] Strang | Approximation in the Finite Element Method[END_REF]. Then, there exists C > 0 such that: for any d ≤ min λ0 R , Θ, η , for any ε > 0, and for any h ∈ H, we have:

f -u d ǫh l,Ω,R n ≤ C h m ′ -m d m-l-n/2 √ ε + h m ′ -m + √ 2 f m ′ , e Ω,R n d m-l + d -l+n/2 ε 1/2 ,
where f is introduced in (25).

Remark 3.7. Let us note that we could use the relation between d and h to find a non-uniform triangulation adapted to the given data. But in this paper, we assume that we have no a priori knowledge of the dataset distribution. Therefore, a uniform triangulation seems to be the best tradeoff (see [START_REF] Arcangéli | Multidimensional minimizing splines. Theory and applications[END_REF]).

Experimental results.

We conclude the paper by presenting several numerical experiments on both synthetic datasets and real data. The experiments have been performed on a 2.21-GHz Athlon with 1.00 GB of RAM.

In this part, we give a discretization of the variational problem [START_REF] Clément | Approximation by Finite Element Functions Using Local Regularisation[END_REF]. In the following, Ω is a rectangle, and we limit ourselves to the case n = 2 and m = 2 so that the Sobolev's embedding

H 2 (Ω, R 2 ) C 0 ( Ω, R 2 ) holds.
We keep the notations introduced in Sections 2 and 3. We solve (10) on a finite element space (V h )

2 ⊂ H 2 Ω, R 2 whose generic finite element is the C 1 rectangular Bogner-Fox-Schmit (BFS) element. If (K, P K , Σ K ) is the generic C 1 BFS element, the domain Ω is triangulated by means of the rectangles K. It can be proved that problem [START_REF] Clément | Approximation by Finite Element Functions Using Local Regularisation[END_REF] uncouples with respect to each component in R 2 .

Let (v q ) q=1,2 be the components of v ∈ H 2 Ω, R 2 and (w q i ) q=1,2 be the components of [START_REF] Clément | Approximation by Finite Element Functions Using Local Regularisation[END_REF] can therefore be stated by:

w i ∈ R 2 , ∀i ∈ {1, • • • , N }. Problem
                     Search for u ǫ = (u q ǫ ) q=1,2 ∈ H 2 (Ω, R 2 ) such that ∀v = (v q ) q=1,2 ∈ H 2 (Ω, R 2 ), ∀q = 1, 2, N i=1 u q ǫ (a i )v q (a i ) + ǫ(u q ǫ , v q ) 2,Ω,R = N i=1 w q i v q (a i ). (26)
We solve (26) in V h for q = 1, 2. Let M h be the dimension of V h and P h j j=1,••• ,M h be the basis functions. If we denote by u h,q ǫ the solution of problem (26) approximated in (V h ) 2 , we can decompose u h,q ǫ into the basis

P h j j=1,••• ,M h :                ∀q = 1, 2, ∃ α q j j=1,••• ,M h ∈ R, u h,q ǫ = M h j=1 α q j P h j . (27)
For q = 1, 2, taking successively

v q = P h l , l = 1, • • • , M h in (26), the studied problem is rephrased by:                Search for α q ∈ R M h such that, ∀l = 1, • • • , M h , N i=1 M h j=1 α q j P h j (a i )P h l (a i ) + ǫ M h j=1 α q j P h j , P h l 2,Ω,R = N i=1 w q i P h l (a i ). (28)
Denoting by A h and R h the matrices defined by:

A h = P h j (a i ) 1≤i≤N, 1≤j≤M h ∈ M N ×M h (R), R h = P h j , P h i 2,Ω,R 1≤i≤M h , 1≤j≤M h ∈ M M h ×M h (R),
we are led to solve two uncoupled sparse linear subsystems that can be written in the form:

A h T A h + ǫR h α q = A h T w q , ∀q = 1, 2, (29)
with w q ∈ R N , ∀q = 1, 2. It can be easily demonstrated that these systems admit a unique solution.

The matrix that intervenes in both systems is symmetric, positive definite, and independent of q.

Only the second member differs in both systems.

In each test, the parameter ǫ has been set to 10 -6 . Also, in order to provide a quantitative measure of the degree of approximation, we have computed for each example the relative error χ defined as follows:

χ =       N i=1 u h ǫ (a i ) -w i 2 2 N i=1 w i 2 2       1 2 . ( 30
)
We now propose several examples to illustrate the efficiency of the proposed method.

Experimental tests on synthetic datasets

In the first experiment, we are concerned with the function f defined on Ω = (0, 1) × (0, 1) by plots of each component of the vector field, Fig. 2. The finite element space V h is constructed on the triangulation from the Bogner-Fox-Schmit rectangle of class C 1 and the corresponding approximant u h ǫ related to A is computed. The results are presented in Fig. 34. We have used a mesh refinement for the visualization: each square has been subdivided into 10 × 10 equal squares. As shown, the approximant proves to be satisfactory: the relative error is equal to 0.0003946990. To measure the accuracy of the method, we have then successively removed (randomly) 25%, 50%, 75% and 85% of the data and have computed the relative error only on the missing data. The obtained results are displayed in Table 1. The linear regression associated with the data X = log(d) and Y = log(error) gives as slope: 3.42. This coefficient is very close to the coefficient given by Theorem 3.6 : 3.5.

f (x, y) = exp -(3x -1) 2 -(3y -
In the second experiment, we focus on Franke's function defined on Ω = (0, 1) × (0, 1) by: f (x, y) =0.75 exp -0.25 (9x -2) Also, this figure is complemented by the plots of each component of the vector field, Fig. 6. The finite element space V h is constructed on the triangulation from the Bogner-Fox-Schmit rectangle of class C 1 and the corresponding approximant u h ǫ related to A is computed. The results are presented in Fig. 78. We have used a mesh refinement for the visualization: each square has been subdivided into 20 × 20 equal squares. The relative error is equal to 0.0035175393. To measure the accuracy of the method, we have then successively removed (randomly) 25%, 50%, 75% and 85% of the data and have computed the relative error only on the missing data. The obtained results are displayed in Table 2.

Datasets of wind velocity fields

We focus on data of wind velocity fields that have been acquired in the Mozambique Channel. The Mozambique Channel is a portion of the Indian Ocean located between the island of Madagascar and Southeast Africa (see Fig. 9). These data have been measured by the satellite ERS2. They have been smoothed and monthly averaged by the CERSAT laboratory (IFREMER-France).

No data have been recorded on the land so we have set them to 0. In the first experiment, the data were acquired in January. In January, the monsoon is present and induces a north wind in the north part of the channel. The south-east wind present in the south part of the channel corresponds to trade winds (see [START_REF] Monier | Dynamic of the currents in the Mozambique canal[END_REF] for more details). The data are regularly distributed on the domain Ω. The set A consists of 528 points. The triangulation of the domain Ω is made up of 15 × 15 rectangles and so M h = 4 × 16 2 = 1024. The given set of vectors {w i } i=1,••• ,N is depicted in Fig. 10. Also, this figure is complemented by the plots of each component of the vector field, Fig. 11. The finite element space V h is constructed on the triangulation from the Bogner-Fox-Schmit rectangle of class C 1 and the corresponding approximant u h ǫ is computed. The results are presented in Fig. 12-13. We have used a mesh refinement for the visualization: each rectangle has been subdivided into 20 × 20 rectangles. The relative error is equal to 0.0001650361 which is considered as a very good result in this context. We removed 25% of the data (randomly) and computed the relative error at all data locations. It reached 0.076156. For higher percentages, as the number of data is rather small, we do not get as good results as for the synthetic data: the relative error is in 10 -1 .

The second example concentrates upon the same kind of data. The processed data were acquired in August. In August, the monsoon does not exist and the wind field is dominated by trade winds. The set A consists of 528 points. The triangulation of the domain Ω is made up of 15×15 rectangles and so M h = 4 × 16 2 = 1024. The given set of vectors {w i } i=1,••• ,N is depicted in Fig. 14. Also, this figure is complemented by the plots of each component of the vector field, Fig. 15. The finite element space V h is constructed on the triangulation from the Bogner-Fox-Schmit rectangle of class

1 ) 2 .

 12 To apply our method, we have generated a fixed set of points A = {a 1 , • • • , a N } regularly distributed on Ω and have set: ∀i = 1, • • • , N , w i = ∇f (a i ). The set A consists of 1681 points. The triangulation of the domain Ω is made up of 8 × 8 equal squares and so M h = 4 × 9 2 = 324. The fixed set of points A is depicted on Fig. 1, left column, and the given vector field {w i } i=1,••• ,N on Fig. 1, right column. The reference frame is centred on the bottom left corner. The x-direction corresponds to the column direction while the y-direction corresponds to the row direction and points upwards. Also, this figure is complemented by the

Figure 1 :

 1 Figure 1: On the left, the regularly distributed fixed set of points A. On the right, the given vector field. Note that the Matlab function used to display the field scales the arrows to fit within the grid.

Figure 2 :

 2 Figure 2: On the left, plot of the first component of the vector field {w i } i=1,••• ,N . On the right, plot of the second component.

Figure 3 :

 3 Figure 3: On the left, plot of the first component of the obtained vector field. On the right, plot of the second component.

Figure 4 :

 4 Figure 4: Plot of the obtained vector field.

2 - 2 .

 22 0.25 (9y -2) 2 + 0.75 exp -(9x + 1) 2 /49 -(9y + 1) 2 /10 +0.5 exp -0.25 (9x -7) 2 -0.25 (9y -3) 2 -0.2 exp -(9x -4) 2 -(9y -7) As previously done, we have generated a fixed set of points A = {a 1 , • • • , a N } regularly distributed on Ω and have set ∀i = 1, • • • , N , w i = ∇f (a i ). The set A consists of 1681 points. The triangulation of the domain Ω is made up of 8 × 8 equal squares and so M h = 4 × 9 2 = 324. The fixed set of points A is depicted on Fig. 5, left column, and the given vector field {w i } i=1,••• ,N on Fig. 5, right column. The reference frame is centred on the bottom left corner. The x-direction corresponds to the column direction while the y-direction corresponds to the row direction and points upwards.

Figure 5 :

 5 Figure 5: On the left, the regularly distributed fixed set of points A. On the right, the given vector field. Note that the Matlab function used to display the field scales the arrows to fit within the grid.

Figure 6 :

 6 Figure 6: On the left, plot of the first component of the vector field {w i } i=1,••• ,N . On the right, plot of the second component.

Figure 7 :

 7 Figure 7: On the left, plot of the first component of the obtained vector field. On the right, plot of the second component.

Figure 8 :

 8 Figure 8: Plot of the obtained vector field every 3 rows and columns.

Figure 9 :

 9 Figure 9: Depiction of the Mozambique Channel.

Figure 10 :

 10 Figure 10: Plot of the given vector field {w i } i=1,••• ,N .

Figure 11 :

 11 Figure 11: On the left, plot of the first component of the vector field {w i } i=1,••• ,N . On the right, plot of the second component.

Figure 12 :

 12 Figure 12: On the left, plot of the first component of the obtained vector field. On the right, plot of the second component.

Figure 13 :

 13 Figure 13: Plot of the obtained vector field (every six rows and columns). The blue diamond indicates the point to which the vector is attached.

Figure 14 :

 14 Figure 14: Plot of the given vector field {w i } i=1,••• ,N .

Figure 15 :

 15 Figure 15: On the left, plot of the first component of the vector field {w i } i=1,••• ,N . On the right, plot of the second component.

Figure 16 :

 16 Figure 16: On the left, plot of the first component of the obtained vector field. On the right, plot of the second component.

Figure 17 :

 17 Figure 17: Plot of the obtained vector field (every six rows and columns). The blue diamond indicates the point to which the vector is attached.

  endowed with the usual norm denoted by || • || m,Ω,R n and the associated semi-norm denoted by | • | m,Ω,R n ).

	Proof .
	Let us take f such that ||f || A,m,Ω,R n = 0. It implies that |f | m,Ω,R n = 0 and taking into account the connectedness of Ω, it yields f ∈ P m-1 (Ω, R n ). As the set A contains a P m-1 -unisolvent subset, we deduce that f ≡ 0. It is now clear that || • || A,m,Ω,R n is a norm on H m (Ω, R n ) associated with a scalar product. We now prove the equivalence of the norm || • || A,m,Ω,R n with the norm || • || m,Ω,R n . First, we have:

Table 1 :

 1 Relative error when some data are randomly removed. (Fixed triangulation). now propose to illustrate numerically the theoretical convergence result provided in Theorem 3.5. To do that, we consider ε(d) = d 5 . The theoretical error related to fu d

		Test 1 on synthetic dataset
	Percentage of remaining data Relative error
		75		0.000471
		50		0.000590
		25		0.002217
		15		0.004729
	We ǫh 0,Ω,R 2	should be
	of order O(d 7/2 ).
	d	ε	Error	% of data
	0.0442	1.68 10-7	9.8281 10-4	100
	0.0630 9.913 10-7 8.7284 10-4	75
	0.0718 1.9115 10-6	0.0014	50
	0.0954 7.887 10-6	0.0099	25
	0.1063 1.3548 10-5	0.0165	15
	0.1978 3.0276 10-4	0.1496	5
	0.3815	0.0081	0.6248	1

Table 2 :

 2 Relative error when some data are randomly removed. (Fixed triangulation).