Florent Foucaud
email: foucaud@labri.fr

Tero Laihonen

Aline Parreau
email: aline.parreau@ujf-grenoble.fr

An improved lower bound for (1, ≤ 2)-identifying codes in the king grid *

Keywords: Identifying codes, king grid, density, discharging method

We call a subset C of vertices of a graph G a (1, ≤ ℓ)-identifying code if for all subsets X of vertices with size at most ℓ, the sets {c ∈ C|∃u ∈ X, d(u, c) ≤ 1} are distinct. The concept of identifying codes was introduced in 1998 by Karpovsky, Chakrabarty and Levitin. Identifying codes have been studied in various grids. In particular, it has been shown that there exists a (1, ≤ 2)-identifying code in the king grid with density 3 7 and that there are no such identifying codes with density smaller than 5 12 . Using a suitable frame and a discharging procedure, we improve the lower bound by showing that any (1, ≤ 2)-identifying code of the king grid has density at least 47 111 .

Introduction

Let G = (V, E) be a simple undirected graph with vertex set V and edge set E. The closed neighbourhood of a vertex v ∈ V , which consists of the vertex itself and all the adjacent vertices, is denoted by N [v]. In addition, we write N [X] = v∈X N [v] for X ⊆ V . In [START_REF] Honkala | Codes for identification in the king lattice[END_REF], a subset C ⊆ V is called a (1, ≤ ℓ)-identifying code if

N [X] ∩ C = N [Y ] ∩ C
for any two distinct subsets X ⊆ V and Y ⊆ V of size at most ℓ (one of them can be the empty set). The elements of a code C are called codewords.

Identifying codes were introduced by Karpovsky, Chakrabarty and Levitin in [START_REF] Karpovsky | On a new class of codes for identifying vertices in graphs[END_REF] and can be applied to locate objects in sensor networks [START_REF] Ray | Robust location detection in emergency sensor networks[END_REF]. A network is modelled by a graph and a sensor can check its closed neighbourhood. It gives an alarm if it detects at least one of the sought objects there. Suppose we have a (1, ≤ ℓ)-identifying code C in the graph, and we place the sensors to the vertices corresponding to the codewords of C. Then, knowing the set of alarming sensors A ⊆ C, we can determine where the objects are (assuming that there are at most ℓ of them). Indeed, we have A = N [X] ∩ C for some subset X of size at most ℓ, and because these are unique, we can determine X -the set of vertices where the objects are.

Of course, we would like to use as few sensors as possible, so our aim is to find identifying codes with smallest possible cardinality. For infinite grids, we define below the measure density for this purpose.

Identifying codes have been considered, for example, in the following infinite grids: the square grid, the triangular grid, the king grid and the hexagonal mesh [START_REF] Ben-Haim | Exact minimum density of codes identifying vertices in the square grid[END_REF][START_REF] Charon | Identifying codes with small radius in some infinite regular graphs[END_REF][START_REF] Cohen | Improved identifying codes for the grid[END_REF][START_REF] Cohen | Bounds for codes identifying vertices in the hexagonal grid[END_REF][START_REF] Cohen | On codes identifying vertices in the twodimensional square lattice with diagonals[END_REF][START_REF] Cranston | A new lower bound on the density of vertex identifying codes for the infinite hexagonal grid[END_REF][START_REF] Honkala | On identifying codes in the triangular and square grids[END_REF][START_REF] Honkala | On identifying codes in the hexagonal mesh[END_REF][START_REF] Honkala | On identification in the triangular grid[END_REF]. For more papers concerning the topic of identification, see [12]. In this paper, we focus on (1, ≤ 2)-identifying codes in the king grid. The king grid has vertex set Z 2 and two vertices are adjacent if and only if the Euclidean distance between them is at most √ 2. Hence, the closed neighbourhood of a vertex consists of nine vertices (see Figure 1).

Figure 1: The infinite king grid with a (1, ≤ 2)-identifying code of density 3 7 .

In [START_REF] Honkala | Codes for identification in the king lattice[END_REF], a construction of a (1, ≤ 2)-identifying code in the infinite king grid with density 3 7 is given (see Figure 1). Recently, in [START_REF] Pelto | New bounds for (r, ≤ 2)-identifying codes in the infinite king grid[END_REF], a lower bound of 5 12 was obtained. In this paper, we show that our approach which utilizes certain frames, improves this lower bound to 47 111 .

Preliminaries

We call code any subset of Z 2 . As usual, the density D(C) of a code C ⊆ Z 2 is defined by:

D(C) = lim sup n→∞ |C ∩ Q n | |Q n |
where

Q n = {(x, y) ∈ Z 2 | |x| ≤ n, |y| ≤ n}.
The next theorem gives the previously known bounds for the optimal density of a (1, ≤ 2)identifying code in the king grid.

Theorem 1 ( [START_REF] Cohen | On codes identifying vertices in the twodimensional square lattice with diagonals[END_REF][START_REF] Pelto | New bounds for (r, ≤ 2)-identifying codes in the infinite king grid[END_REF]). Let d be the optimal density of a (1, ≤ 2)-identifying code in the king grid. Then 0.416 = 5 12 ≤ d ≤ 3 7 = 0.4286.... To prove bounds on the density of a code C, one method is to use a finite subset X of Z 2 . For any vertex v ∈ Z 2 , we denote by v + X or X + v the set {v + u | u ∈ X}, and by n(X, C) the average number of codewords of C in the translations of X: n(X, C) = lim sup n→∞ v∈Qn |(v+X)∩C| |Qn|

. If we have some knowledge about n(X, C), we can get results on the density of C with the following proposition.

Proposition 2. Let C be a code of Z 2 and X a nonempty finite subset of Z 2 , then:

D(C) = n(X, C) |X| . Proof. Let h be an integer such that X ⊆ Q h . Let n be a positive integer. A vertex of (v + X) ∩ C, for v ∈ Q n+h is either in C ∩ Q n , or in Q n+2h \ Q n . Among all the sets (v + X) ∩ C, with v ∈ Q n+h , each vertex of C ∩ Q n is counted exactly |X| times and each vertex of Q n+2h \ Q n is counted at most |X| times.
Hence we have:

|X||C ∩ Q n | ≤ v∈Q n+h |(v + X) ∩ C| ≤ |X|(|C ∩ Q n | + |Q n+2h \ Q n |).
By dividing every term by |Q n |, we obtain:

|X| • |C ∩ Q n | |Q n | ≤ v∈Q n+h |(v + X) ∩ C| |Q n | ≤ |X| • |C ∩ Q n | |Q n | + |X| • |Q n+2h \ Q n | |Q n | .
Since h is a fixed integer, we have: lim sup n→∞

|Q n+2h \Qn| |Qn| = 0. As n → ∞ we obtain: lim sup n→∞ v∈Q n+h |(v + X) ∩ C| |Q n | = |X| • D(C).
To end the proof, one can easily verify that lim sup

n→∞ v∈Q n+h |(v+X)∩C| |Qn| = n(X, C).
The identifying code problem can be seen as a covering problem. Indeed, a code C is a (1, ≤ ℓ)identifying code if and only if, for any two distinct subsets of vertices, X and Y , of size at most ℓ, the symmetric difference N [X]∆N [Y ] is covered by at least one element of C. The following theorem, from [START_REF] Cohen | On codes identifying vertices in the twodimensional square lattice with diagonals[END_REF], reduces (in the case of the king grid) the family of sets that a code has to cover to be a (1, ≤ 2)-identifying code.

Theorem 3 ([5], Theorem 2.1). A code C is a (1, ≤ 2
)-identifying code in the king grid if and only if, for all (x, y) ∈ Z 2 , the sets 1. {(x, y), (x + 3, y), (x, y + 3), (x + 3, y + 3)}, 2. {(x, y), (x + 1, y), (x + 2, y)} 3. {(x, y), (x, y + 1), (x, y + 2)} each contain at least one codeword of C.

Note that Theorem 3 together with Proposition 2 directly implies that the density of a (1, ≤ 2)identifying code of the king grid is at least 1 3 . To improve this lower bound (and the one of Theorem 1), we will use, as the set X of Proposition 2, the set of Figure 2 (without loss of generality, one can set the left-lower corner of this set to be the origin). We call frame any set isomorphic to the set of Figure 2. The aim of this paper is to give a lower bound on the average number of codewords in a frame. More precisely, we will use a discharging procedure to show that the average number of codewords in a frame is at least 5 + 3 37 . Hence, we will need to consider the neighbourhood of a frame. There is a natural bijection between all the frames and the set Z 2 . We will also consider the frames lattice over the set of all frames, where the distance between two frames will be the distance in the king grid between the two corresponding vertices. As an example, the 2-ball of a frame F is the set of frames

{F + (x, y) | (x, y) ∈ Z 2 , max{|x|, |y|} ≤ 2}.
We will often study a frame together with all the vertices at distance at most two of some vertex of the frame. Therefore, we will use the notation Xy with X∈A-H and y∈a-h, for a vertex on line X and row y with respect to the coordinates of Figure 2. The four vertices in positions Cc, Cf, Fc and Ff are called the corners of the frame. A side of a frame is a set of four vertices of the frame lying in the same line or column.

The conditions of Theorem 3 can be reformulated using our terminology. A code C is a (1, ≤ 2)identifying code in the king grid, if and only if for each frame F : In the following, we assume that C is a (1, ≤ 2)-identifying code of the king grid. Given an integer k, a k-frame (resp. k + -frame) is a frame containing exactly k (resp. at least k) vertices of C. Lemma 4. Let F be a frame, then F is a 4 + -frame. If F is a 4-frame, then all its codewords are in the corners of F . Proof. It is easily observed that any frame F can be partitioned into four sets of three consecutive vertices on a line or column. Hence, by Condition 2, F contains at least four codewords.

Condition 1: At least one corner of F is a codeword of C.
Assume that F contains exactly four codewords. By Condition 1, one of the codewords must be in a corner, say in position Cc. Then, by Condition 2, there must be another codeword on line C. If it was on column d or e, one could partition the remaining vertices into three sets of three consecutive vertices on a line or column and F would be a 5 + -frame, a contradiction. So there is a codeword in position Cf. Using similar arguments, the two other codewords must be in positions Fc and Ff. Observation 5. If a frame has three non-corner codewords on two of its sides, then it is a 6 +frame.

Proof. Assume first that the three non-corner codewords lie on opposite sides. Without loss of generality, we can assume there are in positions Cd, Ce and Fd. By Condition 1, there must be one codeword in a corner. Assume it is in column c. Then, by Condition 2, there must be another codeword in column c and one in column f, leading to a 6 + -frame. If the corner codeword is in column f, then there would be two codewords in column f and one in column c, leading again to a 6 + -frame.

We assume now that the three non-corner codewords are on adjacent sides. Without loss of generality, we can assume they are not in the set S composed by line C and column f. By Condition 2, there are at least two codewords in S. By Condition 1, there is a codeword in a corner. If there is a codeword in position Fc then F is a 6 + -frame. Otherwise, there is a corner codeword in S and then there must be three codewords in S, leading to a 6 + -frame. Observation 6. If a frame F has a corner codeword c and two non-corner codewords on one of the sides of F that c does not belong to, then F is a 6 + -frame.

Proof. Without loss of generality, we assume that there are three codewords in positions Cc, Fd and Fe. By Condition 2, there must be another codeword in column c, another codeword in line C and one among positions Df, Ef and Ff. Hence, F is a 6 + -frame. Lemma 7. Let F be a 4-frame. Then the four frames F + {(0, 2), (0, -2), (2, 0), (-2, 0)} are 6 + -frames.

Proof. By Lemma 4, the four codewords of F are its four corners. By Condition 1 applied on F + (1, 0), there is a codeword on either position Cg or position Fg. Hence by Observation 5, F + (2, 0) is a 6 + -frame. The claim is obtained by symmetry.

Lemma 7 gives a short proof of a result first given in [START_REF] Pelto | New bounds for (r, ≤ 2)-identifying codes in the infinite king grid[END_REF]. The proof is a warm-up for our main result.

Theorem 8 ([13]

). The density of any (1, ≤ 2)-identifying code of the king grind is at least 5 12 .

Proof. Let C be a (1, ≤ 2)-identifying code of the king grid. We use the discharging method to show that the average number of codewords in a frame is at least 5. Since the size of each frame is 12, the result will follow by Proposition 2.

In the beginning, each k-frame has charge k. Then, each 6 + -frame F gives charge 1 4 to each 4-frame among frames F +{(0, 2), (0, -2), (2, 0), (-2, 0)}. By Lemma 4, there exist only 4 + -frames. By Lemma 7, each 4-frame receives charge 1 and each 6 + -frame gives at most charge 1 away. So after the discharging process, each frame has at least charge 5. Each frame only gives charge to vertices at distance at most 2, hence the average number of codewords of C in a frame is at least 5 and we are done.

In the following, we will prove our main result: Theorem 9. The density of any (1, ≤ 2)-identifying code of the king grid is at least 47 111 = 0.423. To prove this theorem, we will use a similar technique than the one of the proof of Theorem 8, but giving a final charge of more than 5 to each frame. Hence, a charge of 5 can be seen as the reference value for the charge of a frame. We say that the charge excess of a k-frame is k -5. The charge excess within a subset S of frames is the sum of the charge excesses of all the frames of S.

Structural properties of code C

We now prove some results on the structure of C in the viewpoint of the frames. We call 4benefactor a 6 + -frame F having a 4-frame among frames F + {(0, 2), (0, -2), (2, 0), (-2, 0)}.

Lemma 10. Let F be a 4-benefactor. Then F has 6 + -frames in each corner position of its 2-ball being at distance 2 from the 4-frames among F + {(0, 2), (0, -2), (2, 0), (-2, 0)}. Moreover, if F is a 6-frame, there is a unique 4-frame among frames F + {(0, 2), (0, -2), (2, 0), (-2, 0)}. If F is a 7 + -frame, there are at most two 4-frames among frames F + {(0, 2), (0, -2), (2, 0), (-2, 0)}.

Proof. The first part of the claim is a direct consequence of Lemma 7.

There is at least one 4-frame among F +{(0, 2), (0, -2), (2, 0), (-2, 0)}, without loss of generality, we can assume that F + (-2, 0) is a 4-frame. Then F + (2, 0) cannot be a 4-frame since then, no corner of F is a codeword, contradicting Condition 1 on F . Suppose now that there is another 4-frame among F + {(0, 2), (0, -2)}, say F + (0, 2). As before, F + (0, -2) cannot be a 4-frame. By Condition 1 applied on F , there is a codeword in position Ff. By Condition 1 applied on F + (-1, 0) and F + (0, 1), there are two codewords among positions Ce, Fe, Ec and Ef. Hence F is a 7 + -frame and we are done. The 6 + -frames in the corner positions of a 4-benefactor 6 + -frame described in the previous lemma are called the co-benefactors of F . The next lemma is valid for other cases by smmetry. Lemma 11. Let F be a 4-benefactor 6-frame oriented as in Figure 3. Then:

F 6 + 6 + 4 X X X X X X X X X Z X Z X Z Z Z Y Y X X X
• Frame F has at least two 6 + -frames F 1 , F 2 in its 2-ball (in addition to its co-benefactors).

• If F has no 6 + -frame in a Y -or Z-position, then F has at least three 6 + -frames in its 2-ball (in addition to its co-benefactors), or one of its co-benefactors is a 7 + -frame and F + (-1, 0) is not a 6 + -frame.

Proof. In this proof, we will often use Conditions 1 and 2 without explicitly referring to them. By Condition 1 applied on F + (-1, 0), there is a codeword in either Ce or Fe. Without loss of generality, we may assume that there is a codeword in position Ce (the other case follows by symmetry).

Since F is not a 7 + -frame, there is no codeword in position Fe and only one codeword among positions Dc and Ec.

Assume there is a codeword in position Ec and no codeword in position Dc. Due to Condition 2, there is a codeword in De. Hence, by Observation 5, F + (-1, 1) is a 6 + -frame. If F + (-1, 0) is a 6 + -frame, we are done. Otherwise, F + (2, 1) is a 6 + -frame and we are also done.

Assume now that there is no codeword in position Ec. Then by Condition 2, there is a codeword in position Dc.

If there is a codeword in position De, F +(-1, 0) is a 6 + -frame and, by Observation 6, F +(-1, 1) is a 6 + -frame and we are done. Hence we may suppose that there is no codeword in position De. This implies that there is a codeword in position Df, and no codeword in position Ef.

If there is a codeword in position Cf, by Observation 5, F + (0, 1) is a 6 + -frame. Moreover, by Observation 6, F + (1, 0) is a 6 + -frame and we are done. So, we may assume there is no codeword in position Cf, which implies that there is a codeword in position Ff.

If there is a codeword in position Gd, by Observation 5, F + (-1, -1) and F + (1, -2) are 6 + -frames and we are done. Hence we may assume that there is no codeword in position Gd.

If there is a codeword in position Fg, F + (1, 0) is a 6 + -frame and by Observation 5, F + (2, 0) is a 6 + -frame and we are done. So we may assume that there is no codeword in position Fg.

If there is a codeword in position Db, F +(-1, 0) is a 6-frame and by Observation 5, F +(-2, -1) is 6 + -frame. If F + (-1, -1) is a 6 + -frame we are done, otherwise, by Observation 5, F + (0, -2) is a 6 + -frame. Hence we may assume that there is no codeword in position Db.

If there is no codeword in position Bf, then there is a codeword in position Bc, and by Observation 5, F + (-2, 1) is a 6 + -frame. Then, there must be a codeword among positions Bd and Be, implying that F + (-1, 1) is a 6 + -frame. Hence we may assume that there is a codeword in position Bf.

If there is a codeword in position Be, by Observation 5, F + (1, 1) and F + (2, 2) are 6 + -frames and we are done. Hence we assume that there is no codeword in position Be.

If there is no codeword in position Eg, by Observation 5, F + (2, -1) and F + (2, 2) are 6 +frames. Now, if there is a codeword in Bd, then by Observation 5, F + (1, 2) is a 6 + -frame and we are done. Otherwise, there is a codeword in Bc and Ba, and F + (-2, 1) is a 6-frame, and we are done. So, assume that there is a codeword in position Eg.

If there is a codeword in position Dg, then F + (1, 0) is a 6 + -frame and by Observation 5, F + (2, -1) is a 6 + -frame. So we assume there is no codeword in position Dg, which implies that there is a codeword in position Gg by Condition 1.

If there is a codeword in position Gf, by Observation 5, F + (0, -2) and F + (2, -1) are 6 +frames. Now, if there is a codeword in position Ge, F + (1, -1) is a 6 + -frame and we are done. Otherwise, there is a codeword in position Gb and F + (-2, -1) is a 6 + -frame. So we assume that there is no codeword in position Gf, which implies that there is a codeword in position Ge.

If there is a codeword in position Gb, F + (-1, -1) and F + (-2, -1) are 6-frames and we are done. So we assume that there is no codeword in position Gb.

If there is a codeword in position Bd, by Observation 5, F + (1, 2) is a 6 + -frame. Moreover F + (1, 1) is a 6 + -frame so we are done. So we assume there is no codeword in position Bd. Then there is a codeword in both positions Ba and Bc, implying that F + (-2, 1) is a 6-frame.

If there is a codeword in position Cg, F + (1, 0) is a 6-frame and we are done. So, we assume there is no codeword in Cg. Then, by Observation 5, F + (2, 1) is a 6-frame. Now, if F + (1, -2) is a 6 + -frame we are done. Otherwise, there is no codeword in position Hd and exactly one codeword among positions He and Hf. Since F + (0, -2) is a 5 + -frame, there is a codeword in position Hc. By Condition 1, there is a codeword in position Ha. This implies that F + (-2, -2) is a 7-frame and finishes the case analysis.

A frame is called 1-poor (resp. 2-poor) if it is a 5-frame having at most one 6-frame in its 2-ball (resp. no 6 + -frame at distance 1 and at most two 6-frames at distance 2). Lemma 12. Let F be a 5-frame, then F has at least one 6 + -frame in its 2-ball. Moreover, one of the following properties holds:

• F has a 6-frame at distance 1 and another 6 + -frame in its 2-ball,

• F has a total charge excess of at least 3 within its 2-ball, with at least two 6 + -frames,

• F is 1-poor and, up to symmetry, the configuration around F is depicted on Figure 5, • F is 2-poor and, up to symmetry, the configuration around F is depicted on Figure 6. Proof. Let F be 5-frame. We make a case analysis according to the configurations of F depicted on Figure 4. Note that cases A 1 to A 1 0 represent cases where F has only one corner codeword, cases B 1 to B 4 are those cases where F has exactly two corner codewords on the same side, cases C 1 to C 3 are the cases where F has three corner codewords, and case D is the case where F has four corner codewords. The cases which are symmetric to those of Figure 4 will follow from the same arguments. Moreover, note that if F had exactly two corner codewords on opposite corners, by Condition 2, each side of F would contain an additional codeword and F would be a 6 + -frame.

A 1 A 2 A 3 A 4 A 5 A 6 A 7 A 8 A 9 A 10 B 1 B 2 B 3 B 4 C 1 C 2 C 3 D
Case A 1 . By Observation 5, F + (-1, 0) and F + (0, 1) both are 6 + -frames and we are done.

Case A 2 . By Observation 5, F + (-1, 0) is a 6 + -frame. By Observation 6, F + (0, 1) is a 6 + -frame and we are done.

Case A 3 . By Observation 6, F + (-1, 0) and F + (0, 1) both are 6 + -frames and we are done.

Case A 4 . By Observation 5, F + (0, 1) is a 6 + -frame. Now, by Condition 2, there is a codeword on either position Dd or Ed. In either case, by Observation 5, F + (-1, 2) or F + (1, -1) respectively, is a 6 + -frame and we are done.

Case A 5 . By Observation 6, F + (0, 1) is a 6 + -frame. Now, by Condition 2, there is a codeword on either position Dd or De. In either case, by Observation 5, F + (-1, -1) or F + (2, 1) respectively, is a 6 + -frame and we are done.

Case A 6 . By Observation 6, F + (0, 1) is a 6 + -frame. Now, by Condition 2, there is a codeword in either position Dd or Ed. In both cases, by Observation 5, F + (-1, -1) or F + (1, -1) respectively, is a 6 + -frame and we are done.

Case A 7 . By Observation 5, F +(0, 1) is a 6 + -frame. If F +(-1, -1) is a 6 + -frame, we are done. Otherwise, from the assumptions there are already four non-corner codewords in F + (-1, -1), hence F + (-1, -1) is a 5-frame and, by Condition 1, all its further non-corner positions contain no codeword. In particular, there is no codeword in both positions Dd and Gd. By Condition 1, there is a codeword on at least one of the positions Dg and Gg. Hence, by Observation 5, F + (2, -1) is a 6 + -frame and we are done.

Case A 8 . Assume there is no codeword in position Dd. Then, F + (-1, 0) and F + (2, 1) are 6 + -frames and we are done. So, assume there is a codeword in position Dd.

If there is a codeword in position Ed, by Observation 5, F +(-2, 0) and F +(1, 0) are 6 + -frames. Hence we may assume there is no codeword in position Ed, and by symmetry, no codeword in position De either.

If there is a codeword in position Ee, by Observation 5, F + (1, 1) and F + (-1, -1) are 6 + -frames. Hence we may assume there is no codeword in position Ee.

If there is a codeword in position Be, by Observation 5, F + (1, 1) and F + (2, 2) are 6 + -frames. Hence we may assume there is no codeword in position Be, and by symmetry, no codeword in position Eb either. By Condition 1, this implies that there is a codeword in position Bb.

If there is a codeword in position Eg, by Observation 5, F + (1, 0) and F + (2, 1) are 6 + -frames. Hence we may assume there is no codeword in position Eg, and by symmetry, no codeword in position Ge either.

If there is no codeword in position Bd, by Observation 5, F + (0, 2) and F + (-2, 1) are 6 +frames. If F + (1, 1) is a 6 + -frame, we are done. Otherwise, by Observation 5, F + (2, 1) is a 6 + -frame and we are done too. Hence we may assume there is a codeword in position Bd, and by symmetry, another codeword in position Db.

If there is a codeword in position Bc, F + (0, 1) is a 6-frame and by Observation 5, F + (-1, 1) is a 6 + -frame, so we are done. Hence we may assume there is no codeword in position Bc, and by symmetry, no codeword in position Cb either.

If there is a codeword in position Gb, F + (-1, -1) and F + (-1, -2) are 6 + -frames. Hence we may assume there is no codeword in position Gb, and by symmetry, no codeword in position Bg either.

If there is no codeword in position Gf, there is, by Condition 1, a codeword in position Gc and then F + (-1, -1) is a 6-frame. Now, if there is no codeword in position Fg, F + (1, -1) is a 4-frame and by Lemma 7, F + (1, 1) is a 6 + -frame and we are done. Otherwise, by Observation 5, F + (1, -2) is a 6 + -frame. Hence we may assume that there is a codeword in position Gf, and by symmetry, another codeword in position Fg.

If there is a codeword in position Gg, F +(1, -1) is a 6-frame and by Observation 5, F +(1, -2) is a 6 + -frame, so we are done. Hence we may assume there is no codeword in position Gg.

If there is a codeword in position Gc, F + (0, -1) and F + (-1, -1) is a 6 + -frame, so we are done. Hence we may assume there is no codeword in positions Gc and Cg.

If F + (-1, 2) is a 5-frame, by Condition 2, there is no codeword in position Ae and a codeword on exactly one of the positions Ac and Ad. Then there must be a codeword in position Af because F + (0, 2) has a fifth codeword. But then F + (2, 2) is a 6 + -frame. Hence, there is a 6 + -frame among F + (-1, 2) and F + (2, 2). By symmetry, there is also a 6 + -frame among F + (-2, 1) and F + (-2, -2).

If F + (-1, -2) is a 6 + -frame, we are done. Otherwise, there is no codeword in position Hd, which implies that there is a codeword in position Hg, as well as one codeword among positions He and Hf. By symmetry, if F + (2, 1) is a 5-frame, there a codeword in position Gh, as well as a codeword among positions Eh and Fh. Then, F + (2, -2) is a 6 + -frame and we are done.

Case A 9 . By Condition 2, there are at least two codewords among positions Dd, De, Ed and Ee. If there are two codewords on the same column or line, (resp. d, D, e, E), then by Observation 5, resp. F + (1, 0), F + (0, -1), F + (-1, 0), F + (0, 1), is a 6 + -frame. Hence, if there are three codewords among positions Dd, De, Ed and Ee, there are two 6 + -frames at distance 1 of F and we are done.

Therefore, we assume that there are exactly two codewords among positions Dd, De, Ed and Ee. Then by Condition 2, they cannot lie on a same line or column. If there are codewords in positions De and Ed, by Observation 5, F + (1, -1) and F + (-1, 1) are 6 + -frames and we are done. Hence, we assume that there are codewords in positions Dd and Ee.

If there is a codeword in position Eb, by Observation 5, F + (-2, 1) is a 6 + -frame, and by Condition 2, F + (-1, 1) is a 6 + -frame as well. Hence we may assume there is no codeword in position Eb.

If there is a codeword in position Cb, by Observation 5, F + (-2, 0) is a 6 + -frame, and by Condition 2, F + (-1, 0) is a 6 + -frame as well. Hence we may assume there is no codeword in position Cb.

If there is a codeword in position Gd, F + (1, -1) is necessarily a 6 + -frame, and by Observation 5, F + (1, -2) is a 6 + -frame. Hence we may assume there is no codeword in position Gd.

By symmetry, we can assume that there is no codeword in any of the positions Be, Bc and Dg. Using Condition 2, this implies that there are codewords in positions Db and Bd.

If there is a codeword in position Gg, F + (1, -1) is necessarily a 6 + -frame, and by Observation 5, F + (2, -1) is a 6 + -frame. Hence we may assume there is no codeword in position Gg.

If there is a codeword in position Ge, F + (1, -1) is necessarily a 6 + -frame. Now, if there is a codeword in position Fb, F +(-1, -1) is necessarily a 6-frame too. Otherwise, either F +(-2, -1), F + (-1, 1) or F + (-2, 2) is a 6 + -frame, and we are done. Hence we may assume there is no codeword in position Ge and by symmetry, there is no codeword in position Eg either.

If there is a codeword in position Fb (resp. in position Bf), F + (-1, 0) (resp. F + (0, 1)) is necessarily a 6-frame. If there is a codeword in both positions, we are done. Hence by symmetry, without loss of generality, we can assume that there is no codeword in position Fb. This implies that F + (-2, -1) is a 6 + -frame. If there is a 6 + -frame at distance 1 of F , we are done. Otherwise, F + (1, 2) is a 6 + -frame. If F + (-2, 2) is a 6 + -frame, we are done. Otherwise, F + (-2, -1) and F + (1, 2) are both necessarily 7 + -frames and we are done too. If there is a codeword in position De or Ee, then F + (-1, 0) is a 6 + -frame. If there is a codeword in position De (resp. Ee), by Observation 5, F + (1, -1) (resp. F + (2, -1)) is a 6 + -frame and we are done.

Case C 1 . By Observation 5, both F + (-1, 0) and F + (0, 1) are 6 + -frames.

Case C 2 . By Observation 5, F + (0, 1) is a 6 + -frame. By Condition 1, there is a codeword among positions Bc and Bf. Then, by Observation 5, F + (0, 2) is a 6 + -frame.

Case C 3 . Assume first that there is no codeword in position De. Then, there must be codewords in both positions Dd and Ee. Then, by Observation 5, F + (1, 1) and F + (-1, -1) are both 6 + -frames. Hence we can assume that there is a codeword in position De, and, by symmetry in position Ed.

If there is a codeword in position Dd, by Observation 5, F +(0, 2) and F +(0, -1) are 6 + -frames and we are done. Hence we may assume there is no codeword in position Dd.

If there is a codeword in position Ee, by Observation 5, F + (0, 1) and F + (2, 0) are 6 + -frames and we are done. Hence we may assume there is no codeword in position Ee.

If there is a codeword in position Bc or Bf, by Observation 5, F + (0, 2) is a 6 + -frame and F + (0, 1) must be a 6 + -frame too, so we are done. Hence we may assume there is no codeword in positions Bc and Bf and, by symmetry, on positions Cb and Fb.

If there is a codeword in position Be, then F + (0, 1) is a 6 + -frame. Now, if either F + (-1, 1) or F + (-2, 2) is a 6 + -frame, we are done. Otherwise, there is no codeword in position Ca and hence there is a codeword in position Fa by Condition 1.This implies that F + (-2, -1) is a 6 + -frame, so we are done. Hence we may assume there is no codeword in position Be and, by symmetry, in position Eb.

If there is a codeword in position Eg, by Observation 5, F +(2, 1) and F +(2, -2) are 6 + -frames. If F + (1, 0) is a 6 + -frame, we are done. Otherwise, there is no codeword in position Dg, and a codeword among positions Cg and Fg. If there is a codeword in position Fg, by Observation 5, F + (1, -1) is a 6 + -frame. If there is a codeword in position Cg, by Observation 5, F + (2, 0) is a 6 + -frame. Hence we may assume there is no codeword in position Eg and, by symmetry, in position Ge.

If there is a codeword in position Gc, by Observation 5, F + (0, -2) and F + (-2, -1) are 6 + -frames. If there are only two 6-frames and no additional 6 + -frame in the 2-ball of F , F is 2-poor, then, using the same techniques than previously, one can check that the codewords around F are fixed as in Figure 6 and we are done. Otherwise, there is a total charge excess of at least 3 within the 2-ball of F and we are done too. Hence we may assume there is no codeword in position Gc and, by symmetry, in position Cg.

We note that F + (1, -1) is a 6 + -frame. If F + (1, -1) is a 6-frame and if it is the only 6-frame within the 2-ball of F , then, one can check that F is 1-poor and that the codewords around F are fixed as in Figure 5, and so, we are done. Otherwise, either F + (1, -1) is a 7 + -frame and then F + (1, -2) is a 6 + -frame, or F + (1, -1) is a 6-frame and there exists another 6 + -frame within the 2-ball of F . In both cases we are done.

Case D. By Observation 5, F + (0, 1) is a 6 + -frame. By Condition 1, there is a codeword among positions Bc and Bf. Then, by Observation 5, F + (0, 2) is a 6 + -frame. Lemma 13. If F is a 1-poor frame, then F has a 6-frame F 6 at distance 1 such that:

• F 6 has at most two 1-poor frames in its 2-ball and each of them is at distance 1 of F 6 , • F 6 has at least two 6 + -frames at distance 2,

• F 6 is not a 4-benefactor.

Proof. By Lemma 12, if F is a 1-poor frame, the neighbourhood of F is fixed and depicted in Figure 5 (up to symmetry). The only 6 + -frame in the neighbourhood of F is F 6 = F + (1, -1) and contains exactly six codewords. The first point of the lemma follows from the fact that a 1-poor frame must have a 6-frame in a corner of its 1-ball, however by using the left part of Figure 5 the positions marked by "X" on the right part of the figure cannot host a 1-poor frame. The second point follows using Observation 5 on the left part of Figure 5. The last point follows by noting that a 4-benefactor must have two non-codeword corners in the same column or line, which is not the case here.

A frame playing the role of F 6 in the previous lemma is called a 1-poor-benefactor. Lemma 14. If F is a 2-poor frame, then F has two 6-frames at distance 2. One of them, say F 6 , has the following properties:

• F 6 is neither a 4-benefactor nor a 1-poor-benefactor,

• F 6 has at least four 6 + frames in its 2-ball, 4-benefactor. Without loss of generality, we can assume that F 6 is oriented as in Figure 3. Then F 6 gives always charge 2α to the 5-frames in Z-positions, except in Rule 1a where F 6 + (0, 1) receives only charge α, but in this case, F 6 + (0, 1) has three 6 + -frames in its 2-ball, so F = F 6 + (0, 1). Frame F 6 + (-1, 0) has three 6 + -frames in its 2-ball, so F = F 6 + (-1, 0). Thus, we can assume that F is on a Y -position. Hence, F ′ 6 is also a 4-benefactor. In this case, by Rules 1 and 2, both F 6 and F ′ 6 give charge 3α 2 to F , and we are done. Lemma 18. After the application of the discharging rules, each 4-frame has charge at least 5 + 3α.

Proof. By Lemma 7, a 4-frame F has four 4-benefactors. By Rules 1, 2 and 3, each 4-benefactor gives charge 3α+1 4 to F . Hence, F receives charge 1 + 3α and ends with charge 5 + 3α.

After the application of our discharging rules, by Lemmas 15, 16, 17 and 18, each frame has charge at least 5 + 3α = 188 37 , therefore, the average number of codewords in each frame is at least 188 37 (this is due to the fact that each frame gives charge to vertices at distance at most 2). There are twelve vertices in each frame, hence, by Proposition 2, we obtain Theorem 9.

Condition 2 :

 2 Each set of three consecutive vertices on a side of F contains at least one codeword of C.

Figure 2 :

 2 Figure 2: Notation for vertices in the neighbourhood of a frame.

Figure 3 : 2 -

 32 Figure 3: 2-ball around a 4-benefactor 6-frame.

Figure 4 :

 4 Figure 4: All the possibilities for a 5-frame.

Case A 10 .Case B 4 .

 104 Assume that there are no codewords among positions Ed and Ee. Then, by Condition 2, there are codewords on both positions Dd and De, and by Observation 5, F + (0, -1) and F + (2, 1) are 6 + -frames. Now, assume that there is a codeword in position Ed or Ee. Then, by Condition 2, F + (0, 1) is a 6 + -frame. If there is a codeword in position Ed, by Observation 5, F + (1, -1) is a 6 + -frame. If there is a codeword in position Ee, by Observation 5, F + (1, 1) is a 6 + -frame. In both cases we are done.Case B 1 . By Observation 5, F + (-1, 0) is a 6 + -frame. By Condition 1, there is a codeword among positions Cb and Fb. Then, by Observation 5, F + (-2, 0) is a 6 + -frame and we are done.Case B 2 . By Observation 5, F + (-1, 0) is a 6 + -frame. By Condition 2, there is a codeword among positions Dd and Ed. If there is a codeword in position Dd (resp. Ed and not in Dd), by Observation 5, F + (1, 1) (resp. F + (1, -1)) is a 6 + -frame and we are done.Case B 3 . By Observation 5, F + (-1, 0) is a 6 + -frame. If F + (1, 1) is a 6 + -frame, we are done. Otherwise, F + (1, 1) is a 5-frame and it has a unique corner codeword in position Ed. By Condition 2, this implies that there is a codeword in position Bc, and then by Observation 5, F + (0, 2) is a 6 + -frame and we are done. Assume first that there are no codewords among positions De and Ee. Then by Condition 2, there must a codeword on both positions Dd and Ed. Then by Observation 5, F + (1, 0) and F + (-2, 0) are 6 + -frames.

Figure 5 :

 5 Figure 5: Neighbourhood of a 1-poor frame in the vertex lattice (on the left) and in the frame lattice (on the right).

Figure 6 :

 6 Figure 6: Neighbourhood of a 2-poor frame in the vertex lattice (on the left) and in the frame lattice (on the right).

This research is supported by the ANR project IDEA, under contract ANR-08-EMER-007, 2009-2011.

• F 6 has only one 2-poor frame in its 2-ball.

Proof. By Lemma 12, if F is a 2-poor frame, the neighbourhood of F is fixed and depicted in Figure 6 (up to symmetry). Let F 6 = F + (-2, -1). The first point of the lemma follows from the fact that, by Lemma 12, a 1-poor benefactor must have a 6-frame in its 2-ball (and only this 6-frame in its 2-ball) and that a 4-benefactor must have two non-codeword corners in the same line or column (hence, the only possibility would be that F + (-2, 1) is a 4-frame which is not true). The second point follows by using Observation 5 on the left part of Figure 6. For the third point, observe that the only possible locations for another 2-poor frame in the 2-ball of F 6 are those marked by an "X" on the right part of Figure 6: F + (-3, 1) and F + (-1, 1). However, both have two 6 + -frames in their 2-ball that do not match the codeword configuration of a 2-poor frame as given on the left part of Figure 6.

A frame playing the role of F 6 in the previous lemma is called a 2-poor-benefactor.

The discharging procedure

We are now ready to describe the discharging procedure which leads to our result. In the beginning, each k-frame has charge k. Let α = 1 37 . We apply the following rules:

1. Let F be a 4-benefactor 6-frame without any 7 + -frame among its co-benefactors. Assume F is oriented as in Figure 3 (other cases follow by symmetry). We consider two subrules:

(a) If there is a 6 + -frame in a Y -position of Figure 3, say in position F + (-1, 1) (the other case will be covered by symmetry), and no other 6 + -frames in Y -and Z-positions of Figure 3, then F gives charge 3α+1 4 to the 4-frame F + (-2, 0), charge α to frame F + (0, 1), charge 2α to the other 5-frames in Z-positions of Figure 3, charge 3α 2 to the other 5-frame in Y -position of Figure 3 2. Let F be a 4-benefactor 6-frame with a 7 + -frame among its co-benefactors. Assume F is oriented as in Figure 3, and that F + (-2, 2) is a 7 + -frame (other cases follow by symmetry).

Then F gives the same charges than in Rule 1b, except for the frame F + (-1, 1) which receives no charge from F .

3. A 4-benefactor 7 + -frame F gives charge 3α+1 4 to the 4-frames among F +{(0, 2), (0, -2), (2, 0), (-2, 0)}, charge 3α to frames at distance 1 and charge 2α to frames at distance 2.

4. A 1-poor-benefactor, recall that it is a 6-frame, gives charge 3α to the 1-poor frames in its 2-ball, 2α to the other 5-frames at distance 1 and charge α to the other 5-frames at distance 2.

5. A 2-poor-benefactor, recall that it is a 6-frame, gives charge 2α to the unique 2-poor frame in its 2-ball, 2α to the other 5-frames at distance 1 and charge α to the other 5-frames at distance 2.

6. Other 6-frames give charge 2α to 5-frames at distance 1 and charge α to 5-frames at distance 2.

7. Other 7 + -frames give charge 3α to 5-frames at distance 1 and charge 2α to 5-frames at distance 2.

We note that the rules are not ambiguous. Indeed, by Lemmas 13 and 14, a 6-frame can be either a 4-benefactor, a 1-poor benefactor, a 2-poor benefactor or not a benefactor at all, but never two at the same time.

Lemma 15. After the application of the discharging rules, each 7 + -frame has charge at least 5 + 3α.

Proof. Let F be a 7 + -frame. It is sufficient to show that F gives at most charge 2 -3α = 71 37 . If F is a 4-benefactor, then by Lemma 10, it has at most two 4-frames among the frames F + {(0, 2), (0, -2), (2, 0), (-2, 0)}.

If F has two 4-frames among F + {(0, 2), (0, -2), (2, 0), (-2, 0)} then it has at least three cobenefactors. Hence, F has at most eleven 5-frames at distance 2. By Rule 3, F gives at most charge 2

37 . If F has only one 4-frame among F + {(0, 2), (0, -2), (2, 0), (-2, 0)}, then it has two cobenefactors and hence, at most thirteen 5-frames at distance 2. By Rule 3, F gives at most charge 3α+1 4 + 8 • 3α + 13 • 2α = 60 37 . Finally, if F is not a 4-benefactor, then by Rule 7, it gives at most charge 8•3α+16•2α = 56 37 . Lemma 16. After the application of the discharging rules, each 6-frame has charge at least 5 + 3α.

Proof. Let F be a 6-frame. It is sufficient to show that F gives at most charge 1 -3α = 34 37 . Assume first that F is a 4-benefactor without any 7 + -frame among its co-benefactors. If there were no other 6 + -frames than the co-benefactors of F within its 2-ball, F would give, according to Rule 1, charge 3α+1 4

2 + 14α = 37 37 . However, by Lemma 11, F has at least two non-co-benefactor 6 + -frames in its 2-ball. Moreover by Lemma 11 as well, either one of those is on a Y -or Z-position, or there is a third non-co-benefactor 6 + -frame in the 2-ball of F . We distinguish two cases. If there is an extra 6 + -frame on a Y -position and one on an X-position (Rule 1a applies), then we save charge at least 3α 2 + α from the extra frames and charge α from frame F + (0, 1). Otherwise (Rule 1b applies), the total charge saved on the extra 6 + -frames is at least 3α. Indeed, we save at least 2 • 3α 2 if the two extra 6 + -frames are both on Y -positions, 2 • 2α if they are both on Z-positions, 2α + 3α 2 for a Z-position and a Y -position, 2α + α for a Z-position and an X-position, and 3 • α if there are three extra 6 + -frames. In all cases, F at most 34 37 of charge.

If F is a 4-benefactor with a 7 + -frame among its co-benefactors, by Lemma 11, there are two 6 + -frames in the 2-ball of F in addition to its two co-benefactors and they are at distance 2 from F . Hence, F has at most eleven 5-frames at distance 2. By Rule 2, F gives at most 3α+1 4

+ 5 • 2α + 3α 2 + 12α = 33.5 37 . If F is a 1-poor-benefactor, by Lemma 13, F has at most two 1-poor frames in its 2-ball and each of them is at distance 1 of F . Moreover, there are at least two 6 + -frames at distance 2 of F . We may assume that F has two 1-poor frames in its 2-ball since by Rule 4, F would give away more charge in this case. Now, by Rule 4, F gives at most 2 • 3α + 6 • 2α + 14 • α = 32 37 . If F is a 2-poor-benefactor, by Lemma 14, F has only one 2-poor frame and at least four 6 + -frames in its 2-ball. We may assume that these frames are at distance 2 of F since by Rule 5, F would give away more charge in this case. By Rule 5, F gives at most 2α + 8 • 2α + 11 • α = 29 37 . Finally, if F is not a benefactor, then by Rule 6, it gives at most charge 8 • 2α + 16 • α = 32 37 . Lemma 17. After the application of the discharging rules, each 5-frame has charge at least 5 + 3α.

Proof. It is enough to prove that each 5-frame receives charge 3α.

We first note that by our discharging rules, each 6-frame F 6 gives at least charge α to each 5-frame of its 2-ball, except in Rule 2 where F 6 does not give anything to one 5-frame F ′ at distance 1 of F 6 . However, F ′ has a 4-benefactor 7 + -frame at distance 1, which, by Rule 3, gives charge 3α to F ′ .

Let F be a 5-frame. By the previous paragraph, we can consider that F receives at least charge α from each 6-frame in its 2-ball. By Rules 3 and 7, F receives at least charge 2α from each 7 + -frame in its 2-ball. Hence, if F has a total charge excess of at least 3 within its 2-ball, with two 6 + -frames, F receives at least charge 3α and we are done. Otherwise, by Lemma 12, F is either 1-poor, 2-poor or has a 6-frame at distance 1 and another one in its 2-ball.

If F is 1-poor, by Lemma 13, it has a 1-poor-benefactor in its 2-ball which, by Rule 4, gives charge 3α to F .

If F is 2-poor, by Lemma 14, it has a 2-poor-benefactor in its 2-ball which, by Rule 5, gives charge 2α to F . Moreover, by Lemma 14, F has another 6 + -frame in its 2-ball which gives charge α to F .

Finally, suppose F has a 6-frame F 6 at distance 1 and another one, F ′ 6 , in its 2-ball. If F receives charge 2α from F 6 , we are done. Otherwise, by our discharging rules F 6 is necessarily a


F De Dinechin

J.-M Muller

B Pasca

A Plesco

An FPGA architecture for solving the Table Maker's Dilemma

Keywords: table maker's dilemma, floating-point arithmetic, correct rounding, elementary functions, FPGA

Solving the Table Maker's Dilemma, for a given function and a given target floating-point format, requires testing the value of the function, with high precision, at a very large number of consecutive values. We give an algorithm that allows for performing such computations on a very regular architecture, and present an FPGA implementation of that algorithm.

I. THE TABLE MAKER'S DILEMMA

The IEEE 754-2008 standard for Floating-Point arithmetic [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF] recommends (yet does not dictate) that some elementary functions should be correctly rounded. That is, given a rounding function •, (e.g., round to nearest even, or round to ±∞), when evaluating function f at the floatingpoint number x, the system should always return •(f (x)).

Building a fast correctly rounded library for some target floating-point (FP) format requires preliminarily solving a problem called the table maker's dilemma. This requires very large computations which may use environments and formats totally different from the target environment and format. In this paper, we suggest performing these computations on an FPGA.

For the sake of simplicity, in this paper, we assume that the target rounding function is the IEEE 754 standard double-precision round to nearest even. On the target environment, to compute f (x) in a given format, where x is a FP number, we must first compute an approximation to f (x) with a given accuracy [START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF], which we round to the nearest FP number in the considered target format. The problem is the following: find what the accuracy of the approximation must be to make sure that the obtained result be always equal to the "exact" f (x) rounded to the nearest FP number. To solve that problem we have to locate, for each considered target floating-point format and for each considered function f , the hardest to round (HR) points. Defining a midpoint as the exact middle of two consecutive floating-point numbers, an HR point is an input floating-point numbers x such that f (x) is very close to a midpoint, without being exactly a midpoint.

This work is partly supported by the TaMaDi project of the french Agence Nationale de la Recherche (ANR). The donation of a DK-DEV-4GX530N board by the Altera University Program is also gratefully acknowledged.

Assuming (after some re-normalization) that x and f (x) are between 1 and 2, that the target floating-point format is a binary format of precision p, we need to find the largest possible value of m such that there exist a FP number x that satisfies:

• f (x) is not exactly equal to a midpoint;

• the binary representation of f (x) has the form m bits

1.xxxxx • • • xxx p bits 1000000 • • • 000000 xxx • • • or m bits 1.xxxxx • • • xxx p bits 0111111 • • • 111111 xxx • • • ;
Two different algorithms have been suggested for dealing with this problem. The first one, first presented in Lefèvre PhD dissertation [START_REF] Lefèvre | Moyens arithmétiques pour un calcul fiable[END_REF], [START_REF]New results on the distance between a segment and Z 2 . Application to the exact rounding[END_REF], allowed Lefèvre and Muller to publish the first tables of HR points for the most common functions in double-precision/binary64 FP arithmetic [START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF]. We will call it the L-algorithm. The second one, the SLZ algorithm, was introduced by Stehlé, Lefèvre and Zimmermann [START_REF] Stehlé | Worst cases and lattice reduction[END_REF], [START_REF]Searching worst cases of a one-variable function[END_REF]. It was used to find the HR points for the exponential function in decimal64 arithmetic [START_REF] Lefèvre | Worst cases for the exponential function in the IEEE 754r decimal64 format[END_REF]. SLZ has a better asymptotic complexity than the L-algorithm, however, when the target format is the double precision/binary64 format of the IEEE 754 standard, they require similar computational times: weeks of computation for all input exponents (hours of computation for one input exponent), using massive parallelism.

The problem with these algorithms does not only lie in this huge computation time. In both cases we have a very complex algorithm, implemented in a very complex program, that runs for weeks and just outputs one result: what confidence can we have in that result? The HR points are the main weak point of a library such as CRLibm [START_REF] Daramy-Loirat | CR-LIBM, a library of correctly-rounded elementary functions in double-precision[END_REF]: each function of that library comes with a theorem of the form "if the HR points have been rightly computed then the function always outputs a correctly rounded result". Hence, our major goal here is to design a very simple, very regular algorithm (therefore suited for FPGA implementation): if it outputs the same results as the L-algorithm, then this will give much confidence in these results.

The method we are going to suggest here has an even worse asymptotic complexity than the L-and SLZalgorithms. And yet, due to its simplicity, the hidden constant in the complexity term is so small that, still with double precision/binary64 as a target, our method will require similar times on a single FPGA.

Let us now present our method. Defining u = 2 1-p , we will compute the values of f (1) mod

2 1-p , f (1 + u) mod 2 1-p , f (1 + 2u) mod 2 1-p , f (1 + 3u) mod 2 1-p , . . . , f (2) mod 2 1-p ,
with a given, predetermined, accuracy 2 -µ (with µ larger than p and-for probabilistic reasons [START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF]-less than 2p). Each time we find a value f (1 + ku) mod 2 1-p extremely close to 2 -p (i.e., whose leading bits are of the form 011111111 • • • or 10000000 • • • ), we output the value of k for further testing using multipleprecision software. The major difficulty here is that there are 2 p-1 values f (1 + ku) mod 2 1-p and p is fairly large (a typical value is 53 for the double precision/binary64 format). Therefore, the computation of these values must be done very quickly.

To do this, we will approximate f by some polynomial P (with an accuracy of approximation significantly better than 2 -µ ), and compute the successive values P (1 + ku) mod 2 1-p using a modulo 2 1-p adaptation of the wellknown tabulated differences method [START_REF] Knuth | The Art of Computer Programming[END_REF].

II. THE TABULATED DIFFERENCES METHOD

Let P be a polynomial of degree n and x 0 a real. We define x k = x 0 + ku for k > 0, and we wish to compute the successive values P (x 0 ), P (x 1 ), P (x 2 ), P (x 3 ), . . . . The tabulated differences method is based on the fact that if we define the following "discrete derivatives":

• P (1) (x) = P (x + u) -P (x); • P (2) (x) = P (1) (x + u) -P (1) (x); • . . . • P (n) (x) = P (n-1) (x + u) -P (n-1) (x); then P (n) is a constant C.
This leads to the following algorithm.

Initialization: compute P (x 0 ), P (x 1 ), P (x 2 ), . . . , P (x n ), and deduce from these values all the possible partial discrete derivatives, of which we keep the initial vector at point x 0 : P (n) (x 0 ) = C, P (n-1) (x 1 ), P (n-2) (x 2 ), . . . , P (2) (x n-2 ), P (1) (x n-1 ), P (x n ).

Iteration:

The following recurrence computes the value of this vector at point x k+1 = x k + u out of the vector at 2) (x k+n-1 ) = P (2) (x k+n-2 ) + P (3) (x k+n-2 ) P (1) (x k+n ) = P (1) (x k+n-1 ) + P (2) (x k+n-1 ) P (x k+n+1 ) = P (x k+n ) + P (1)

point x k .                  P (n-1) (x k+2 ) = P (n-1) (x k+1 ) + C P (n-2) (x k+3 ) = P (n-2) (x k+2 ) + P (n-1) (x k+2 ) . . . . . . . . . P (
(x k+n ) (1)
P (x2) P (xn-1) P (xn) P (x1)

P (1) (x0)
P (1) (xn-1) P (1) (x1)

P (x0) P (2) (x1) P (n-1) (x0) P (n) (x0) = C C C P (1) (xn) P (xn+1)
P (n-1) (x3)
Iteration

P (n-1) (x1) P (n-1) (x2)
Initialization Figure 1. The tabulated difference method

The computations in (1) are done modulo 2 1-p : they are simple fixed-point additions, with the bits of weight ≥ 2 1-p being ignored. Notice that the n additions in (1) are straightforwardly pipelined, hence, once the initialization is done computing a new value P (x k ) takes the time of one addition.

It should also be noted that the problem is embarrassingly parallel: Although the iteration itself is intrinsically sequential, the full domain of an elementary function in doubleprecision may be split into arbitrarily many sub-domains, and we may perform initialization/iteration processes in parallel for each sub-domain. We indeed aim at processing hundreds of sub-intervals in parallel within a single FPGA.

Besides, the initial values P (x 0 ), P (x 0 + u), P (x 0 + 2u), . . . , P (x 0 + nu) cannot be computed exactly in practice. They are correct within some rounding error, and we will see in the following that when performing (1), these errors accumulate quite quickly. After some value of k, say k max this accumulated error becomes unacceptable, and we have to invoke the initialization process again, with x 0 replaced by x 0 + k max u. In other words, the size of a sub-interval is dictated by an error analysis that will be the subject of Section III.

The initialization process requires n + 1 polynomial evaluations, and n(n+1)/2 subtractions. There are two possible ways of performing it: 1) on the FPGA itself, or 2) in software on the host computer. In both cases we must choose k max so that the initialization time is totally overlapped by the iterations [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF].

The initialization process first involves evaluating the polynomial in n + 1 points using a classical multiplicationbased scheme (typically Horner's). Modern FPGAs contain up to several thousand small multipliers that could be used for this purpose, but designing an architecture for this initialization would add a lot to the FPGA design effort. In the sequel of the paper, we therefore choose the simpler second approach. It also has the advantage of exploiting the computing power of the host processor. However, there is a price to pay: due to limited bandwidth between the host and the FPGA, iteration (1) must run for a much longer k max . We will see in next section that this entails a significantly wider data-path, hence more resource consumption for the iteration hardware, possibly cancelling the benefits of saving the initialization hardware. This question remains to study quantitatively.

Let us now formalize the dependency between k max and the datapath width.

III. ERROR ANALYSIS

Let us bound the difference between the value computed F (x k ) and the true value of the function f (x k ). We may first decompose this error as follows:

F (x k ) -f (x k ) = (F (x k ) -P (x k )) + (P (x k ) -f (x ))
The second term is the approximation error, and the Remez approximation algorithm will allow us to keep it as small as needed. Let us focus on the first term, the rounding error.

Let δ(i)(x 0 + ku) be the overall rounding error in the initial evaluation of P (i) (x 0 + ku). Notice that the additions in (1) are performed in fixed-point, ignoring the outgoing carries: they do not induce any error, yet they propagate the initial rounding errors on the P (i) (x 0 ). Let be a bound on the errors on P (n) (x 0 ) = C, P (n-1) (x 0 + u), P (n-2) (x 0 + 2u), . . . , P (1) (x 0 + (n -1)u). We easily find

δ (n-1) (x 0 + ku) ≤ δ (n-1) (x 0 + (k -1)u) + , so that δ (n-1) (x 0 + ku) ≤ (k + 1) .
From that, we deduce

δ (n-2) (x 0 + ku) ≤ δ (n-2) (x 0 + (k -1)u) +δ (n-1) (x 0 + (k -1)u) ≤ δ (n-2) (x 0 + (k -1)u) + k , so that δ (n-2) (x 0 + ku) ≤ (1 + 2 + 3 + • • • + k) = k(k + 1) 2 .
Similarly,

δ (n-3) (x 0 + ku) ≤ 1 + 2(2+1) 2 + 3(3+1) 2 + • • • + (k-1)k 2 = (k-1)(k)(k+1)
6

.

An elementary induction shows that the bound on the error of the computed value of P (x 0 + ku) satisfies

δ (0) (x 0 + ku) ≤ (k -n + 2) • • • (k -1)(k)(k + 1) n! . (2)
IV. A CASE STUDY: THE EXPONENTIAL FUNCTION All parameters in the method are function-dependent, so we cannot give a general performance result (although it should not vary much with the function). Hence, we give here some figures related to the exponential function on the input interval [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF][START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF]. We take f (x) = 1 2 e x to normalise the output to [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF][START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF].

We first split the input interval into 2 m sub-intervals, each of size 2 -m , and we will compute one approximation polynomial on each sub interval, using the Remez algorithm. The trade-off here is between the degree of the obtained polynomials (the smaller, the better for the subsequent evaluation) and the number of Remez polynomial to compute. A good choice here is m = 15: on each of the 2 15 = 32768 sub-intervals, a polynomial of degree 4 approximates f (x) with an accuracy better than 2 -90 , and the Sollya tool is able to compute all these polynomial and formally validate their accuracy [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF] in about 3 hours. Now we must choose k max , the length of an evaluation run between reinitializations. Having decided that reinitializations are performed on the host processor, we now have to take into account the limits on 1/ computing power of the host processor and 2/ data bandwidth between processor and FPGA. A larger k max means fewer initializations, but larger data-path, hence slower operation and less parallelism. Note also that if we have P parallel iteration cores, the host must serve them all.

Our current trade-off is to take k max = 2 20 . Equation (2) with degree n = 4 tells us that the error is smaller than 2 76 . For our target accuracy of 2 -85 modulo 2 52 , we need to have 85 -52 = 33 valid bits at the end of the computation. A datapath width of 33 + 76 = 109 bits ensures this accuracy. Assuming M = 2 8 parallel iterations on the FPGA, the host must be able to compute one initialization every 2 20 /2 8 = 4096 FPGA cycles. FPGA cycles are typically 10 times slower than processor cycles, so the host has roughly 40,000 cycles to compute each initialization. Efficient multiple-precision libraries such as GMP and MPFR make this possible. Host-FPGA bandwidth, in this scenario, is not a problem.

V. OUR DESIGN

Field-programmable gate arrays (FPGAs) are a natural technology for implementing this type of algorithm:

• for a given set of input parameters we need to perform a big, one-off computation. Once completed we can reconfigure the FPGA for a different set of input parameters.

• the implemented method is based on binary additions for which FPGAs are very efficient.

Nevertheless, implementing such a complex, multiparametrized design using a hardware description language (HDL) is a tedious and error-prone task. To overcome this, we used the FloPoCo framework [START_REF]FloPoCo, Floating-Point Core generator[END_REF]. We implemented in FloPoCo an architecture generator that inputs the many parameters and generates an FPGA-specific implementation of the circuit in the form of a portable, human-readable VHDL file.

A. Functional model 1) TaMaDi Core: The core component of our design is the polynomial evaluator based on the tabulated differences method. The architecture of this component is depicted in Figure 2. Its main entities are the n -1 adders chained together which are used to evaluate the vector of discrete derivatives.

Each computation starts with the component receiving a '1' value on the Initialize line together with a unique interval identifier on the Interval bus. During the next N + 1 clock cycles, the values of the initialization vector C = P (n) (x 0 ), ..., P (x 0 + nu) are received in sequence on the DataIn bus, and fill the pipeline. A counter is used to keep track of k in evaluating P (x 0 + ku). The output of the n th adder feeds a pattern detector unit, implemented as wide ANDs. A value of '1' at the output of the pattern detector signals that the value present on the output Counter bus, together with the interval identifier, points to one HR case. The component raises the Ready line to '1' when it has finished the allowed number of iterations and needs a new reinitialization.

The architecture of the TaMaDi Core (Core) is perfectly suited to FPGA hardware. Adders benefit from the fast carrychains which allow the simple ripple carry adder (RCA) scheme to be implemented efficiently. The wide AND of the pattern detector may also take advantage of these fast carry-chains.

Table I presents area and timing post place-and-route results of the TaMaDi Core on modern FPGAs from Xilinx [START_REF]Virtex-5 FPGA User Guide, Xilinx[END_REF], [START_REF]Virtex-6 FPGA User Guide[END_REF] and Altera [START_REF]Stratix IV Device Handbook[END_REF], [START_REF] Stratix | [END_REF] 2) TaMaDi Cluster: Multiple TaMaDi Cores may be assembled in a larger component named the TaMaDi Cluster, whose architecture is depicted in Figure 3.

The presented system has several parameters:

• M the number of TaMaDi Cores in the system. The maximum value of this parameter depends on the size of one Core and the size of the FPGA. The practical value for this parameter also depends on the bandwidth between FPGA and the host system, processing element datapath width and degree and the reinitialization interval.

• size of the input and output FIFOs, also depends on bandwidth, M and processing element characteristics. • size of processing element output FIFO. Their dimension can be as small as one element. However, for good performance their size depends on the probability of finding HR cases in that interval and the output bandwidth. The TaMaDi Cluster is connected to the host system (or the next hierarchical level) my means of two FIFOs. Data is fed by the host system to ClusterInFIFO while this FIFO is not full. Each element on this FIFO has the structure depicted below:

The ClusterInFIFO element contains the necessary infor- Core is ready to process new information (signaled '1' on the corresponding output Ready port) the input FIFO is popped one element. The uppermost w iid bits of information containing the interval ID are fed to the processing element together with a value on '1' on the corresponding Initialize input pin. The lowermost (N + 1)p w bits are loaded into a N + 1-level shift-register in order be serialized in chunks of p w bits. During the next N +1 clock cycles, the shift-register feeds the TaMaDi Core as the pipeline starts.

P (n-1) (x0 + u) P (n-2) (x0 + 2u) P (X0 + nu)
When the TaMaDi Cores signals the detection of a HR case, the information concerning this case (counter value and interval identifier, totaling w k + w iid bits) is pushed into the corresponding Core output FIFO.

The data from the CoreFIFOs is then placed in the ClusterOutFIFO whenever this FIFO is not full. Simple priority encoders on both inputs and outputs manage the access to the Cluster input and output FIFOs.

3) TaMaDi System: The TaMaDi Cluster has low resource count and fast clock speeds for modest number of Cores (up to 32). Unfortunately, this multiplexed data dispatch architecture scales badly due to several issues: (1) size of the multiplexers and priority encoder-decoder circuitry [START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF] the long lines between the dispatcher (shiftregister in our implementation) and the computing cores (TaMaDi Cores). By grouping multiple TaMaDi Cores into a small number of Clusters (where the size of the Cluster remains reasonably small say 16), we could potentially use the same dispatching architecture at a macro-level.

However, when filling up a large FPGA chip, a new problem arises: connections between the dispatcher and the Clusters become very long, introducing large delays. To overcome this, we have used a different, credit-based, dispatcher, depicted in Figure 5. It allows us to pipeline the communication lines to Clusters, thus breaking the long delays into several shorter ones. To keep track of data in flight on these long lines, a sender has a credits counter, initialized to the size of the receiver FIFO. The sender sends only when its credit counter is positive, decrements it on each send, and blocks when is is zero. The receiver sends a credit back to the sender at each read from its FIFO.

A detailed view of this sender/receiver interface is presented in Figure 6. The receiver part of the interface is tightly coupled to the TaMaDi Cluster.

A global dispatcher circuit manages the communication with the host. It has two functionalities [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF] dispatching reinitialization data to the TaMaDi Cluster (Figure 5 left) and (2) retrieval of HR cases from these modules (Figure 5 right).

4) Full Prototype: For prototyping purposes we use the Stratix IV GX development board featuring a Stratix IV GX EP4SGX530KH40C2 FPGA. The board communicates with a host PC by means of a PCI Express 2.0 8x interface that can provide up to 3.4 GB/s full-duplex.

The Altera PCI Express hard IP together with the PLDA EZDMA2 IP [START_REF]PLDA EZDMA2 DMA for PCI Express Hard IP[END_REF] offered in the PLDA reference design ensure a simple FIFO interface for our pipelined creditbased Dispatcher Interface (Figure 7). The PLDA host driver offers a high-level API interface for feeding and retrieval of information from the DMA FIFOs by means of multiple DMA channels (2 in our case).

B. Bandwidth requirement

In this section we will compute the bandwidth requirement of the entire TaMaDi System. First, we need to compute the bandwidth of TaMaDi Cluster depicted in Figure 3.

We use the following notations: k max the number of iterations between re-initializations w k = log 2 (k max ) , width in bits of the iteration counters. η the maximum number of intervals to be processed by the system w iid = log 2 (η) , width in bits of interval identifiers ξ the probability of finding one HR case The input bandwidth for one TaMaDi Core:

f
B in Core = ((N + 1)w dp + w iid )

f k max + N + 1 (3)
and the output bandwidth of one Core is:

B out Core = ξ( log 2 (k max ) + w iid ) f k max + N + 1 (4)
The Core bandwidth is B Core = B in Core +B out Core . Considering that a TaMaDi Cluster has M such Cores, the total bandwidth requirement for a Cluster is B Cluster = M • B Core . A TaMaDi System is composed out of K Clusters, therefore requiring a bandwidth equal to:

B System = K • B Cluster = K • M • B Core .
Table II presents the dependency between the parameters of a Core, its area and the required bandwidth for keeping it busy at 100MHz. A larger bandwidth requirement leads to more pressure on the I/Os but a smaller Core size, which allows fitting more in one single FPGA.

For a system comprising of 100 TaMaDi Cores, each requiring a bandwidth of 5.42 MBit/s the bandwidth requirement is approximatively 0.5 Gb/s which seems to be reasonably within our available bandwidth potential. Nevertheless, this configuration would require us generating more than 57 TBytes of reinitialization data (some of which can indeed be generated on-the-fly) compared to a more manageable 1.9 TBytes required for a w dp = 120.

For a 200-Core system B System = 3.4M bits/s, which can easily be provided by FPGA platforms connected through to the host system through the PCIE bus, Ethernet interface and even USB2.0.

C. Performance estimation

We are currently using the Altera StratixIV development kit based on an EP4SGX530KH40C2 FPGA to prototype our system. This gives us an environment for estimating the performance and scalability of our architecture. However, as presented in this section the amount of computation needed for an elementary function, on one exponent value is as the order of tens of hours. We therefore envision mapping this architecture on even larger FPGAs, and even multi-FPGA based systems.

In this section we provide a performance estimation for the case of the exponential function for double precision p = 52 (we consider one input exponent). Considering the 2 20 iterations until having to reinitialize the Core, the total number of intervals to process is about 4.3 • 10 9 .

Table III shows the dependency between system frequency, number of Cores and task completion.

On our current FPGA prototyping system we conservatively estimate to be able to pack 200 PE which yields an realistic execution time of approximatively 50H.

D. Reality Check

We have tested the real performance of different configurations of our proposed systems on two FPGAs, our prototyping StratixIV development-kit and a large StratixV FPGA. The purpose of these tests was to show the performance of our solutions and to determine the degree of scalability of the proposed architectures (both the simpledispatcher and the credit-based dispatcher solution). The results are shown in Table IV First, the results obtained validate that, once the number of Cores exceeds a certain threshold (32 for StratixIV), the credit-based dispatch offers a more attractive solution. It is expected that, as the number of Cores scales up, the creditbased dispatch will continue to function at frequencies over 150 MHz, as the critical path of this system is in the adders of the TaMaDi Core.

Secondly, the results obtained on the StratixV FPGA prove that the increased capacity linearly improves the task completion time.

However, when reading this table one should consider that, as the size of the FPGA increases linearly, the time needed to compile the project on the FPGA increases at best polynomially. In other words, if for the 16 cores StratixIV design, compilation took some tens of minutes on a fast server, the StratixV designs took tens of hours to compile.

A solution to improve the compile/execution time ratio is to use a multi-FPGA based system, comprising of multiple similar FPGAs, such as the multi-FPGA prototyping board DN7020K10 from Dini Group [START_REF]Uncle of Monster' Altera Stratix IV ASIC Prototyping Engine[END_REF], comprising 16 Altera StratixIV FPGAs. The TaMaDi System would be compiled once, then replicated on these FPGAs. For simplicity one FPGA will also contain a dispatcher interface and will be connected to the host system. We estimate that one such system would complete the execution of one exponent in less than 2 hours.

All in all, depending on the available FPGA, the order of magnitude of the time required to process one input exponent is between a few hours and two days. Although the double precision/binary64 format has 2046 possible exponents, we do not need to perform such a calculation for every exponent: the exponential of a number larger than 710 is an overflow, and if |u| ≤ 2 -54 , then e u correctly rounded to that format is equal to 1. The most up-to-date implementation of the L-algorithm takes 45 hours to process one input exponent on a fairly recent FP core (AMD Opteron 2.19 GHz). Since these algorithm are very different and are run on very different machines, we suggest using both of them, which allows one to get much confidence in the obtained results.

CONCLUSION

We have suggested an algorithm and an FPGA architecture that make it possible to find hardest-to-round points for elementary functions in double precision. This requires huge computations, but they are done once for all, and allow one to design efficient libraries or hardware for elementary function evaluation. The achieved performance is slightly better than the one obtained using Lefèvre's L-algorithm but the real gain is not there: it lies in the fact that if, with a completely different method that runs on a completely different hardware, we obtain the same results, this gives much confidence in these results.

Figure 2 .

 2 Figure 2. Polynomial Evaluator based on the tabulated differences method

Figure 4 .

 4 Figure 4. Structure of one element in the ClusterInFIFO

Figure 5 .

 5 Figure 5. Global system dispatcher interface

Figure 6 .

 6 Figure 6. Credit based pipelined communication interface for one Cluster

Figure 7 .

 7 Figure 7. Global system architecture

Table I POST

 I PLACE-AND-ROUTE RESULTS OF THE TAMADI CORE PE

	Datapath Degree width N	FPGA	Frequency	Area LUTs REGs
			StratixIV	237 MHz	584	750
	120	4	StratixV Virtex5	359 MHz 262 MHz	585 640	750 646
			Virtex6	332 MHz	640	646

  circuit frequency K the number of TaMaDi Clusters M number of TaMaDi Cores within a Cluster N approximation polynomial degree w dp the datapath width on the TaMaDi Core

	HOST PC PCIExpress
	Reference Design	128 Avalon ST Interface Altera Hard IP EZDMA IP	FPGA
		DMA FIFOs
	Pipelined Credit Based
		Dispatcher Interface
	registers	...	...
	Interface	...	Interface
	TaMaDi		TaMaDi
	Cluster		Cluster

Table II DEPENDENCY

 II BETWEEN TAMADI CORE PARAMETERS, ITS AREA AND THE NECESSARY BANDWIDTH/CORE FOR A STRATIXIV. SIMILAR RESULTS HOLD FOR OTHER FPGAS

		Parameters		Area		Bandwidth
	N	w dp	kmax	w iid	LUTs REGs	(Mbit/s)
	4	120	2 20	32	584	750	0.017
	4	81	8,192	39	400	531	5.42