Nicolas Brisebarre

George Constantinides
email: g.constantinides@imperial.ac.uk

Miloš Ercegovac

Silviu-Ioan Filip

Matei Istoan
email: m.istoan@imperial.ac.uk

Jean-Michel Muller

A High Throughput Polynomial and Rational Function Approximations Evaluator

We present an automatic method for the evaluation of functions via polynomial or rational approximations and its hardware implementation, on FPGAs. These approximations are evaluated using Ercegovac's iterative E-method adapted for FPGA implementation. The polynomial and rational function coefficients are optimized such that they satisfy the constraints of the E-method. We present several examples of practical interest; in each case a resource-efficient approximation is proposed and comparisons are made with alternative approaches.

Introduction

We aim at designing a system able to approximate (in software) and then evaluate (in hardware) any regularenough function. More precisely, we try to minimize the sup norm of the difference between the function and the approximation in a given interval.

For particular functions, ad hoc solutions such as CORDIC [START_REF] Volder | The CORDIC computing technique[END_REF] or some specific tabulate-and-compute algorithms [START_REF] Wong | Fast hardware-based algorithms for elementary function computations using rectangular multipliers[END_REF] can be used. For low precision cases, table-based methods [START_REF] Sarma | Faithful bipartite ROM reciprocal tables[END_REF][START_REF] Schulte | Approximating elementary functions with symmetric bipartite tables[END_REF][START_REF] De Dinechin | Multipartite table methods[END_REF] methods are of interest. However, in the general case, piecewise approximations by polynomial or rational functions are the only reasonable solution. From a theoretical point of view, rational functions are very attractive, mainly because they can reproduce function behaviors (such as asymptotes, finite limits at ±∞) that polynomials do not satisfy. However, for software implementation, polynomials are frequently preferred to rational functions, because the latency of division is larger than the latency of multiplication. We aim at checking if rational approximations are of interest in hardware implementations. To help in the comparison of polynomial and rational approximations in hardware we use an algorithm, due to Ercegovac [START_REF] Ercegovac | A general method for evaluation of functions and computation in a digital computer[END_REF][START_REF]A general hardware-oriented method for evaluation of functions and computations in a digital computer[END_REF], called the E-method, that makes it possible to evaluate a degree-𝑛 polynomial, or a rational function of degree-𝑛 numerator and denominator at a similar cost without requiring division.

The E-method solves diagonally-dominant linear systems using a left-to-right digit-by-digit approach and has a simple and regular hardware implementation. It maps the problem of evaluating a polynomial or rational function into a linear system. The linear system corresponding to a given function does not necessarily satisfy the conditions of diagonal dominance. For polynomials, changes of variables allow one to satisfy the conditions. This is not the case for rational functions. There is however a family of rational functions, called E-fractions, that can be evaluated with the E-method in time proportional to the desired precision. One of our aims is, given a function, to decide whether it is better to approximate it by a polynomial or by an E-fraction. Furthermore, we want to design approximations whose coefficients satisfy some constraints (such as being exactly representable in a given format). We introduce algorithmic improvements with respect to [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF] for computing E-fractions. We 1 present a circuit generator for the E-method and compare its implementation on an FPGA with FloPoCo polynomial designs [START_REF] De Dinechin | On fixed-point hardware polynomials[END_REF] for several examples of practical interest. Since FloPoCo designs are pipelined (unrolled), we focus on an unrolled design of the E-method.

An Overview of the E-method

The E-method evaluates a polynomial 𝑃 𝜇 (𝑥) or a rational function 𝑅 𝜇,𝜈 (𝑥) by mapping it into a linear system. The system is solved using a left-to-right digit-by-digit approach, in a radix 𝑟 representation system, on a regular hardware. For a result of 𝑚 digits, in the range (-1, 1), the computation takes 𝑚 iterations. The first component of the solution vector corresponds to the value of 𝑃 𝜇 (𝑥) or 𝑅 𝜇,𝜈 (𝑥). Let 𝑅 𝜇,𝜈 (𝑥) = 𝑃 𝜇 (𝑥)

𝑄 𝜈 (𝑥) = 𝑝 𝜇 𝑥 𝜇 + 𝑝 𝜇-1 𝑥 𝜇-1 + • • • + 𝑝 0 𝑞 𝜈 𝑥 𝜈 + 𝑞 𝜈-1 𝑥 𝜈-1 + • • • + 𝑞 1 𝑥 + 1
where the 𝑝 𝑖 's and 𝑞 𝑖 's are real numbers. Let 𝑛 = max{𝜇, 𝜈}, 𝑝 𝑗 = 0 for 𝜇 + 1 𝑗 𝑛, and 𝑞 𝑗 = 0 for 𝜈 + 1 𝑗 𝑛. According to the E-method 𝑅 𝜇,𝜈 (𝑥) is mapped to a linear system 𝐿 : A × y = b:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 -𝑥 0 • • • 0 𝑞 1 1 -𝑥 0 • • • 0 𝑞 2 0 1 -𝑥 • • • 0 . . . . . . . . . . . . . . . 0 𝑞 𝑛-1 1 -𝑥 𝑞 𝑛 • • • 0 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑦 0 𝑦 1 𝑦 2 . . . . . . 𝑦 𝑛-1 𝑦 𝑛 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑝 0 𝑝 1 𝑝 2 . . . . . . 𝑝 𝑛-1 𝑝 𝑛 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (1)
so that 𝑦 0 = 𝑅 𝜇,𝜈 (𝑥). Likewise, 𝑦 0 = 𝑃 𝜇 (𝑥) when all 𝑞 𝑖 = 0.

The components of the solution vector y = [𝑦 0 , 𝑦 1 , . . . , 𝑦 𝑛 ] 𝑡 are computed, digit-by-digit, the most-significant digit first, by means of the following vector iteration:

w (𝑗) = 𝑟 × [ w (𝑗-1) -Ad (𝑗-1) ] , (2)
for 𝑗 = 1, . . . , 𝑚, where 𝑚 is the desired precision of the result. The term w (𝑗) is the vector residual in iteration 𝑗 with w (0) = [𝑝 0 , 𝑝 1 , . . . , 𝑝 𝑛 ] 𝑡 . The solution y is produced as a sequence of digit vectors:

d (𝑗-1) = [𝑑 (𝑗-1) 1 , . . . , 𝑑 (𝑗-1) 𝑛
] 𝑡 -a digit vector obtained in iteration 𝑗 -1 and used in iteration 𝑗. After 𝑚 iterations,

𝑦 𝑘 = ∑ 𝑚 𝑗=1 𝑑 (𝑗)
𝑘 𝑟 -𝑗 . The digits of the solution components 𝑦 0 , 𝑦 1 , . . . , 𝑦 𝑛 are computed using very simple scalar recurrences. Note that all multiplications in these recurrences use 𝑚 × 1 multipliers and that division required by the rational function is not explicitly performed.

𝑤

(𝑗) 𝑖 = 𝑟 × [ 𝑤 (𝑗-1) 𝑖 -𝑞 𝑖 𝑑 (𝑗-1) 0 -𝑑 (𝑗-1) 𝑖 + 𝑑 (𝑗-1) 𝑖+1 𝑥 ] , (3)
𝑤 (𝑗) 0 = 𝑟 × [ 𝑤 (𝑗-1) 0 -𝑑 (𝑗-1) 0 + 𝑑 (𝑗-1) 1 𝑥 ] (4)
and

𝑤 (𝑗) 𝑛 = 𝑟 × [ 𝑤 (𝑗-1) 𝑛 -𝑑 (𝑗-1) 𝑛 -𝑞 𝑛 𝑑 (𝑗-1) 0 ] . (5)
Initially, d (0) = 0. The radix-𝑟 digits 𝑑 (𝑗) 𝑖 are in the redundant signed digit-set 𝐷 𝜌 = {-𝜌, . . . , 0, 1, . . . 𝜌} with 𝑟/2 𝜌 𝑟 -1. If 𝜌 = 𝑟/2, 𝒟 𝜌 is called minimally redundant, and if 𝜌 = 𝑟 -1, it is maximally redundant. The choice of redundancy is determined by design considerations. The radix of computation is 𝑟 = 2 𝑘 so that internally radix-2 arithmetic is used. The residuals, in general, are in a redundant form to reduce the iteration time. Since the target is an FPGA technology which provides fast carry chains, we have non-redundant residuals. The digits 𝑑 to a single signed digit, following [START_REF]A general hardware-oriented method for evaluation of functions and computations in a digital computer[END_REF][START_REF] Ercegovac | Digital Arithmetic[END_REF]:

𝑑 (𝑗) 𝑖 = 𝑆(𝑤 (𝑗) 𝑖 ) = ⎧ ⎨ ⎩ sign (𝑤 (𝑗) 𝑖 ) × ⌊⃒ ⃒ ⃒𝑤 (𝑗) 𝑖 ⃒ ⃒ ⃒ + 1 2 ⌋ , if ⃒ ⃒ ⃒𝑤 (𝑗) 𝑖 ⃒ ⃒ ⃒ 𝜌, sign (𝑤 (𝑗) 𝑖 ) × ⌊⃒ ⃒ ⃒𝑤 (𝑗) 𝑖 ⃒ ⃒ ⃒ ⌋ , otherwise.
The selection is performed using a low-precision estimate ̂ 𝑤

(𝑗) 𝑖 of 𝑤 (𝑗)
𝑖 , obtained by truncating 𝑤 (𝑗) 𝑖 to one fractional bit.

Since the matrices considered here have 1s on the diagonal, a necessary condition for convergence is

∑ 𝑗̸ =𝑖 |𝑎 𝑖,𝑗 | < 1. Specifically, ⎧ ⎨ ⎩ ∀𝑖, |𝑝 𝑖 | 𝜉, ∀𝑖, |𝑥| + |𝑞 𝑖 | 𝛼, |𝑤 (𝑗) 𝑖 -̂ 𝑤 (𝑗) 𝑖 | ∆/2. (6)
where the bounds 𝜉, 𝛼, and ∆ satisfy [START_REF]A general hardware-oriented method for evaluation of functions and computations in a digital computer[END_REF]:

𝜉 = 1 2 (1 + ∆), 0 < ∆ < 1, 𝛼 (1 -∆)/(2𝑟) (7)
for maximally redundant digit sets used here. While the constraints (7) may seem restrictive, for polynomials, scaling techniques make it possible to satisfy them. However, this is not the case for all rational functions.

To remove this limitation the authors of [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF] have suggested the derivation of rational functions, called simple E-fractions, which are products of a power of 2 by a fraction that satisfies [START_REF]A general hardware-oriented method for evaluation of functions and computations in a digital computer[END_REF]. In this work we make further improvements to the rational functions based on E-fractions.

Outline of the paper

In Section 2, we discuss the effective generation of simple E-fractions, whose coefficients are exactly representable in a given format. Section 3 presents a hardware implementation of the E-method that targets FPGAs. In Section 4 we present and discuss some examples in various situations. We also present a comparison with FloPoCo implementations.

Effective computation of simple E-fractions

We show how to compute a simple E-fraction with fixed-point or floating-point coefficients. A first step (see Section 2.1), yields a simple E-fraction approximation with real coefficients to a function 𝑓 . In [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF], linear programming (LP) is used. Here, we use faster tools from approximation theory. This allows us to quickly check how far the approximation error of this E-fraction is from the optimal error of the minimax approximation (obtained using the Remez algorithm [START_REF] Cheney | Introduction to Approximation Theory[END_REF][START_REF] Powell | Approximation theory and methods[END_REF]), and how far it is from the error that an E-polynomial, with the same implementation cost, can yield. If this comparison suggests that it is more advantageous to work with an E-fraction, we use the Euclidean lattice basis reduction approach from [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF] for computing E-fractions with machine-number coefficients. We introduce in Section 2.2.2 a trick that improves its output.

Real approximation step

Let 𝑓 be a continuous function defined on [𝑎, 𝑏]. Let 𝜇, 𝜈 ∈ N be given and let R 𝜇,𝜈 (𝑥) = {𝑃/𝑄 : 𝑃 = ∑ 𝜇 𝑘=0 𝑝 𝑘 𝑥 𝑘 , 𝑄 = ∑ 𝜈 𝑘=0 𝑞 𝑘 𝑥 𝑘 , 𝑝 0 , . . . , 𝑝 𝜇 , 𝑞 0 , . . . , 𝑞 𝜈 ∈ R}. The aim is to compute a good rational fraction approximant 𝑅 ∈ R 𝜇,𝜈 (𝑥), with respect to the supremum norm defined by ‖𝑔‖ = sup 𝑥∈[𝑎, 𝑏] |𝑔(𝑥)|, to 𝑓 such that the real coefficients of 𝑅 (or 𝑅 divided by some fixed power of 2) satisfy the constraints imposed by the E-method.

As done in [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF], we can first apply the rational version of the Remez exchange algorithm [11, p. 173] to get 𝑅 ⋆ , the best possible rational approximant to 𝑓 among the elements of R 𝜇,𝜈 (𝑥). This algorithm can fail if 𝑅 ⋆ is degenerate or the choice of starting nodes is not good enough.

To bypass these issues, we develop the following process. It can be viewed as a Remez-like method of the first type, following ideas described in [11, p. 96-97] and [START_REF] Reemtsen | Modifications of the first Remez algorithm[END_REF]. It directly computes best real coefficient E-fractions with magnitude constraints on the denominator coefficients. If we remove these constraints, it will compute the minimax rational approximation, even when the Remez exchange algorithm fails.

We first show how to solve the problem over 𝑋, a finite discretization of [𝑎, 𝑏]. We apply a modified version (with denominator coefficient magnitude constraints) of the differential correction (DC) algorithm introduced in [START_REF] Cheney | Two new algorithms for rational approximation[END_REF]. It is given by Algorithm 1. System ( 8) is an LP problem and can be solved in practice very efficiently using a simplex-based LP solver. Convergence of this EDiffCorr procedure can be shown using an identical argument to the convergence proofs of the original DC algorithm [START_REF] Dua | Further remarks on the differential correction algorithm[END_REF][START_REF] Barrodale | The differential correction algorithm for rational ℓ∞-approximation[END_REF] (see Appendix A).

𝛿 ← max 𝑥∈𝑋 |𝑓 (𝑥) -𝑅(𝑥)| 4: find 𝑅 new = 𝑃 new /𝑄 new = ∑ 𝜇 𝑘=0 𝑝 ′ 𝑘 𝑥 𝑘 1 + ∑ 𝜈 𝑘=1 𝑞 ′ 𝑘 𝑥 𝑘 such that the expression max 𝑥∈𝑋 { |𝑓 (𝑥)𝑄 new (𝑥) -𝑃 new (𝑥)| -𝛿𝑄 new (𝑥) 𝑄(𝑥) } (8) subject to max 1 𝑘 𝜈 |𝑞 ′ 𝑘 | 𝑑, is minimized 5: 𝛿 new ← max 𝑥∈𝑋 |𝑓 (𝑥) -𝑅 new (𝑥)| 6: 𝑅 ← 𝑅 new 7: until |𝛿 -𝛿 new | < 𝜀 Algorithm 2 E-fraction Remez algorithm Input: 𝑓 ∈ 𝒞([𝑎, 𝑏]), 𝜇, 𝜈 ∈ N, finite set 𝑋 ⊆ [𝑎, 𝑏] with |𝑋| > 𝜇 + 𝜈, threshold 𝜀 > 0, coefficient magnitude bound 𝑑 > 0 Output: approximation 𝑅 ⋆ (𝑥) = ∑ 𝜇 𝑘=0 𝑝 ⋆ 𝑘 𝑥 𝑘 1 + ∑ 𝜈 𝑘=1 𝑞 ⋆ 𝑘 𝑥 𝑘 of 𝑓 over [𝑎, 𝑏] s.t. max 1 𝑘 𝜈 |𝑞 ⋆ 𝑘 | 𝑑 //
Compute best E-fraction approximation over 𝑋 using a // modified version of the differential correction algorithm 1: 𝑅 ⋆ ← EDiffCorr(𝑓, 𝜇, 𝜈, 𝑋, 𝜀, 𝑑)

2: 𝛿 ⋆ ← max 𝑥∈𝑋 |𝑓 (𝑥) -𝑅 ⋆ (𝑥)| 3: ∆ ⋆ ← max 𝑥∈[𝑎,𝑏] |𝑓 (𝑥) -𝑅 ⋆ (𝑥)| 4: while ∆ ⋆ -𝛿 ⋆ > 𝜀 do 5: 𝑥 new ← argmax 𝑥∈[𝑎,𝑏] |𝑓 (𝑥) -𝑅 ⋆ (𝑥)| 6: 𝑋 ← 𝑋 ∪ {𝑥 new } 7:
𝑅 ⋆ ← EDiffCorr(𝑓, 𝜇, 𝜈, 𝑋, 𝜀, 𝑑)

8: 𝛿 ⋆ ← max 𝑥∈𝑋 |𝑓 (𝑥) -𝑅 ⋆ (𝑥)| 9: ∆ ⋆ ← max 𝑥∈[𝑎,𝑏] |𝑓 (𝑥) -𝑅 ⋆ (𝑥)| 10: end while
To address the problem over [𝑎, 𝑏], Algorithm 2 solves a series of best E-fraction approximation problems on a discrete subset 𝑋 of [𝑎, 𝑏], where 𝑋 increases at each iteration by adding a point where the current residual term achieves its global maximum.

Our current experiments suggest that the speed of convergence for Algorithm 2 is linear. We can potentially decrease the number of iterations by adding to 𝑋 more local extrema of the residual term at each iteration. Other than its speed compared to the LP approach from [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF], Algorithm 2 will generally converge to the best E-fraction approximation with real coefficients over [𝑎, 𝑏], and not on a discretization of [𝑎, 𝑏].

Once 𝑅 ⋆ is computed, we determine the least integer 𝑠 such that the coefficients of the numerator of 𝑅 ⋆ divided by 2 𝑠 fulfill the first condition of [START_REF] Ercegovac | A general method for evaluation of functions and computation in a digital computer[END_REF]. It gives us a decomposition 𝑅 ⋆ (𝑥) = 2 𝑠 𝑅 𝑠 (𝑥). 𝑅 𝑠 is thus a rescaled version of 𝑅. We take 𝑓 𝑠 = 2 -𝑠 𝑓 to be the corresponding rescaling of 𝑓 . The magnitude bound 𝑑 is usually equal to 𝛼 -max(|𝑎|, |𝑏|), allowing the denominator coefficients to be valid with respect to the second constraint of [START_REF] Ercegovac | A general method for evaluation of functions and computation in a digital computer[END_REF].

Both Algorithm 1 and 2 can be modified to compute weighted error approximations, that is, work with a norm of the form ‖𝑔‖ = max 𝑥∈[𝑎,𝑏] |𝑤(𝑥)𝑔(𝑥)|, where 𝑤 is a continuous and positive weight function over [𝑎, 𝑏]. This is useful, for instance, when targeting relative error approximations. The changes are minimal and consist only of introducing the weight factor in the error computations in lines 3, 5 of Algorithm 1, lines 2, 3, 5, 8, 9 of Algorithm 2 and changing [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF] with

max 𝑥∈𝑋 { 𝑤(𝑥) |𝑓 (𝑥)𝑄 new (𝑥) -𝑃 new (𝑥)| -𝛿𝑄 new (𝑥) 𝑄(𝑥) } .
The weighted version of the DC algorithm is discussed, for instance, in [START_REF] Dudgeon | Recursive filter design using differential correction[END_REF].

Lattice basis reduction step

Our goal is to compute a simple E-fraction

̂ 𝑅(𝑥) = ∑ 𝜇 𝑗=0 ̂ 𝑝 𝑗 𝑥 𝑗 1 + ∑ 𝜈 𝑗=1 ̂ 𝑞 𝑗 𝑥 𝑗 ,
where ̂ 𝑝 𝑗 and ̂ 𝑞 𝑗 are fixed-point or floating-point numbers [START_REF] Ercegovac | Digital Arithmetic[END_REF][START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF], that is as close as possible to 𝑓 𝑠 , the function we want to evaluate. These unknown coefficients are of the form 𝑀 2 𝑒 , 𝑀 ∈ Z:

• for fixed-point numbers, 𝑒 is implicit (decided at design time);

• for floating-point numbers, 𝑒 is explicit (i.e., stored). A floating-point number is of precision

𝑡 if 2 𝑡-1 𝑀 < 2 𝑡 -1.
A different format can be used for each coefficient of the desired fraction. If we assume a target format is given for each coefficient, then a straightforward approach is to round each coefficient of 𝑅 𝑠 to the desired format. This yields what we call in the sequel a naive rounding approximation. Unfortunately, this can lead to a significant loss of accuracy. We first briefly recall the approach from [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF] that makes it possible to overcome this issue. Then, we present a small trick that improves on the quality of the output of the latter approach. Eventually, we explain how to handle a coefficient saturation issue appearing in some high radix cases.

Modeling with a closest vector problem in a lattice

Every fixed-point number constraint leads to a corresponding unknown 𝑀 , whereas each precision-𝑡 floating-point number leads to two unknowns 𝑀 and 𝑒. A heuristic trick is given in [START_REF] Brisebarre | Efficient polynomial 𝐿 ∞ approximations[END_REF] to find an appropriate value for each 𝑒 in the floating-point case: we assume that the coefficient in question from ̂ 𝑅 will have the same order of magnitude as the corresponding one from 𝑅 𝑠 , hence they have the same exponent 𝑒. Once 𝑒 is set, the problem is reduced to a fixed-point one.

Then, given 𝑢 0 , . . . , 𝑢 𝜇 , 𝑣 1 , . . . , 𝑣 𝜈 ∈ Z, we have to determine 𝜇 + 𝜈 + 1 unknown integers 𝑎 𝑗 (= ̂ 𝑝 𝑗 2 -𝑢𝑗 ) and 𝑏 𝑗 (= ̂ 𝑞 𝑗 2 -𝑣𝑗 ) such that the fraction

̂ 𝑅(𝑥) = ∑ 𝜇 𝑗=0 𝑎 𝑗 2 𝑢𝑗 𝑥 𝑗 1 + ∑ 𝜈 𝑗=1 𝑏 𝑗 2 𝑣𝑗 𝑥 𝑗
is a good approximation to 𝑅 𝑠 (and 𝑓 𝑠 ), i.e., ‖ ̂ 𝑅 -𝑅 𝑠 ‖ is small. To this end, we discretize the latter condition in 𝜇 + 𝜈 + 1 points 𝑥 0 < • • • < 𝑥 𝜇+𝜈 in the interval [𝑎, 𝑏], which gives rise to the following instance of a closest vector problem, one of the fundamental questions in the algorithmics of Euclidean lattices [START_REF] Nguyen | The LLL Algorithm -Survey and Applications, ser. Information Security and Cryptography[END_REF]: we want to compute 𝑎 0 , . . . , 𝑎 𝜇 , 𝑏 1 , . . . , 𝑏 𝜈 ∈ Z such that the vectors

𝜇 ∑ 𝑗=0 𝑎 𝑗 𝛼 𝑗 - 𝑛 ∑ 𝑗=1 𝑏 𝑗 𝛽 𝑗 and r (9)
are as close as possible, where

𝛼 𝑗 = [2 𝑢𝑗 𝑥 𝑗 0 , . . . , 2 𝑢𝑗 𝑥 𝑗 𝜇+𝜈 ] 𝑡 , 𝛽 𝑗 = [2 𝑣𝑗 𝑥 𝑗 0 𝑅 𝑠 (𝑥 0 ), . . . , 2 𝑣𝑗 𝑥 𝑗 𝜇+𝜈 𝑅 𝑠 (𝑥 𝜇+𝜈 )]
𝑡 and r = [𝑅 𝑠 (𝑥 0 ), . . . , 𝑅 𝑠 (𝑥 𝜇+𝜈 )] 𝑡 . It can be solved in an approximate way very efficiently by applying techniques introduced in [START_REF] Lenstra | Factoring polynomials with rational coefficients[END_REF] and [START_REF] Babai | On Lovász' lattice reduction and the nearest lattice point problem[END_REF]. We refer the reader to [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF][START_REF] Brisebarre | Efficient polynomial 𝐿 ∞ approximations[END_REF] for more details on this and how the discretization 𝑥 0 , . . . , 𝑥 𝜇+𝜈 should be chosen. We thus propose to fix the problematic values of 𝑏 𝑗 to the closest value to the allowable limit that does not break the E-method magnitude constraints. The change is minor in (9); we just move the corresponding vectors in the second sum on the left hand side of (9) to the right hand side with opposite sign. The resulting problem can also be solved using the tools from [START_REF] Brisebarre | An Efficient Method for Evaluating Polynomial and Rational Function Approximations[END_REF][START_REF] Brisebarre | Efficient polynomial 𝐿 ∞ approximations[END_REF]. This usually gives a valid simple E-fraction ̂ 𝑅 of very good quality.

Higher radix problems

Coefficient saturation issues get more pronounced by increasing the radix 𝑟. In such cases, care must also be taken with the approximation domain: the

|𝑞 𝑗 | upper magnitude bound 𝛼 -max(|𝑎|, |𝑏|) can become negative, since 𝛼 = (1 -∆)/(2𝑟) → 0 as 𝑟 → ∞.
To counter this, we use argument and domain scaling ideas presented in [START_REF] Brisebarre | Functions approximable by E-fractions[END_REF]. This basically consists in approximating 𝑓 (𝑥) = 𝑓 (2 𝑡 𝑦), for 𝑦 ∈ [2 -𝑡 𝑎, 2 -𝑡 𝑏] as a function in 𝑦. If 𝑡 > 0 is large enough, then the new |𝑞 𝑗 | bound 𝛼 -max(|2 -𝑡 𝑎|, |2 -𝑡 𝑏|) will be 0.

A hardware implementation targeting FPGAs

We now focus on the hardware implementation of the E-method on FPGAs. This section introduces a generator capable of producing circuits that can solve the system A • y = b, through the recurrences of Equations ( 3)-( 5).

The popularity of FPGAs is due to their ability to be reconfigured, and their relevance in prototyping as well as in scientific and high-performance computing. They are composed of large numbers of small look-up tables (LUTs), with 4-6 inputs and 1-2 outputs, and each can store the result of any logic function of their inputs. Any two LUTs on the device can communicate, as they are connected through a programmable interconnect network. Results of computations can be stored in registers, usually two of them being connected to the output of each LUT. These features make of FPGAs a good candidate as a platform for implementing the E-method, as motivated even further below.

A minimal interface

An overview of the hardware back-end is presented in Figure 1. Although not typically accessible to the user, its interface is split between the functional and the performance specification. The former consists of the radix 𝑟, the input and output formats, specified as the weights of their most significant (MSB) and least significant (LSB) bits, the parameter ∆ (input by the user) and the coefficients of the polynomials 𝑃 𝜇 (𝑥) and 𝑄 𝜈 (𝑥) (coming from the rational approximation). Having 𝑚𝑠𝑏 𝑖𝑛 as a parameter is justified by some of the examples of Section 4, where even though the input 𝑥 belongs to [-1, 1], the maximum value it is allowed to have is smaller, given by the constraints ( 6) and [START_REF]A general hardware-oriented method for evaluation of functions and computations in a digital computer[END_REF]. This allows for optimizing the datapath sizes. The 𝑚𝑠𝑏 𝑜𝑢𝑡 can be computed during the polynomial/rational approximation phase and passed to the generator.

The circuit generator is developed inside the FloPoCo framework [START_REF] De Dinechin | Designing custom arithmetic data paths with FloPoCo[END_REF], which facilitates the support of classical parameters in the performance specification, such as the target frequency of the circuit, or the target FPGA device. It also means that we can leverage on the automatic pipelining and test infrastructure present in the framework, alongside the numerous existing arithmetic operators.

w (j-1) i d (j-1) 0 d (j-1) i 𝜌x . . . x 0 -x . . . (-𝜌)x d (j-1) i+1 KCM ×𝑞 𝑖 MUX Bitheap ≪ w (j) i + - - + CU i

Implementation details

An overview of the basic iteration, based on Equation (3), is presented in Figure 2. As this implementation is targeted towards FPGAs, several optimizations can be applied. First, the multiplication 𝑑 (𝑗-1) 0

• 𝑞 𝑖 can be computed using the KCM technique for multiplying by a constant [START_REF] Chapman | Fast integer multipliers fit in FPGAs (EDN 1993 design idea winner)[END_REF][START_REF] Wirthlin | Constant coefficient multiplication using look-up tables[END_REF], with the optimizations of [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF], that extend the method for real constants. Therefore, instead of using dedicated multiplier blocks (or of generating partial products using LUTs), we can directly tabulate the result of the multiplication 𝑑 (𝑗-1) 0

• 𝑞 𝑖 , at the cost of one LUT per output bit of the result. This remains true even for higher radices, as LUTs on modern devices can accommodate functions of 6 Boolean inputs.

A second optimization relates to the term 𝑑 𝑖+1 . The multiplications by the negative values in the digit set come at the cost of just one bitwise operation and an addition, which are combined in the same LUT by the synthesis tools.

Finally, regarding the implementation of the CUs, the multi-operand addition of the terms of Equation ( 3) is implemented using a bitheap [START_REF] Brunie | Arithmetic core generation using bit heaps[END_REF]. The alignments of the accumulated terms and their varied sizes would make for a wasteful use of adders. Using a bitheap we have a single, global optimization of the accumulation. In addition, managing sign extensions comes at the cost of a single additional term in the accumulation, using a technique from [START_REF] Ercegovac | Digital Arithmetic[END_REF]. The sign extension of a two's complement fixed point number sxx . . . xx is performed as: 00...0sxxxxxxx + 11...110000000 = ss...ssxxxxxxx

The sum of the constants is computed in advance and added to the accumulation. The final shift comes at just the cost of some routing, since the shift is by a constant amount.

Modern FPGAs contain fast carry chains. Therefore, we represent the components of the residual vector 𝑤 using two's complement, as opposed to a redundant representation. The selection function uses an estimate of 𝑤 𝑖 with one fractional bit allowing a simple tabulation using LUTs.

Iteration 0, the initialization, comes at almost no cost, and can be done through fixed routing in the FPGA. This is also true for the second iteration, as simplifying the equations results in 𝑤 can be pre-computed and stored directly. This not only saves one iteration, but also improves the accuracy, as 𝑤 (1) 𝑖 and 𝑑

(1) 𝑖 can be pre-computed using higher-precision calculations. Going one iteration further, we can see that most of the computations required for 𝑤 (2) 𝑖 can also be done in advance, except, of course, those involving 𝑥.

Figure 3 shows an unrolled implementation of the E-method, that uses the CUs of Figure 2 as basic building blocks. At the top of the architecture, a scaling can be applied to the input. This step is optional, and the scale factor (optional parameter in the design) can either be set by the user, or computed by the generator so that, given the input format, the parameter ∆ and the coefficients of 𝑃 and 𝑄, the scaled input satisfies the constraints ( 6) and [START_REF]A general hardware-oriented method for evaluation of functions and computations in a digital computer[END_REF]. The multiplications between 𝑥 and the possible values of the digits 𝑑 (𝑗) 𝑖 are done using the classical shift-and-add technique, a choice justified by the small values of the constants and the small number of bits equal to 1 in their representations. At the bottom of Figure 3, the final result 𝑦 0 is obtained in two's complement representation. Again, this step is also optional, as users might be content with having the result in the redundant representation.

There is one more optimization that can be done here due to an unrolled implementation. Because only the 𝑑 (𝑗) 0 digits are required to compute 𝑦 0 , after iteration 𝑚 -𝑛 we can compute one less element of w (𝑗) and d (𝑗) at each iteration. This optimization is the most effective when the number of required iterations 𝑚 is comparable to 𝑛, in which case the required hardware is reduced to almost half.

Error Analysis

To obtain a minimal implementation for the circuit described in Figure 3, we need to size the datapaths in a manner that guarantees that the output 𝑦 0 remains unchanged, with respect to an ideal implementation, free of potential rounding errors. To that end, we give an error-analysis, which follows [START_REF] Ercegovac | A general method for evaluation of functions and computation in a digital computer[END_REF]Ch. 2.8]. For the sake of brevity, we focus on the radix 2 case.

In order for the circuit to produce correct results, we must ensure that the rounding errors do not influence the selection function: 𝑆( ̃ 𝑤

(𝑗) 𝑖 ) = 𝑆(𝑤 (𝑗) 𝑖 ) = 𝑑 (𝑗)
𝑖 , where the tilded terms represent approximate values. In [START_REF] Ercegovac | A general method for evaluation of functions and computation in a digital computer[END_REF], the idea is to model the rounding errors due to the limited precision used to store the coefficients 𝑝 𝑗 and 𝑞 𝑗 inside the matrix 𝐴 as a new error matrix E A = (𝜀 𝑖𝑗 ) 𝑛×𝑛 . With the method introduced in this paper, the coefficients are machine representable numbers, and therefore incur no additional error. What remains to deal with are errors due to the limited precision of the involved operators. The only one that could produce rounding errors is the multiplication 𝑑 (𝑗-1) 0

• 𝑞 𝑖 . We know that 𝑑 (𝑗-1) 0 1 (the case 𝑑 (𝑗-1) 0 = 0 is clearly not a problem), so the LSB of 𝑑 (𝑗-1) 0

• 𝑞 𝑖 is at least that of 𝑞 𝑖 , if not larger. If the output precision satisfies 𝑙𝑠𝑏 𝑜𝑢𝑡 𝑙𝑠𝑏 𝑞𝑖 (which is usually the case), we perform this operation on its full precision, so we do not require any additional guard bits for the internal computations. If this assumption does not hold, based on [START_REF] Ercegovac | A general method for evaluation of functions and computation in a digital computer[END_REF], we obtain the following expression for the rounding errors introduced when computing w (𝑗) inside Equation ( 2), denoted with 𝜀 (𝑗)

w :

𝜀 (𝑗) w = 2 • (𝜀 (𝑗-1) w + 𝜀 𝑐𝑜𝑛𝑠𝑡_𝑚𝑢𝑙𝑡 + E A • d (𝑗-1) ).
We can thus obtain an expression for 𝜀 (𝑚) w , the error vector at step 𝑚, where 𝑚 is the bitwidth of 𝑦 0 and 𝜀 𝑐𝑜𝑛𝑠𝑡_𝑚𝑢𝑙𝑡 are the errors due to the constant multipliers. Since 𝜀 (0)

w = 0, 𝜀 (𝑚) 𝑤0 = 2 𝑚 • (𝜀 𝑐𝑜𝑛𝑠𝑡_𝑚𝑢𝑙𝑡 + ‖E A ‖ • 𝑚 ∑ 𝑗=1 𝑑 (𝑗) 0 • 2 -𝑗 ),
where ‖E A ‖ is the matrix 2-norm. We use a larger intermediary precision for the computations, with 𝑔 extra guard bits. Therefore, we can design a constant multiplier for which 𝜀 𝑐𝑜𝑛𝑠𝑡_𝑚𝑢𝑙𝑡 ε𝑐𝑜𝑛𝑠𝑡_𝑚𝑢𝑙𝑡 2 -𝑚-𝑔 . Also,

‖E A ‖ max 𝑖 𝑛 ∑ 𝑗=1 |𝜀 𝑖𝑗 | and 𝑚 ∑ 𝑗=𝑖 𝑑 (𝑗) 0 • 2 -𝑗 < 1,
hence we can deduce that for each 𝑤 (𝑚) 𝑖

we have 𝜀

(𝑚) 𝑤𝑖 ε(𝑚) 𝑤𝑖 2 𝑚 (2 -𝑚-𝑔 + 𝑛 • 2 -𝑚-𝑔 ).
In order for the method to produce correct results, we need to ensure that ε(𝑚) 𝑤 ∆/2, therefore we need to use 𝑔 2 + log 2 (2(𝑛 + 1)/∆) additional guard bits. This also takes into consideration the final rounding to the output format.

Examples, Implementation and Discussion

In this section, we consider fractions with fixed-point coefficients of 24, 32, and 48 bits: these coefficients will be of the form 𝑖/2 𝑚 , with -2 𝑚 𝑖 2 𝑚 , where 𝑤 = 24, 32, 48. The target approximation error in each case is 2 -𝑚 , i.e., ∼ 5.96 • 10 -8 , 2.33 • 10 -10 , 3.55 • 10 -15 respectively.

Examples. All the examples are defined in the first column of Table 1. When choosing them we considered:

• Functions useful in practical applications. The exponential function (Example 2) is a ubiquitous one. Functions of the form log 2 (1 + 2 ±𝑘𝑥 ) (as the one of Example 3) are useful when implementing logarithmic number systems.The erf function (Example 4) is useful in probability and statistics, while the Bessel function 𝐽 0 (Example 5) has many applications in physics.

• Functions that illustrate the various cases that can occur: polynomials are a better choice (Example 3); rational approximation is better (Examples 1, 2, Example 4 if 𝑟 8 and Example 5 if 𝑟 = 2). We also include instances where the approximating E-fractions are very different from the minimax, unconstrained, rational approximations with similar degrees in the numerator and denominator (Examples 1 and 2).

All the examples start with a radix 2 implementation after which higher values of 𝑟 are considered. Table 1 displays approximation errors in the real coefficient and fixed-point coefficient E-fraction cases. Notice in particular the lattice-based approximation errors, which are generally much better than the naive rounding ones. We also give some complementary comments. Example 1. The type (4, 4) rational minimax unconstrained approximation error is 4.59 • 10 -16 , around 5 orders of magnitude smaller than the E-fraction error. A similar difference happens in case of Example 2, where the type [START_REF] Sarma | Faithful bipartite ROM reciprocal tables[END_REF][START_REF] Sarma | Faithful bipartite ROM reciprocal tables[END_REF] unconstrained minimax approximation has error 2.26 • 10 -16 .

Example 2. In this case, we are actually working with a rescaled input and are equivalently approximating exp(2𝑥), 𝑥 ∈ [0, 7/128]. Also, for 𝑟 = 8, the real coefficient E-fraction is the same as the E-polynomial one (the magnitude constraint for the denominator coefficients is 0).

Example 3. Starting with 𝑟 = 8, we have to scale both the argument 𝑥 and the approximation domain by suitable powers of 2 for the E-method constraints to continue to hold (see end of Section 2.1).

Example 4. As with the previous example, for 𝑟 = 16, 32 we have to rescale the argument and interval to get a valid E-polynomial.

Example 5. By a change of variable, we are actually working with 𝐽 0 (2𝑥 -1/16), 𝑥 ∈ [0, 1/16]. If we consider 𝑟 16, the 48 bits used to represent the coefficients were not sufficient to produce an approximation with error below 2 -48 .

Implementation. We have generated the corresponding circuits for each of the examples, and synthesized them. The target platform is a Xilinx Virtex6 device xc6vcx75t-2-ff484, and the toolchain used is ISE 14.7. The resulting circuit descriptions are in an easily readable and portable VHDL. For each of the examples compare against a state of the art polynomial approximation implementation generated by FloPoCo (described in [START_REF] De Dinechin | On fixed-point hardware polynomials[END_REF]). FloPoCo [START_REF] De Dinechin | Designing custom arithmetic data paths with FloPoCo[END_REF] is an open-source arithmetic core generator and library for FPGAs. It is, to the best of our knowledge, one of the only alternatives capable of producing the functions chose for comparison. Table 2 presents the results.

At the top of Table 2, for Example 1, we show the flexibility of the generator: it can achieve a wide range of latencies and target frequencies. The examples show how the frequency can be scaled from around 100MHz to 300MHz, at the expense of a deeper pipeline and an increased number of registers. The number of registers approximately doubles each time the circuit's period is reduced by a factor 2. This very predictable behavior should help the end user make an acceptable trade-off in terms of performance to required resources. The frequency cap of 300MHz is not something inherent to the E-method algorithm, neither to the implementation. Instead it comes from pipeline performance issues of the bitheap framework inside the FloPoCo generator. We expect that once this bottleneck is fixed, our implementations will reach much higher target frequencies, without the need for any modifications in the current implementation.

Discussion. Examples 1 and 2 illustrate that for functions where classical polynomial approximation techniques, like the one used in FloPoCo, manage to find solutions of a reasonably small degree, the ensuing architectures also manage to be highly efficient. This shows, as implementations produced by FloPoCo (with polynomials of degree 6 in both cases) are twice (if not more) as efficient in terms of resources.

However, this is no longer the case when E-fractions can provide a better approximation. This is reflected by Examples 3 to 5, where we obtain a more efficient solution, by quite a large margin in some cases.

For Example 5, Table 2 does not present any data for the FloPoCo implementation as they do not currently support this type of function.

There are a few remarks to be made regarding the use of a higher radix in the implementations of the E-method. Example 4 is an indication that the overall delay of the architecture reaches a point where it can no longer benefit from increasing the radix. The lines of Table 2 marked with an asterisk were generated with an alternative implementation for the CUs, which uses multipliers for computing the 𝑑 (𝑗-1) 𝑖+1 • 𝑥 products. This is due to the exponential increase of the size of multiplexers with the increase of the radix, while the equivalent multiplier only increases linearly. Therefore, there is a crossover point from which it is best to use this version of the architecture, usually at radix 8 or 16. Finally, the effects of truncating the last iterations become the most obvious when the maximum degree 𝑛 is close to the number of required iterations 𝑚 in radix 𝑟. This effect can be observed for Example 3 and 4, where there is a considerable drop in resource consumption between the use of radix 8 and 16, and 16 and 32, respectively.

Summary and Conclusions

A high throughput system for the evaluation of functions by polynomials or rational functions using simple and uniform hardware is presented. The evaluation is performed using the unfolded version of the E-method, with a latency proportional to the precision. An effective computation of the coefficients of the approximations is given and the best strategies (choice of polynomial vs rational approximation, radix of the iterations) investigated. Designs using a circuit generator for the E-method inside the FloPoCo framework are developed and implemented using FPGAs for five different functions of practical interest, using various radices. As it stands, and as our examples show for FPGA devices, the E-method is generally more efficient as soon as the rational approximation is significantly more efficient than the polynomial one. From a hardware standpoint, the results show it is desirable to use the E-method with high radices, usually at least 8. The method also becomes efficient when we manage to find a balance between the maximum degree 𝑛 of the polynomial or E-fraction and the number of iterations required for converging to a correct result, which we can control by varying the radix. An open-source implementation of our approach will soon be available online at https://github.com/sfilip/emethod.

A Convergence of the E-fraction differential correction algorithm

The E-fraction approximation problem can be formulated in terms of finding functions of the form 𝑅(𝑥) = 𝑃 (𝑥)/𝑄(𝑥) :=

∑ 𝜇 𝑗=0 𝛼𝑗𝑥 𝑗 1 + ∑ 𝜈 𝑗=1 𝛽𝑗𝑥 𝑗 , (10)
that best approximate 𝑓 , where the denominator coefficients are bounded in magnitude, i.e., max

1 𝑘 𝜈 |𝛽 𝑘 | 𝑑 > 0, (11)
and 𝑄(𝑥) > 0, ∀𝑥 ∈ 𝑋. We denote this set as ℛ𝜇,𝜈 (𝑋). In this case, the differential correction algorithm is defined as follows: given 𝑅 𝑘 = 𝑃 𝑘 /𝑄 𝑘 ∈ ℛ𝜇,𝜈 (𝑋) satisfying ( 10) and [START_REF] Cheney | Introduction to Approximation Theory[END_REF], find 𝑅 = 𝑃/𝑄 that minimizes the expression 𝑀 over 𝑋 for any rational function satisfying [START_REF] Ercegovac | Digital Arithmetic[END_REF] and [START_REF] Cheney | Introduction to Approximation Theory[END_REF]

where 𝑧 𝑘 is an element of 𝑋 minimizing 𝑄 𝑘 (𝑥) 𝑄 𝑘+1 (𝑥)

. Since 𝛿 𝑘 converges to 𝛿 > 𝛿, [START_REF] Reemtsen | Modifications of the first Remez algorithm[END_REF] It follows that ∏ 𝑥∈𝑋 𝑄 𝑘 (𝑥) diverges as 𝑘 → ∞. This contradicts the normalization condition [START_REF] Cheney | Introduction to Approximation Theory[END_REF].

  (𝑗) 𝑖 are selected so that the residuals |𝑤 (𝑗) 𝑖 | remain bounded. The digit selection is performed by rounding the residuals 𝑤 (𝑗) 𝑖

Figure 2 :

 2 Figure 2: The basic Computation Unit (CU).

(𝑗- 1 )

 1 𝑖+1 • 𝑥, from Equation (3). Since 𝑑 (𝑗-1) 𝑖+1 ∈ {-𝜌, . . . , 𝜌}, we can compute the products 𝑥 • 𝜌, 𝑥 • (𝜌 -1), 𝑥 • (𝜌 -2), . . . , only once, and then select the relevant one based on the value of 𝑑 (𝑗-1)

Figure 3 :

 3 Figure 3: The E-method circuit generator.

  max 𝑥∈𝑋 { |𝑓 (𝑥)𝑄(𝑥) -𝑃 (𝑥)| -𝛿 𝑘 𝑄(𝑥) 𝑄 𝑘 (𝑥) } , subject to max 1 𝑘 𝜈 |𝛽𝑗| 𝑑, where 𝛿 𝑘 = max𝑥∈𝑋 |𝑓 (𝑥) -𝑅 𝑘 (𝑥)|. If 𝑅 = 𝑃/𝑄 is not good enough, continue with 𝑅 𝑘+1 = 𝑅. The convergence properties of this variation are similar to those of the original differential correction algorithm. Let 𝛿 * = inf 𝑅∈ℛ𝜇,𝜈 (𝑋) ‖𝑓 -𝑅‖, where ‖𝑔‖ = max𝑥∈𝑋 |𝑔(𝑥)|. Theorem 1. If 𝑄 𝑘 (𝑥) > 0, ∀𝑥 ∈ 𝑋 and if 𝛿 𝑘 ̸ = 𝛿, then 𝛿 𝑘+1 < 𝛿 𝑘 and 𝑄 𝑘+1 (𝑥) > 0, ∀𝑥 ∈ 𝑋. Proof. Since 𝛿 𝑘 ̸ = 𝛿 * , there exists 𝑅 = 𝑃 /𝑄 ∈ ℛ𝜇,𝜈 (𝑋) s.t. ‖𝑓 -𝑅‖ = 𝛿 < 𝛿 𝑘 . Let Δ = min𝑥∈𝑋 ⃒ ⃒ 𝑄(𝑥) ⃒ ⃒ and 𝑀 = 1 + 𝑑 max𝑥∈𝑋 ∑ 𝜈 𝑘=1 ⃒ ⃒ 𝑥 𝑘 ⃒ ⃒ so that |𝑄(𝑥)|

  Algorithm 1 E-fraction EDiffCorr algorithm Input: 𝑓 ∈ 𝒞([𝑎, 𝑏]), 𝜇, 𝜈 ∈ N, finite set 𝑋 ⊆ [𝑎, 𝑏] with |𝑋| > 𝜇 + 𝜈, threshold 𝜀 > 0, coefficient magnitude bound 𝑑 > 0 𝑘 𝑥 𝑘 of 𝑓 over 𝑋 s.t. max 1 𝑘 𝜈 |𝑞 𝑘 | 𝑑 // Initialize the iterative procedure (𝑅 = 𝑃/𝑄)

	Output: approximation 𝑅(𝑥) =	∑ 𝜇 𝑘=0 𝑝 𝑘 𝑥 𝑘 1 + ∑ 𝜈 𝑘=1 𝑞

1: 𝑅 ← 1 2: repeat 3:

  While we generally obtain integer 𝑎 𝑗 and 𝑏 𝑗 which correspond to a good approximation, the solution is not always guaranteed to give a valid simple E-fraction. What happens is that, in many cases, some of the denominator coefficients in 𝑅 𝑠 are maximal with respect to the magnitude constraint in (6) (recall that the second line in (6) can be restated as |𝑞 𝑗 | 𝛼 -max(|𝑎|, |𝑏|)). In this context, the corresponding values of |𝑏 𝑗 | are usually too large.

	radix 𝑟
	input format (𝑚𝑠𝑏 𝑖𝑛 , 𝑙𝑠𝑏 𝑖𝑛 )
	output format (𝑚𝑠𝑏 𝑜𝑢𝑡 , 𝑙𝑠𝑏 𝑜𝑢𝑡 )	Circuit generator	.vhdl
	0 < ∆ < 1
	(𝑝 𝑗 ) 0 𝑗<𝜇 and (𝑞 𝑗 ) 0 𝑗<𝜈	FPGA	frequency
	Functional specification	Performance specification
	Figure 1: Circuit generator overview.
	2.2.2 A solution to a coefficient saturation issue

Table 1 :

 1 Approximation errors in the real coefficient and fixed-point coefficient E-fraction cases

			Function Type of error	Δ	𝑟	(𝜇, 𝜈) 𝑚	Real coefficient E-fraction error	Lattice-based error	Naive rounding error
	Ex. 1	√	1 + (9𝑥/2) 4 , 𝑥 ∈ [0, 1/32] absolute	1 8	2 4 8	(4, 4) 32	5.22 • 10 -11 6.32 • 10 -11 8.25 • 10 -11	5.71 • 10 -11 4.93 • 10 -10 1.11 • 10 -10	1.11 • 10 -9 7 • 10 -11 1.78 • 10 -9
	Ex. 2		exp(𝑥), 𝑥 ∈ [0, 7/64] relative	1 8	2 4 8	(3, 3) (4, 4) (5, 0)	32	1.64 • 10 -10 10 -12 1.16 • 10 -12	1.94 • 10 -10 1.11 • 10 -12 1.39 • 10 -12	3.24 • 10 -10 1.91 • 10 -11 1.74 • 10 -11
	Ex. 3	log 2 (1 + 2 -16𝑥 ), 𝑥 ∈ [0, 1/16] absolute	1 2	2 4,8,16	(5, 5) 24 (5, 0)	1.98 • 10 -8 2.04 • 10 -8	2.33 • 10 -8 2.64 • 10 -8	4.37 • 10 -7 4.22 • 10 -7
	Ex. 4		erf(𝑥), 𝑥 ∈ [0, 1/32] absolute	1 8	2 4 8,16,32 (5, 0) (4, 4) (4, 4)	48	2.92 • 10 -17 3.44 • 10 -17 1.34 • 10 -15	3.43 • 10 -17 4.23 • 10 -17 1.64 • 10 -15	1.67 • 10 -16 1.13 • 10 -16 2.7 • 10 -15
	Ex. 5	𝐽0(𝑥), 𝑥 ∈ [-1/16, 1/16] relative	1 2	2 4,8	(4, 4) 48 (6, 0)	2.15 • 10 -17 1.23 • 10 -17	2.37 • 10 -15 2.37 • 10 -15	2.49 • 10 -15 2.53 • 10 -15

Table 2 :

 2 Synthesis results for a Xilinx Virtex6 device

	Design Approach radix	Resources LUT reg.	Performance cycles@period(ns)
				7,880	0	1@94.3
			2	7,966 7,299	1,523 2,689	11@9.6 17@5.7
				6,786	5,202	36@3.7
				4,871	0	1@57.9
		Ours	4	4,768 4,600	988 1,583	7@12.3 11@6.9
				4,853	3,106	22@3.8
	Ex. 1			4,210	0	1@44.4
				3,875*	0	1@62.2*
			8	5,307*	309	5@18.4*
				5,184*	499	8@10.4*
				4,707* 1,027	15@5.8*
				994	0	1@29.5
		FloPoCo	-	1,032	138	7@6.7
				1,147	335	19@5.3
			2	6,820	0	1@88.5
	Ex. 2	Ours	4 8	6,356 5,042	0 0	1@68.0 1@39.0
		FloPoCo	-	3,024	0	1@41.1
			2	2,944	0	1@67.0
			4	2,742	0	1@35.1
	Ex. 3	Ours	8 16	2,582 2,856 1,565*	0 0 0	1@33.1 1@31.2 1@29.0*
		FloPoCo	-	3,622	0	1@55.7
			2	19,564	0	1@139.6
	Ex. 4	Ours	4 8	23,052 21,179* 15,388*	0 0 0	1@92.5 1@131.5* 1@250.7*
			16	12,878*	0	1@76.9*
			32	3,909*	0	1@86.7*
		FloPoCo	-	20,494	0	1@139.9
			2	19,423	0	1@368.1
	Ex. 5	Ours	4 8	13,642 18,653	0 0	1@70.3 1@58.6
		FloPoCo	-	-	-	-

  . By the definition of 𝑃 𝑘+1 and 𝑄 𝑘+1 , we havemax𝑥∈𝑋 { |𝑓 (𝑥)𝑄 𝑘+1 (𝑥) -𝑃 𝑘+1 (𝑥)| -𝛿 𝑘 𝑄 𝑘+1 (𝑥) 𝑄 𝑘 (𝑥)This chain of inequalities tells us that𝛿 𝑘 𝑄 𝑘+1 (𝑥)/𝑄 𝑘 (𝑥) Δ(𝛿 𝑘 -𝛿)/𝑀, ∀𝑥 ∈ 𝑋(12)resulting in 𝑄 𝑘+1 (𝑥) > 0, ∀𝑥 ∈ 𝑋 and{|𝑓 (𝑥)𝑄 𝑘+1 (𝑥) -𝑃 𝑘+1 (𝑥)| -𝛿 𝑘 𝑄 𝑘+1 (𝑥)} /𝑄 𝑘 (𝑥) < 0, ∀𝑥 ∈ 𝑋, leading to |𝑓 (𝑥) -𝑃 𝑘+1 (𝑥)/𝑄 𝑘+1 (𝑥)| < 𝛿 𝑘 , ∀𝑥 ∈ 𝑋.The conclusion now follows by the definition of 𝛿 𝑘 . 𝛿 𝑘 → 𝛿 * as 𝑘 → ∞.Proof. By the previous theorem, we know that {𝛿 𝑘 } is a monotonically decreasing sequence bounded below by 0. It is thus convergent. Denote its limit by 𝛿 and assume that 𝛿 > 𝛿 * . Then, there exists 𝑅 = 𝑃 /𝑄 ∈ ℛ𝜇,𝜈 (𝑋) s.t.Let Δ and 𝑀 be as in the previous theorem. By the aforementioned chain of inequalities, we get 𝛿 𝑘 -Δ 𝑀 (𝛿 𝑘 -𝛿) 𝑄 𝑘 (𝑧 𝑘 ) 𝑄 𝑘+1 (𝑧 𝑘 ) ,

									}
		max𝑥∈𝑋	{	𝑄 𝑘 (𝑥) ⃒ 𝑓 (𝑥)𝑄(𝑥) -𝑃 (𝑥) ⃒	⃒ ⃒ -𝛿 𝑘 𝑄(𝑥)	}
		max𝑥∈𝑋	{[⃒ ⃒ ⃒ ⃒ 𝑓 (𝑥) -	𝑃 (𝑥) 𝑄(𝑥)	⃒ ⃒ ⃒ ⃒ -𝛿 𝑘	]	𝑄 𝑘 (𝑥) 𝑄(𝑥)	}
			max𝑥∈𝑋	{	(𝛿 -𝛿 𝑘 )	𝑄 𝑘 (𝑥) 𝑄(𝑥)	}
		= -(𝛿 𝑘 -𝛿) min𝑥∈𝑋	𝑄(𝑥) 𝑄 𝑘 (𝑥)
					-	Δ 𝑀	(𝛿 𝑘 -𝛿).
	Theorem 2.
					⃦ ⃦ 𝑓 -𝑅	⃦ ⃦ = 𝛿 < 𝛿.
	⃒ ⃒ ⃒ ⃒ 𝑓 (𝑥) -	𝑃 𝑘+1 (𝑥) 𝑄 𝑘+1 (𝑥)	⃒ ⃒ ⃒ ⃒ -𝛿 𝑘 -	Δ 𝑀	(𝛿 𝑘 -𝛿)	𝑄 𝑘 (𝑥) 𝑄 𝑘+1 (𝑥)	,	∀𝑥 ∈ 𝑋
	and
		𝛿 𝑘+1 -

  it follows from[START_REF] Cheney | Two new algorithms for rational approximation[END_REF] that there exists a positive integer 𝑘0 such that𝑄 𝑘 (𝑧 𝑘 ) 𝑄 𝑘+1 (𝑧 𝑘 ) 𝑐 |𝑋| , ∀𝑘 𝑘0.(16)Using (15) and (16), we have 𝑄 𝑘 (𝑧 𝑘 ) 𝑄 𝑘+1 (𝑧 𝑘 )

			∏ 𝑥∈𝑋,𝑥̸ =𝑧 𝑘	𝑄 𝑘 (𝑥) 𝑄 𝑘+1 (𝑥)	𝑐 |𝑋|	(	1 𝑐	) |𝑋|-1	= 𝑐,
	and		∏ 𝑥∈𝑋	𝑄 𝑘+1 (𝑥) 𝑄 𝑘 (𝑥)	1 𝑐	2.
									implies that
				lim 𝑘→∞	𝑄 𝑘 (𝑧 𝑘 ) 𝑄 𝑘+1 (𝑧 𝑘 )	= 0.	(14)
	From (12) we get
			𝑄 𝑘+1 (𝑥) 𝑄 𝑘 (𝑥)	Δ 𝑀	(𝛿 𝑘 -𝛿)/𝛿 𝑘 𝑐𝑄 𝑘 (𝑥),	∀𝑥 ∈ 𝑋,	(15)
	provided that 𝑐	Δ 𝑀 (1 -𝛿/𝛿)	Δ 𝑀 (𝛿 𝑘 -𝛿)/𝛿 𝑘 . Taking
			𝑐 = min	{	1 2	,	Δ 𝑀	(1 -𝛿/𝛿)	}	,

Acknowledgments

This work was partly supported by the FastRelax project of the French Agence Nationale de la Recherche, EPSRC (UK) under grants EP/K034448/1 and EP/P010040/1, the Royal Academy of Engineering and Imagination Technologies.


Mohamed Ait Hamza
email: mohamed.aithamza@edu.uiz.ac.ma

Hicham Lakhtar

Hafssa Tazi

Abdelmajid Moukhli

Odile Fossati-Gaschignard

Lucie Miche

Sevastianos Roussos

Zahra Ferji

Abdelhamid El Mousadik

Thierry Mateille

Mohamed Aït Hamza

Lucie Miché

Sebastianos Roussos

Abdelhamid El Mousadik

Hassan Boubaker
email: hassanboubaker@yahoo.fr

Diversity of nematophagous fungi in Moroccan olive nurseries: Highlighting prey-predator interactions and efficient strains against root-knot nematodes

Keywords: Culture substrate, Ecology, Fungal antagonists, Olive tree, Root-knot nematode

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The olive tree (Olea europaea L.) is widely distributed throughout the Mediterranean Basin.

In Morocco alone, the olive-growing area is estimated at 650,000 ha and produces 120,000 tons of oil and 150,000 tons of table olives per year (Ministry of Agricultural and Marine Fisheries, 2009). Planting material comes from several nurseries distributed throughout the olive-producing areas and olive plantlets are certified free of pathogens (e.g., Verticillium dahlia) and parasites (e.g., plant-parasitic nematodes). Nevertheless, standard health practices are not applied in all nurseries and seasonal and informal nurseries coexist.

Several plant-parasitic nematode species (PPN) have been found to be associated with olive trees [START_REF] Castillo | Plant-parasitic nematodes attacking olive trees and their management[END_REF][START_REF] Ali | Plant-parasitic nematodes associated with olive tree (Olea europaea L.) with a focus on the Mediterranean Basin: A review[END_REF][START_REF] Hamza | Plant-parasitic nematodes associated with olive tree in Southern Morocco[END_REF], including root-knot nematodes (RKN, Meloidogyne spp.) that were shown to impact olive growth [START_REF] Afshar | Effects of the root-knot nematodes Meloidogyne incognita and M. javanica on olive plants growth in glasshouse conditions[END_REF]. RKN have been reported to be major pests on olive trees, mainly in nurseries where irrigation conditions are favorable to their multiplication [START_REF] Nico | Incidence and population density of plant-parasitic nematodes associated with olive planting stocks at nurseries in southern Spain[END_REF][START_REF] Sanei | Incidence of plant-parasitic nematodes associated with olive planting stocks at nurseries in northern Iran[END_REF]. Moreover, most RKN infestations in olive orchards result from contaminated plant material produced in uncertified nurseries [START_REF] Nico | Incidence and population density of plant-parasitic nematodes associated with olive planting stocks at nurseries in southern Spain[END_REF]. To protect nurseries, vegetative propagation of olive is a crucial step to prevent plant pest dispersal.

While some industrial olive-producing countries use synthetic substrates to cultivate olive trees, other countries use more traditional olive cultivation techniques, planting roots in soil substrates. In countries that have implemented large-scale planting programs (e.g., 1,220,000 ha in the Morocco Green Plan, 2015), the "nursery risk" is a concern, with two main consequences: (i) the weakening of olive plants before their transplantation into orchards; and

(ii) the introduction of pathogens into orchards. The nematode risk is now controlled in nurseries, either by multiplication of resistant varieties [START_REF] Palomares-Rius | Nematode community populations in the rhizosphere of cultivated olive differs according to the plant genotype[END_REF] or by inoculating microbial antagonists in the substrates, such as mycorrhizae [START_REF] Castillo | Protection of olive planting stocks against parasitism of root-knot nematodes by arbuscular mycorrhizal fungi[END_REF].

Thus, in view of the large-scale development of olive cultivation in Morocco, sanitization of the substrates must be a priority. Some previous studies highlighted that certain substrates can be suppressive depending on their composition, either because of microbial antagonisms or of toxicity. For example, the addition of forest residues provides organic matter that could contribute to PPN suppression by increasing soil microflora [START_REF] Rodriguez-Kabana | Biological control of nematodes: soil amendments and microbial antagonists[END_REF].

Composted dry cork, especially in nurseries, could be effective for M. incognita suppression due to the toxicity of the products released (e.g., ammonia, phenolic compounds) [START_REF] Nico | Control of root-knot nematodes by composted agro-industrial wastes in potting mixtures[END_REF]. Moreover, the introduction of microbial antagonists into substrates would strengthen the sustainability of preventive techniques and help introduce biocontrol agents into orchards.

Microbial antagonists of PPN include nematophagous fungi (NF), raising expectations for their use in integrated pest management [START_REF] Waller | From discovery to development: current industry perspectives fo the development on novel methods of helminth control in livestock[END_REF][START_REF] Larsen | Biological control of nematodes parasites in sheep[END_REF][START_REF] Maingi | Control of gastrointestinal nematodes in goats on pastures in South Africa using nematophagous fungi Duddingtonia flagrans and selective anthelmintic treatments[END_REF].

These fungi act by antibiosis, parasitism or predation [START_REF] Imerglik | Recherches préliminaires sur la spécificité du piégage des nématodes par des hyphomycètes prédateurs[END_REF][START_REF] Gaspard | Nematophagous fungi associated with Tylenchulus semipenetrans and the citrus rhizosphere[END_REF]. Fungi trap nematodes using adhesive networks or buttons. Some NF capture their prey with constrictor rings before strangling them and others with hyphae networks that produce sticky compounds. The specific enzymes then digest the nematodes. Recognition and attachment of the mycelium to the nematode cuticle are mainly due to compatible glycoproteins such as lectins [START_REF] Nordbring-Hertz | Nematophagous fungi[END_REF].

Nevertheless, edaphic factors such as soil pH, temperature, moisture and structure influence their efficiency [START_REF] Brown | Principles and practice of nematode control in crops[END_REF]. Fungal diversity can increase in stressed or polluted sites [START_REF] Korol | Recombination variability and evolution[END_REF] in order to adapt to changing environments [START_REF] West | The fourth dimension of life: fractal geometry and allometric scaling of organisms[END_REF]. In Arthrobotrys oligospora, genetic diversity varies according to environmental conditions.

Their strains could generate recombinant genotypes by crossing with native strains, thus enhancing their environmental adaptability and parasitizing ability [START_REF] Zhang | Genetic diversity and recombination in natural populations of the nematode-trapping fungus Arthrobotrys oligospora from China[END_REF].

However, the reproductive strategies of fungi like A. oligospora can change depending on the population. According to [START_REF] Cook | Making greater use of introduced microorganisms for biological control of plant pathogens[END_REF], the search for microorganisms from the rhizospheric soil of a specific crop could lead to the isolation of effective antagonists against pathogens that could be adapted to the plant species as well as to particular environmental conditions.

Considering that olive seedlings in Morocco are generally cultivated in non-sanitized substrates consisting of soil material from various habitat origins, the recovery of native NF from substrates may be of great interest in order to develop RKN biocontrol agents adapted to nursery conditions. In this context, this work aims to (i) evaluate the nematode populations in substrates from olive nurseries in Morocco; (ii) isolate and characterize NF able to control RKN (Meloidogyne spp.); (iii) determine the occurrence and diversity of nematode-associated NF and discuss prey-predator interactions; and (iv) investigate their in vitro predatory potential towards M. javanica juveniles and eggs, and their usefulness as biocontrol agents for olive protection in nurseries and in orchards.

Materials and methods

Site description and olive plantlet sampling

Soil samples were collected in spring 2013 and 2014 from 25 commercial olive nurseries located in the main olive-producing areas in Morocco (Fig. 1 and Table 1): the Jbala, Guerouane, Haouz and Souss regions. The nurseries were selected for their plantlet production, the cultivars grown and the rearing substrates used. In each nursery and for each variety, five olive plantlets (Olea europaea subsp. europaea)

growing in plastic bags were sampled. Information about the origin and the preparation of the growth substrates and about the cultivars was recorded. A total of 305 olive plantlets were collected and maintained in the laboratory and kept under greenhouse conditions (12 h light at 25 °C; 12 h dark at 20°C).

Nematode extraction and quantification

A 250-cm 3 substrate subsample was removed from the rhizosphere of each olive plantlet and used for nematode extraction using the [START_REF] Oostenbrink | Estimating nematode populations by some selected methods[END_REF] elutriation procedure (ISO 23611-4). Free-living nematodes (FLN) and plant-parasitic nematodes (PPN), which specifically exhibit a stylet for plant-cell feeding, were enumerated in 5-cm 3 counting chambers [START_REF] Merny | Les techniques d'échantillonnage des peuplements de nématodes dans le sol[END_REF] under a stereomicroscope (×60 magnification). Nematode population levels were expressed as the number of individuals per dm 3 of fresh substrate. Among PPN, rootknot nematodes belonging to the Meloidogyne genus [START_REF] Mai | Plant-parasitic nematodes: a pictorial key to genera, 5th edn[END_REF] were counted.

Isolation of nematophagous fungi

NF were isolated from solid substrate samples infested with RKN using the soil sprinkling technique [START_REF] Duddington | Notes on the technique of handling predacious fungi[END_REF] as modified by [START_REF] Santos | Detection and ecology of nematophagous fungi from Brazil soils[END_REF]. The direct soil powdering on media was preferable to the modified Baermann method (aqueous soil suspensions) [START_REF] Hernández-Chavarría | A simple modification of the Baermann method for diagnosis of strongyloidiasis[END_REF] because it provided more concentrated material for fungal isolation. Moreover, isolation success was increased by sprinkling cold soil (stored at 4°C) on culture medium at 37 to 40°C [START_REF] Davet | Detection and isolation of soil fungi[END_REF]. Substrate aliquots from each olive plantlet were spread on a tray to be air-dried. One gram was then sprinkled on the surface of Petri dishes containing water-agar (WA 2% w/v) supplemented with antibiotics (0.05% streptomycin-sulphate and 0.05% chloramphenicol). Three replicates were done per olive plantlet. A 1-mL suspension containing approximately 3,000 M. javanica second-stage juveniles (J2) and 10 eggs produced in the laboratory was added as fungal bait, according to the procedure used by [START_REF] Drechsler | Some hyphomycetes parasitic on free-living terricokms nematodes[END_REF]

Morphological and molecular characterization of the fungi

The initial identification of the fungi was based on colony morphology and microscopic characteristics. Slide sub-cultures from pure NF cultures were observed under a dissecting microscope (up to x100 magnification). Genera and species were assigned according to specialized morpho-taxonomical keys [START_REF] Cooke | Nematode-trapping species of Dactylella and Monacrosporium[END_REF][START_REF] Haard | Taxonomic studies on the genus Arthrobotrys Corda[END_REF][START_REF] Barron | The nematode-destroying fungi: Canadian Biological Publications Ltd[END_REF][START_REF] Yu | Taxonomy of nematode-trapping fungi from Orbiliaceae, Ascomycota[END_REF][START_REF] Philip | Nematophagous fungi: Guide by Philip Jacobs, BRICVersion online[END_REF].

Sequence analyses of the ITS (internal transcribed spacer) region in the ribosomal RNA gene cluster were performed to confirm the identity of the NF species. The DNA was extracted from 50-200 mg of mycelium (fresh weight) using the NucleoSpin ® Plant II Genomic DNA Purification Kit (Promega ® ) according to the manufacturer's instructions. PCR reaction was performed according to the method described by [START_REF] White | Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics[END_REF]: the ITS rDNA gene cluster was amplified using the primers ITS1 (5´TCC GTA GGT GAA CCT GCG G 3´) and ITS4 (5´TCC TCC GCT TAT TGA TAT GC 3´). The PCR amplification was carried out using the GeneAmpR PCR System 9700 (Applied Biosystems ® ). Twenty µl of reaction mixture contained 2 µl (10 ng) of template DNA, 1 µl of each ITS1 and ITS4 primer (10 mM), 4 µl PCR buffer, 2.4 µl MgCl2 (25 mM), 0.6 µl dNTPs (10 mM), 0.1 µl BSA (0.1 mg/ml) and 0.2 units of GoTaq ® DNA polymerase.

PCR cycling conditions consisted of an initial denaturation step at 94°C for 3 min, 35 cycles of denaturation at 94°C for 30 sec, annealing at 55°C for 30 sec and elongation at 72°C for 10 min. PCR products were checked for length, quality and quantity by agarose gel electrophoresis (1% (w/v) in 0.5x Tris-Acetate-EDTA (TAE).

PCR products were sequenced from both ends by Eurofins MWG GmbH (Ebersberg, Germany), using the same ITS primers. CHROMAS LITE v2.1.1 (Technelysium Pty Ltd.)

software was used to edit and assemble DNA sequences. BLAST similarity searches were performed in the non-redundant nucleotide database of GenBank [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search programs[END_REF] to identify/verify species or genus affiliation of collected isolates. Sequences were aligned with ITS sequences of reference strains obtained from GenBank. Subsequently, the alignment was used to perform the phylogenetic tree with PHYLOGENY.FR [START_REF] Dereeper | Phylogeny. fr: robust phylogenetic analysis for the non-specialist[END_REF] and Mega 6 [START_REF] Tamura | MEGA4: molecular evolutionary genetics analysis (MEGA) software version 4.0[END_REF] using the neighbor-joining method [START_REF] Saitou | The neighbor-joining method: a new method for reconstructing phylogenetic trees[END_REF] and the Jukes-Cantor correct distance model [START_REF] Jukes | Evolution of protein molecules[END_REF]. The phylogenetic tree was obtained from data using one of three equally parsimonious trees through 1,000 bootstrap replicates [START_REF] Felsenstein | Confidence intervals on phylogenetics:an approach using bootstrap[END_REF] with a heuristic search consisting of 10 random-addition replicates for each bootstrap replicate.

Diversity of the fungal communities

Four diversity indices were calculated to assess NF communities: total number of isolates found in each soil sample (N); richness (S = number of species in the community); Shannon-Wiener local diversity index (H'=-Σ (pi.lnpi), where pi is the proportion of isolates with the species i); and evenness (E = H'/ln(S)), which quantifies the numerical equality of populations in communities.

Pathogenicity analyses

An olive population of the RKN M. javanica (detected in 72% of the nurseries surveyed and dominant in orchards) [START_REF] Ali | Trend to explain the distribution of root-knot nematodes Meloidogyne spp. associated with olive trees in Morocco[END_REF] was reared on RKN-susceptible tomato (cv. Roma) in a greenhouse (12 h light at 25°C; 12 h dark at 20°C). Fungal strains that exhibited trapping, adhesive or encysting organs were sub-cultured on WA in 9-cm-diameter Petri dishes. One week later, 100 second-stage juveniles (J2) of M. javanica (the only free form in the soil)

were washed five times with 0.05% streptomycin-sulfate in sterilized distilled water, introduced into each fungal sub-culture and maintained at 25°C in darkness. Three Petri dishes per fungal strain were considered as replicates. Fungal predation structures were observed and predated/dead nematodes were counted after four days under a microscope (x100 magnification). Similar procedures were used to study predation of M. javanica eggs by specialized egg parasites such as Paecilomyces and Pochonia fungi. Dishes with nematodes but without fungi were considered as control replicates.

Statistical analyses

Mean values were analyzed by one-way ANOVA and Kruskall-Wallis tests were used for all pair-wise multiple comparisons. NF community patterns were explored through a Principal

Component Analysis (PCA) of the diversity indices. Region grouping was tested using Monte-Carlo tests on PCA eigenvalues (randtest, ade4). Calculations were performed and graphics prepared using R language (readxl, base and ade4 packages) (R Development Core Team, 2011;[START_REF] Chessel | The ade4 package. I. One-table methods[END_REF][START_REF] Dray | The ade4 package: implementing the duality diagram for ecologists[END_REF]Dufour, 2007, Wickham, 2016), with a level of significance = 0.05). A rarefaction regression was used to analyze the dependence between root-knot nematodes (Meloidogyne spp.) and nematophagous fungi.

Results

Nematofauna

A non-significant gradient was present in plant-parasitic-nematodes (PPN), with a north-south increase of the population levels (Table 2). A parallel significant gradient was revealed in root-knot nematodes (RKN, Meloidogyne spp.). Free-living nematodes (FLN) were three times more abundant in the Souss olive nurseries than in the other regions. The ratios between FLN and PPN were different between all regions, the highest ratio being found in the Souss region.

Nematophagous species and phylogenetic diversity

Several fungal isolates were recovered from the 305 soil samples examined. Observation of characteristic conidia and traps around dead M. javanica juveniles revealed 149 soil samples positive for NF. Morphological identification using microculture techniques revealed 73 NF strains belonging to 11 genera.

In order to confirm the characterization of the fungal strains, the ITS regions of rDNA were sequenced. Five species (Catenaria anguillulae, Nematoctonus leiosporus, Haptoglossa heterospora, Dactylaria sp. and Monacrosporium microphoides) were excluded from the sequencing because it was not possible to purify them. The BLAST test showed that the ITS sequences of all sequenced strains were at least 99% similar to the corresponding GenBank reference sequences (Table 3). The phylogenetic analysis including ITS sequences of the Moroccan NF isolates and 16 reference sequences of identified close relatives revealed five distinct clusters (Figs. 2 A-E

Diversity patterns

Because of their scarcity (only one nursery surveyed), the substrate samples from the Jbala region were excluded from the dataset prior to running analyses. Richness (S) and local diversity (H') were correlated to the PC1 axis, while the PC2 axis was related to numbers (N) and to evenness (E) of fungal isolates on its positive and negative sides, respectively (Fig. 3A). Region grouping was significant in the whole analysis and on the PC2 axis (Fig. 3B).

Isolates were more numerous in the Souss region, whereas fungal communities were more numerically alike in the Guerouane region. The Haouz region had lower NF richness and diversity (non-significant PC1 coordinates). Indices were affected by the north-south distribution of the nurseries (Figs. 3C-3E-3D).

The percentage of samples with NF was lower in the Haouz region than elsewhere (Table 5).

Both trapping and endoparasitic fungi occurred more often in the Souss region than in the others regions, and no endoparasitic species were found in the Jbala region.

The rarefaction regression established between the occurrence of RKN (Meloidogyne spp.)

and the occurrence of NF (Fig. 4) indicated a significant positive correlation, regardless of the region sampled.

3.4.

In vitro efficiency of the nematophagous strains NF were distinguished according to their ability to kill M. javanica J2s (Fig. 5). Talaromyces assiutensis killed all juveniles in all replicates. The Orbiliaceae species (Arthrobotrys spp., Dreschlerella spp., Monacrosporium spp.) were efficient against M. javanica since they trap 50 to 80% of the J2s using adhesives networks, buttons, constricting rings and hyphae networks. Paecilomyces and Trichoderma strains killed 30 to 50% of the J2s. Fusarium oxysporum strains were less efficient (less than 20% of dead J2s). P. lilacinus and P.

chlamydosporia strains infected all the M. javanica eggs.

Discussion

Our first objective was to evaluate nematode populations in substrates from olive nurseries. In all regions of Morocco, PPN abundance was greater than 1.4 nematodes/dm 3 of soil and FLN abundance greater than 1.9. As a comparison, [START_REF] Hamza | Plant-parasitic nematodes associated with olive tree in Southern Morocco[END_REF] found 0.2 to 5.1 PPN/dm 3 of soil and 0.3 to 4.3 FLN/dm 3 of soil in 23 Souss and Haouz orchards. It can thus be hypothesized that the multiplication of nematode populations may be boosted by acidity and by hydrophilic and non-degraded organic matter [START_REF] Neher | Nematode communities in soils of four farm cropping management systems[END_REF][START_REF] Manlay | Relationships between abiotic and biotic soil properties during fallow periods in the sudanian zone of Senegal[END_REF][START_REF] Ou | Vertical distribution of soil nematodes under different land use types in an aquic brown soil[END_REF][START_REF] Mcsorley | Effect of disturbances on trophic groups in soil nematode assemblages[END_REF]. The usual dominance of FLN in olive soils could be due to the origin of organic substrates (mountain soils, peat and manure) used in most nurseries [START_REF] Castillo | Protection of olive planting stocks against parasitism of root-knot nematodes by arbuscular mycorrhizal fungi[END_REF]. FLN are known to dominate in soil substrates not yet used for agriculture [START_REF] Hillocks | Associations between soilborne pathogens and other soil-inhabiting microorganisms[END_REF]. The FLN dominance in the Souss region, one of the most intensively cultured areas in Morocco, was thus unexpected.

The FLN/PPN ratios were balanced in all regions except in the Souss where most of the soil substrates come from cropped areas (vegetables and citrus fruits). Because of the sandy texture of the soils in the Souss region, culture practices include high amounts of organic matter (especially cattle manure or tomato leaf compost), leading to an increase in FLN populations and to a decrease in PPN populations because organic matter is unsuitable for PPN [START_REF] Clark | Agronomic, economic, and environmental comparison of pest management in conventional and alternative tomato and corn systems in northern California[END_REF][START_REF] Hominick | Nematodes[END_REF][START_REF] Hu | Abundance and diversity of soil nematodes as influenced by different types of organic manure[END_REF]. These mechanisms may explain the FLN/PPN ratio that is twice as high in the Souss region compared to the other regions The high percentages of RKN (Meloidogyne spp.) in Souss olive nurseries may be justified by the dominance throughout the area of vegetable crops that are highly susceptible to these nematodes [START_REF] Sikora | Plant parasitic nematodes in subtropical and tropical agriculture[END_REF][START_REF] Netscher | Les nématodes parasites des cultures maraîchères au Sénégal[END_REF].

Our second objective was to isolate and characterize nematophagous fungi able to control RKN (Meloidogyne spp.). The success of fungal strain isolation from soils may be correlated with soil temperature, moisture and organic matter content [START_REF] Akhtar | Roles of organic soil amendments and soil organisms in the biological control of plant-parasitic nematodes: a review[END_REF][START_REF] Cayrol | La lutte biologique contre les nématodes phytoparasites[END_REF]. The greatest difficulty encountered was in isolating pure NF strains because of the rapid growth of plant-pathogenic fungi and saprophytes such as Fusarium, Alternaria and Aspergillus species. An increased pH of culture media may prevent the development of other microorganisms [START_REF] Gardner | Production of chlamydospores of the nematode-trapping Duddingtonia flagrans in shake flask culture[END_REF]. Moreover, we found that direct cold soil powdering on hot media (37 to 40°C), like in [START_REF] Kelly | Screening for the presence of nematophagous fungi collected from Irish sheep pastures[END_REF], was preferable to aqueous soil suspensions because it provided more concentrated fungal material [START_REF] Davet | Detection and isolation of soil fungi[END_REF][START_REF] Hernández-Chavarría | A simple modification of the Baermann method for diagnosis of strongyloidiasis[END_REF]. Dispersing agents were avoided because they inhibited the growth of some NF strains [START_REF] Davet | Detection and isolation of soil fungi[END_REF].

This survey provided native NF from Morocco for the first time. The large number of soil samples (305) allowed the detection of numerous NF strains, whereas former studies revealed one strain at best [START_REF] Bridge | Soil fungi: diversity and detection[END_REF]. The strains detected as nematophagous possessed different modes of action: adhesive networks, constricting rings, hyphal tips, adhesive conidia and mycotoxins [START_REF] Imerglik | Recherches préliminaires sur la spécificité du piégage des nématodes par des hyphomycètes prédateurs[END_REF][START_REF] Gaspard | Nematophagous fungi associated with Tylenchulus semipenetrans and the citrus rhizosphere[END_REF]. Despite this, ITS rDNA sequences could not be used to identify some NF strains due to their small size (Pochonia chlamydosporia, for example), BLAST tests were useful to confirm the morphological characterization of other strains. The integrated taxonomical analysis (morphological and molecular) of the NF provided a pool of 28 strains belonging to 19 species. The phylogenic analysis revealed that the ITS rDNA gene was able to distinguish species in a genus group (such as Trichoderma) but did not fully discriminate species belonging to the family Orbiliaceae, indicating the taxonomic proximity of Dreschlerella, Arthrobotrys and Monacrosporium species. Therefore, an integrative molecular analysis should be developed with other molecular markers [START_REF] White | Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics[END_REF] in order to improve NF identification.

Our third objective was to determine the occurrence and diversity of nematode-associated NF.

All NF species except A. oligospora were detected in the Souss region where they were more dominant than in the other regions. This high diversity of NF might be due to the multiple habitat origins of the components (mountain, riverbank and field soils, cattle manure, plant compost, etc.) used to make substrates for root olive plantlets. We hypothesize that the microbial richness detected in the Souss region corresponds to the high endemic plant diversity characteristic of the Macaronesian region [START_REF] Médail | Glacial refugia influence plant diversity patterns in the Mediterranean Basin[END_REF][START_REF] Msanda | Biodiversité et biogéographie de l'arganeraie marocaine[END_REF]. This could also be due to higher organic matter concentration in the soils (i.e., saprophytic substrate for fungi that induce the formation and the activity of trapping structures) and to high PPN levels (i.e., parasitic substrates for fungi) [START_REF] Den Belder | Capture of plant-parasitic nematodes by an adhesive hyphae forming isolate of Arthrobotrys oligospora and some other nematode-trapping fungi[END_REF][START_REF] Singh | Evaluation of biocontrol potential of Arthrobotrys oligospora against Meloidogyne graminicola and Rhizoctonia solani in Rice (Oryza sativa L.)[END_REF].

The number of NF isolates globally increased southwards, whereas evenness decreased (PCA2 axis). It is known that different thermal regimes affect soil microbial diversity [START_REF] Bridge | Soil fungi: diversity and detection[END_REF]. The absence of seasonality and the higher minimum temperatures in the Souss region may explain these developments where irrigation counteracts the arid to semiarid climate. The high occurrence of NF and the high evenness detected in the Guerouane region may be linked to a more continental climate. More than half of the 36 P. lilacinus strains came from the Souss region where the NF richness was the highest but the evenness the lowest. In a restricted soil area, the abundance of P. lilacinus may cause the rarity of the other species. For a constant number of species, maximal diversity is achieved when species have an even distribution. The Haouz, Jbala and Guerouane regions were characterized by less PPN and FLN than in the Souss region, explaining the lower NF richness in those regions.

Our fourth objective was to investigate the in vitro predation of the NF against M. javanica juveniles and eggs. Talaromyces assiutensis (strain UIZFSA-31), whose nematode predation was previously unknown, killed all the M. javanica juveniles in four days. The mode of action of T. assiutensis remains unknown but we hypothesize that the strain may produce specific mycotoxins. In vitro trapping tests prove that the Orbiliaceae species (Arthrobotrys spp., Dreschlerella spp., Monacrosporium spp.) were able to trap M. javanica juveniles. The predatory capacity of A. oligospora was similar to data found in the literature [START_REF] Singh | Evaluation of biocontrol potential of Arthrobotrys oligospora against Meloidogyne graminicola and Rhizoctonia solani in Rice (Oryza sativa L.)[END_REF]. The mechanisms involved during predation are well known [START_REF] Imerglik | Recherches préliminaires sur la spécificité du piégage des nématodes par des hyphomycètes prédateurs[END_REF][START_REF] Gaspard | Nematophagous fungi associated with Tylenchulus semipenetrans and the citrus rhizosphere[END_REF]. Recognition and attachment of the mycelium to the cuticle of the RKN juvenile is mainly due to compatible glycoproteins, e.g., lectins [START_REF] Duponnois | Effect of different west african species and strains of Arthrobotrys nematophagous fungi on Meloidogyne species[END_REF].

Trichoderma species are recognizsed as control agents against nematodes, and various mechanisms have been proposed to explain nematode killing, including antibiosis and enzymatic hydrolysis [START_REF] Sivan | Microbial control of plant diseases[END_REF][START_REF] Elad | Biological control of foliar pathogens by means of Trichoderma harzianum and potential modes of action[END_REF]. [START_REF] Thomas | Studies on the parasitism of Globodera rostochiensis by Trichoderma harzianum using low temperature scanning electron microscopy[END_REF] demonstrated direct interactions between T. harzianum and the cyst nematode Globodera rostochiensis: the fungus penetrated eggs in the cysts, leading to the death of the juveniles [START_REF] Sharon | Biological control of the root-knot nematode Meloidogyne javanica by Trichoderma harzianum[END_REF].

Precise information on the mechanisms involved is very limited and this misunderstanding has hindered the selection of active strains and the development of improved biocontrol methods. Paecilomyces lilacinus and Pochonia chlamydosporia are especially powerful nematode egg parasites [START_REF] Cayrol | Etude préliminaire sur les possibilités d'utilisation des champignons parasites comme agents de lutte biologique[END_REF][START_REF] Irving | Variation between strains of the nematophagous fungus, Verticillium chlamydosporium Goddard. II. Factors affecting parasitism of cyst nematode eggs[END_REF], which could explain the relatively low predation rate obtained on juveniles. These species can also act on the movement of infested nematode juveniles via a paralytic toxin [START_REF] Cayrol | Study of the nematicidal properties of the culture filtrate of the nematophagous fungus Paecilomyces lilacinus[END_REF], purified and identified as acetic acid by [START_REF] Djian | Acetic acid: a selective nematicidal metabolite from culture filtrates of Paecilomyces lilacinus (Thom) Samson and Trichoderma longibrachiatum Rifai[END_REF]. This molecule is abundantly produced during fungal growth in liquid medium (ibid). Fusarium oxysporum strains, well-known plant pathogens, exhibited lower predation rates, but some studies revealed that they are partly able to kill nematodes by producing toxic sulfuric heterocycles (fusarenone and moniliformine) [START_REF] Ciancio | Nematicidal effects of some Fusarium toxins[END_REF][START_REF] Cayrol | La lutte biologique contre les nématodes phytoparasites[END_REF].

Moreover, the significant co-occurrence of RKN and NF highlighted by the rarefaction curve suggests a close interaction between prey (nematodes) and predators (fungi) in olive nurseries, as described by Lotka-Volterra models [START_REF] Barbosa | Ecology of predator-prey interactions[END_REF], probably because of parasitism induction by nematodes [START_REF] Jansson | Interactions between nematophagous fungi and plant-parasitic nematodes: attraction, induction of trap formation and capture[END_REF]. Even though the diversity of the NF detected in the four regions was variable at species and population levels, we may expect symmetric dynamics [START_REF] Marrow | Evolutionary instability in predator-prey systems[END_REF]) due to possible co-evolutionary processes involved between nematodes and fungi. Such processes occur between competitive organisms in an ever-changing environment, as described by the "Red Queen" hypothesis [START_REF] Van Valen | A new evolutionary law[END_REF][START_REF] Dawkins | Arms races between and within species[END_REF].

Conclusion

Olive nursery solid substrates are infested with PPN, including Meloidogyne species. Various predatory fungi were able to kill RKN. Consequently, before selecting NF strains as candidates for biocontrol, studies must be extended to a wider range of Meloidogyne species and populations present on olive trees in Morocco in order to verify their specificity.

Indigenous NF strains were recovered from different substrates (different habitats), making it possible to undertake more research in order to understand the specificity of prey-predator interactions with more diverse PPN species. Predation efficiency in different cropping systems and in varied soil environments should also be explored in the future.

Legends for figures and tables

Figure 1 Distribution of the olive nurseries surveyed in Morocco. See Table 1 for more information.

  Olive plants are grown in 2-to-3 liter plastic bags filled with solid substrates from different sources (alluvial sandy soils, forest soils, loamy open-field soils) supplemented with different proportions of sand, peat fertilizer and animal manure. Plants are first grown in plastic greenhouses and then outside. They are watered by sprinklers and fertilized with Osmocote® (Everris Company™).

  ): Purpureocillium lilacinum, Trichoderma, Fusarium oxysporum, Talaromyces and Arthrobotrys-Dreschslerella-Monacrosporium (family Orbiliaceae). The combination of the morphological and the molecular analyses established that the isolated NF belonged to 19 species, eight families and six orders (

Figure 2 Figure 3

 23 Figure2Neighbor-joining tree inferred from ITS rDNA sequences of the nematophagous

Figure 4

 4 Figure 4 Rarefaction curve for root-knot nematodes (Meloidogyne spp.) and nematophagous

Figure 5 .

 5 Figure 5. In vitro efficiency of nematophagous fungal strains: percentage of dead M. javanica

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Neighbor-joining tree inferred from ITS rDNA sequences of the nematophagous

Figure 5 .

 5 Figure 5. In vitro efficiency of nematophagous fungal stains: percentage of dead M. javanica

Table 4

 4 Paecilomyces lilacinusthat represented 36 strains (50% of the strains). Half of the P. lilacinus strains detected were isolated from the Souss region, and then gradually decreased northwards (12 in Haouz, four in Guerouane, two in Jbala). Arthrobotrys brochopaga, A. scaphoides, Monacrosporium thaumasium, P. lilacinus, F. oxysporum, T. harzianum, T. asperellum and T. longibrachiatum were encountered in the four regions. A. oligospora, Dactylaria sp., Haptoglossa heterospora, Monacrosporium microscaphoides, Nematoctunus leiosporus and Talaromyces assiutensis were very rare (one isolate in one region). All the fungal species except A. oligospora were detected in the Souss region. Only eight species were detected in the Jbala region. The Guerouane and Haouz regions hosted 11 and 12 species, respectively.

). Arthrobotrys was the most diversified genus with five species, followed by Trichoderma (three species). The abundance of fungal species was low (less than five strains per species) except for

Table 1

 1 Location and characteristics of the Moroccan olive nurseries surveyed. Number of samples for each cultivar in each geographic region.

Table 2

 2 Average density (number of nematodes/dm 3 of soil) of plant-parasitic nematodes

	(PPN) and free-living nematodes (FLN), percentages of samples infested with root-knot
	nematodes (RKN) and FLN/PPN ratios (a-d indicate significant groups, P < 0.05).

Table 3

 3 BLAST results of ITS rDNA sequences of the nematophagous fungi isolated.

Table 4

 4 Nematophagous fungi associated with the Moroccan olive nurseries surveyed.

Table 5

 5 Functional diversity of the nematophagous fungi (a-d indicate significant groups, P < 0.05).

	639

Table 1 640

 1 Location and characteristics of the Moroccan olive nurseries. Number of samples for each cultivar in each geographic region.

	641
	Geographic	Location	Climate	City	No. of	Main habitat origin	Olive cultivar	No. of samples
	region				nurseries	of the substrates
	Jbala	South-west face	Sub humid climate,	Ouazzane	1	Clay marls,	Picholine marocaine	5
		of the Rif Mountains	temperature from -6°			sand, forest soil	Haouzia	5
			to 32°C			and topsoil	Menara	5
	Guerouane	Sais Plateau, between	More continental	Meknes	4	Yellow sand,	Picholine marocaine
		the Middle Atlas	climate, temperature			topsoil,	Haouzia
		to the south and the Rif	from -10° to 45°C			mature manure	Menara
		Mountains to the north				and local	Arbequina	5
						compost	Arbosana	5
							Picual
							Picholine Languedoc	5
	Haouz	Northern slope	Semiarid climate,	Marrakech	5	Clay marls,	Picholine marocaine
		of the High Atlas	average temperatures			sand,	Arbequina
		Mountains	from -6° to 49.6°C			forest soil,	Haouzia
						mountain soil	Menara
						and topsoil	Picholine Languedoc
							Arbosana	5
				El Kelaa des Sraghna	3	Forest soil	Picholine marocaine
						and topsoil	Picholine Languedoc
							Menara
							Haouzia
				Sidi Abdellah Ghiat	1	Soil, clay and sand	Picholine marocaine	5
	Souss	On the southern slope	Arid constant climate,	Agadir	8	Sand, topsoil	Picholine marocaine
		of the High Atlas	sunshine > 340 days a			and peat moss	Haouzia
		Mountains	year, average				Menara	5
			temperatures from 14°	Khmiss Aït Amira	2	Topsoil, peat	Picholine marocaine
			to 25°C			and manure
				Biougra	1	Peat, soil and perlite	Menara	5

Table 2 642

 2 Average density (number of nematodes/dm 3 of soil) of plant-parasitic nematodes (PPN) and

	643
	644	free-living nematodes (FLN), percentages of samples infested with root-knot nematodes
	645	(RKN) and FLN/PPN ratios (a-d indicate significant groups, P < 0.05).
		Regions	No. of PPN	% of samples	No. of FLN	FLN/PPN
			(/dm 3 of soil)	infested with RKN	(/dm 3 of soil)
		Jbala	1,441	16.0 d	2,220 b	1.54 b
		Guerouane	1,527	20.0 c	1,914 b	1.25 c
		Haouz	2,003	58.1 b	2,081 b	1.04 d
		Souss	2,395	76.4 a	6,194 a	2.59 a
		P-value	0.153	0.000	0.000	0.000

Table 3 646

 3 BLAST results of ITS rDNA sequences of the nematophagous fungi isolated. Distribution of the olive nurseries surveyed in Morocco. See Table1for more information.

	647
	Strains	GenBank reference strains

Acknowledgements

This research was supported by a Ph.D. grant from the "Institut de Recherche pour le Développement" (Marseille, France). It was also funded by the PESTOLIVE project: Contribution of olive history for the management of soil-borne parasites in the Mediterranean Basin (ARIMNet action KBBE 219262), and by the BIONEMAR project: Development of fungal bionematicides for organic production in Morocco (PHC-Toubkal action 054/SVS/13).

Conflict of interest statement

We declare that we have no conflict of interest.

Table 4

Nematophagous fungi associated with the Moroccan olive nurseries surveyed.

Table 5

Functional diversity of the nematophagous fungi (a-d indicate significant groups, P < 0.05).

Region