Simon Pontie

Alban Bourge

Adrien Prost-Boucle

Paolo Maistri

Olivier Muller

Régis Leveugle

Frédéric Rousseau

HLS-Based Methodology for Fast Iterative Development Applied to Elliptic Curve Arithmetic

High-Level Synthesis (HLS) is used by hardware developers to achieve higher abstraction in circuit descriptions. In order to shorten the hardware development time via HLS, we present an adjustment of the Iterative and Incremental Design (IID) methodology, frequently used in software development. In particular, our methodology is relevant for the development of applications with unusual complexity: the method was applied here to the development of large modular arithmetic, commonly used for cryptography applications (e.g., Elliptic Curves). Rapid feedback on circuit characteristics is used to evaluate deep architectural changes in short time, greatly reducing the timeto-market with respect to hand-made designs. In addition, our approach is highly flexible, since the same generic high-level description can be used to produce an entire set of circuits, each with different area/performance trade-offs. Thanks to the proposed approach, any change to the initial specification (e.g., the curve used) is also very fast, while it may require a large effort in the case of hand-made designs.

I . I N T R O D U C T I O N

Current trends are pushing hardware description towards higher abstraction. In this context, High-Level Synthesis (HLS) tools allow for a concise description of circuits [1]. They offer the best trade-off between precision, the asset of lower level description languages such as VHDL or Verilog, and the rapid prototyping often associated with higher level description languages such as C or C++. These higher-level languages are historically used in the software development field, even if the differences between hardware and software design are numerous. Hardware development is a task involving more tools and steps than software development in the sense that the path toward a functional prototype is longer. Moreover, architectural changes often involve a consequent rewriting and testing workload. On the other hand, incremental and iterative methods (also called agile methods) are often applied in the software field to cope with these issues [START_REF] Larman | Iterative and incremental development: A brief history[END_REF]. These methods are of multiple types, but they have a few rules in common including the objective of having a working prototype as fast as possible. The following design iterations aim at improving the existing version without breaking anything. Such methods are not easily applicable to hardware development due to the usage of verbose hardware description languages.

This paper presents the results of an Incremental and Iterative Design (IID) method applied to the design of a generic Elliptic Curve Cryptography (ECC) crypto-coprocessor, based on an HLS generation flow. HLS is commonly used for some mainstream applications, but the specific case of ECC requires architectures based on uncommon operations (e.g., large integer multiplications or reductions modulo a large prime) that are less suitable for current HLS tools optimizations. In this particular case, it is not straightforward to consider a fully automatic HLS flow giving results that can be compared with the state of the art. This case study is focusing on this example to highlight how an HLS flow can still be used and how an IID method can compensate for some of HLS drawbacks.

By combining HLS and IID during development, we were able to iterate quickly in order to test multiple architectural changes until obtaining a satisfactory solution. Additionally, the generated circuit is described in a high-level language. Therefore, it is very versatile as it supports updating its parameters before synthesis. Choosing the supported elliptic curve or how large integers are partitioned remains easy. Finally, software tests during the development of prototypes give quick feedbacks before hardware tests validate the generated circuits.

The paper is organized as follows. We will first describe the motivations of this work. Then, we will present the proposed design flow and how the IID method is used. The results are then discussed in Section IV.

I I . M O T I VAT I O N S

In the first place, the objective was to obtain quickly an efficient but versatile coprocessor for Elliptic Curve Cryptography. This case study was chosen because it is a difficult case for HLS tools, which are usually efficient for classical integer operations. Indeed, very large integers are generally used in cryptography. We will first provide some details about this application to emphasize development choices before summarizing the bases of our methodology.

A. Designing and Updating an ECC Coprocessor

Using Elliptic curves for asymmetric cryptography was initially proposed in 1985 by Miller [START_REF] Miller | Use of elliptic curves in cryptography[END_REF] and Koblitz [START_REF] Koblitz | Elliptic curve cryptosystems[END_REF]. Most elliptic curves used in cryptography are described by quadratic equations on two (or more) variables, where each term is defined over a binary or a prime Galois field. To use such cryptographic mechanism, participants should agree on one particular elliptic curve which defines a related security level. Due to recent progress in mathematical attacks against curves over binary fields [START_REF] Koblitz | A riddle wrapped in an enigma[END_REF], elliptic curves over prime field seem to be a safer (Miller) [START_REF] Miller | Use of elliptic curves in cryptography[END_REF] to NIST [START_REF] Vanstone | Responses to NISTs proposal[END_REF] ISO 14888-3 [START_REF]Information technology -security techniques -digital signatures with appendix -part 3: Discrete logarithm based mechanisms[END_REF] (NIST) [START_REF] Fips | -2. digital signature standard (DSS)[END_REF] SEC-2 [START_REF] Certicom | Sec 2: Recommended elliptic curve domain parameters[END_REF] standard [START_REF] Brainpool | Ecc brainpool standard curves and curve generation[END_REF] Chinese SM2

French RGS curve448 [START_REF] Hamburg | Ed448-goldilocks, a new elliptic curve[END_REF] curve25519 [START_REF] Bernstein | Curve25519: new diffie-hellman speed records[END_REF] IEEE 1363a [START_REF] Assoc | IEEE std 1363a-2004[END_REF] ECDSA in ECDSA in

IEEE 1363a

IEEE 1363 (Koblitz) [START_REF] Koblitz | Elliptic curve cryptosystems[END_REF] ECC ECC starts [START_REF] Assoc | IEEE std 1363-2000, IEEE standard specifications for public-key cryptography[END_REF] starts [START_REF] Assoc | IEEE std 1363a-2004[END_REF] ANSI X9.62 [START_REF] American | Ansi x9.62[END_REF] IEEE 1363 [START_REF] Assoc | IEEE std 1363-2000, IEEE standard specifications for public-key cryptography[END_REF] ECDSA in ECDSA in FIPS186-2

Certicom publishs

German Brainpool frp256v1 [START_REF] Anssi | Mécanismes cryptographiques: Règles et recommandations concernant le choix et le dimensionnement des mécanismes cryptographiques[END_REF] S. Vanstone proposes ECDSA [START_REF]China cryptography administration: SM2 EC recommended parameters[END_REF] 1985 Fig. 1: ECC history.

choice. Therefore, we focus on curves over GF (p) (prime field) in this paper. Figure 1 illustrates the ECC history by underlining some major events. Real life utilizations date back to the early 21st century after governmental and international standardizations, but also recommendations from consortiums [START_REF]Information technology -security techniques -digital signatures with appendix -part 3: Discrete logarithm based mechanisms[END_REF]- [START_REF] Assoc | IEEE std 1363a-2004[END_REF]. This was followed in October 2005 by a new standard with an alternative curve generation method [START_REF] Brainpool | Ecc brainpool standard curves and curve generation[END_REF]. Recently, new curves with specific parameters were proposed [START_REF] Bernstein | Curve25519: new diffie-hellman speed records[END_REF], [START_REF] Hamburg | Ed448-goldilocks, a new elliptic curve[END_REF] and are already used in some applications. It is even possible for anyone to send random data in order to participate in the seed for a new curve generation [START_REF] Lenstra | A random zoo: sloth, unicorn, and trx[END_REF]. Other elliptic curves exist but are rarely used or with less extensive support by tools. Figure 1 illustrates that new curves are regularly proposed. In the end, there is still no unanimous consensus about which curve should be used.

In this context, focusing a hardware development on a specific elliptic curve could lead to the design of an outdated circuit at short term. Our proposition is therefore to use reconfigurable devices (FPGA) to prevent this. The proposed building process is compatible with all elliptic curves over prime field even if each built bitstream is specialized for a specific curve. This permits to efficiently minimize the resource usage, but requires to provide a simple method to update the supported curve. This is possible thanks to an automatic building process as explained in section II-C. The specialization for a selected curve is done during this build. This automatic process should also permit to adapt the ratio performance/area. For instance, from the same input description, we are able to produce very cheap (but slow) crypto processors, or a fast but more resource consuming one.

ECC can be used for authentication, cyphering, or secret exchange protocols. All are based on the same time consuming operation: the scalar multiplication, which is a one-way function as no known sub-exponential algorithm currently exists to reverse it. Since this is the bottleneck of our application, the case study in this paper is focused on this operation. The

Require: P = a point of the elliptic curve O = the infinite point k = a n-bit width scalar = (k n-1 k n-2 ...k 0 ) 2 Ensure: Q = k × P 1: Q ← O 2: for i = n -1 to 0 do 3: R ← 2 × Q point doubling 4:
S ← R + P point addition always done 5:

if k i = 1 then 6: Q ← S Write-back: Q = 2 × Q + P 7: else 8: Q ← R Write-back: Q = 2 × Q 9:
end if Fig. 2: The Double-And-Add Always ECSM algorithm: point doublings in white and point additions in black or grey elliptic curve scalar multiplication (ECSM) is not a classical integer multiplication, as it is defined between a secret scalar (classically a 256-bit integer) and a point of the elliptic curve, whose coordinates are described by large integer values. Several algorithms for the ECSM exist, but they are usually computed by using the same two basic point operations: point addition and point doubling. These operations are not the classical doubling or addition on integer, because their operands are actually points of an elliptic curve. Although we will not consider here the threat of an attacker with physical access to the device, remote attacks to the device are often realistic. If an attacker can communicate with the device, he will be able to measure the computation time. This can be an important information for accelerating cryptanalysis [START_REF] Kocher | Timing attacks on implementations of diffie-hellman, RSA, DSS, and other systems[END_REF], and this is why only constant time algorithms are usually chosen. For the ECSM, the simplest constant time algorithm is the Double-And-Add Always that is illustrated in Figure 2. In order to compute the ECSM kP , the scalar coefficient k is scanned from the most significant to the less significant bit. For each bit, the point accumulator Q is doubled and stored in R. A point addition P + R is always computed and stored in S because of constant time requirement. The write-back in the accumulator Q depends on the bit k i value (black or grey boxes in Figure 2). For a 256-bit scalar, the algorithm will always require 255 point doublings and 255 point additions.

Specific formulas are used to compute coordinates of the addition and doubling resulting points. Performance of these formulas depends on the chosen point representation. For performance considerations, we have selected the Jacobian representation [START_REF] Cohen | Handbook of elliptic and hyperelliptic curve cryptography[END_REF]. The formulas require large integer operations over the prime field GF (p): addition, subtraction, inversion, and multiplication, all modulo p. The prime number p can be different for each elliptic curve. Some curves over prime fields may have specific p values in order to implement more efficient modular reductions [START_REF] Fips | -2. digital signature standard (DSS)[END_REF], but there are also curves with generic primes and thus less efficient reductions. Since we want to support any curve over arbitrary prime fields, the latter approach is thus chosen. In order to minimize the impact of this choice on global performance, the values are represented in the Montgomery domain, which allows us using a Montgomery multiplier and thus avoiding the time consuming reduction operation [START_REF] Montgomery | Modular multiplication without trial division[END_REF]. This multiplier is sequential as the result is computed by splitting the operands into smaller digits which are processed one cycle at a time. The elementary multiplication between two digits is computed in 1 clock cycle and defines the main critical path. The smaller the digit, the higher the operating frequency, but the more cycles the full integer multiplication will require. Evaluating the optimal digit width is therefore a difficult problem.

B. Fast Design Method

Cryptographic choices impose some constraints but due to many freedom degrees, there are always a lot of non-imposed implementation details. These details can have important effects on the overall performance which can be difficult to anticipate. We thus need a method allowing a fast manual exploration of the different solutions with relevant and precise indicators. This means that performance or cost (area) indicators should be quickly available to fix design choices as soon as possible. That is why we should have a functional system in a short time, even if sub-optimal, which we will be able to improve later by incremental and iterative development steps. The core of our method is hence to import and adapt IID from the software field in order to try and rewrite as many circuit descriptions as required to reach the target. This adaptation may not seem straightforward for hardware development but it will eventually allow a fast design process.

C. Using HLS for IID

To adapt an IID method to hardware development, the key idea is to use an HLS flow in order to benefit from its advantages. First, a designer needs precise feedback on the resources and the performance of his circuit. He also needs this information as soon as possible during the design process, in order to react and rewrite or adapt the source code. HLS allows quickly producing a functional circuit as the source code can be compiled and tested like any software. On the other hand, the same source code can also go through an HLS flow, where the tool can provide useful information to the developer, such as the area of the circuit and performance estimation. With such data, the developer can immediately change the source code and both test functionality and estimate performance. One should note that these iterations are not the classical HLS user choices (e.g., loop unrolling) but source code modifications in order to choose for instance a different algorithm.

Secondly, HLS allows for writing complex circuits in a more concise way compared to a classical method with Hardware Description Languages (HDLs). A minor variation in the source code given to an HLS tool can lead to significant architectural changes that would take a large amount of manpower in order to apply them in HDL.

Eventually, higher-level languages allow for a very powerful set of tools (the C preprocessor for instance) giving the possibility to produce versatile circuits. As seen earlier, in the ECC context and because of the rapid evolution of the curve set, such a feature is mandatory. In our case, as an example, the digit width was never fixed during the development. It was kept as a parameter that permits to generate a set of cryptoprocessors for a curve with different area/performance ratios. In the end, the designer can select the circuit satisfying his constraints among the set of generated crypto-processors.

For all these reasons, our proposal is focused on describing a hardware circuit in a high-level language (not a HDL) to feed an HLS tool. Figure 3 gives an overview of the different steps that the source code goes through during the design flow. It is exploited in two independent build chains: the software one, used to quickly test the functionality, and the hardware one, to check the hardware functionality (slower than in software) and produce a FPGA-ready bitstream. The tedious debug of hardware functionality is (at least partially) replaced by a faster software debug.

D. Summary of Needed Operations and Constraints

The top operation is the ECSM. There are hierarchical dependencies between needed operations:

• The scalar multiplication: a special multiplication between a scalar and a point of the elliptic curve. The scalar length is close to the length of the prime p. • Basic point operations: point doubling and point addition.

• Operations over GF (p): multiplication, subtraction, addition and inversion of large integers modulo a large prime number. The typical size of the prime and of the operands is currently 256-bit for a reasonable security level, but it may increase in the future.

The sub-section IV-A presents our design method on a limited part of the crypto-processor: the multiplication over the prime field GF (p), the most important operation used in ECSM. On the other hand, we compare our final circuits to the state of the art at the top-level (the ECSM) in Section IV-B. In addition, three constraints must be respected: • All operations must be computed in constant time to protect the circuit against timing attacks. • The elliptic curve must be a parameter (chosen before synthesis). This permits an easy update when the supported curve of a design is outdated or a new one is chosen. This also permits to generate a set of circuits supporting different curves.

• The digit width must also be a parameter (chosen before synthesis) because the optimal size is difficult to anticipate before the place and route step. Additionally, this permits to generate a set of circuits with different area/performances ratios.

I I I . H L S A N D I I D

In this section the proposed design flow and the IID method are discussed in more details. The tools used for building the ECC coprocessor are then presented.

A. Design Flow

In this paragraph, we describe the method used during hardware design. Figure 4 represents the flowchart of the design phases. The first step is to describe the hardware design using the C language (actually, a subset of ANSI C cf. III-B). This version can be compiled with gcc for instance in order to be executed on a standard general purpose processor. The outputs of the executed program are checked to assert the algorithmic correctness of the program. As long as the program outputs are wrong, there are corrections to include in the C source code. When the modifications are done, the C source code is given to the HLS tool which produces a hardware description of the program in standard HDL. The HLS tool can give some performance information about the circuit. If these performances are not satisfactory, the designer can try improving the C source code. The output of the HLS tool is then given to an HDL simulator for verification and another corrective iteration can occur. Eventually, the designer can use any tool chain, according to the vendor of the targeted FPGA, to generate the final circuit and obtain precise information that can lead to further performance-oriented iterations.

B. Tools

Keeping the developed application up to date with current curves is a major concern. In order to take advantage of newer FPGAs as they are released, the tools used to build the presented circuit are relatively independent of the target. Moreover, this tool chain is Free and Open Source Softwares (FOSS) for the most part. However, the need for a final estimation of the design choices and the FPGA implementation imply the usage of vendor tools so that specific FPGAs can be targeted, and the design evaluated and built for this specific chip. In this paper, a Xilinx chip is targeted, hence, the Xilinx tool suite is used further in this work. Note that any tool suite from other FPGA vendors could have been used instead.

Figure 5 references all the tools used in the studied design flow. This entire flow presented in III-A is encapsulated in an automated environment handled by make. The first step is to modify the C input to adapt the code to the chosen parameters (particularly the curve type). This code is given to gcc for compilation, the produced binary is tested, and the correct outputs are kept as golden references. In parallel, the HLS tool Augh [START_REF] Prost-Boucle | Fast and standalone design space exploration for high-level synthesis under resource constraints[END_REF] generates the circuit in VHDL and reports estimation of the final circuit latency and performance (i.e. frequency). It also generates a TCL script in order to automatize the final circuit generation through Xilinx's Vivado.

Augh is used in this paper for the following reasons:

• it produces a VHDL output very quickly • it is vendor agnostic and produces code for various FPGA families • it performs a design space exploration automatically, the user is not involved in the hardware circuit optimization process • Eventually, Augh is an academic free and open source software, hence results are reproducible. Augh takes a subset of ANSI-C as input and produces VHDL files.

I V. E VA L U AT I O N O F T H E D E S I G N M E T H O D

How the results were obtained is here as much important as what the results are. Hence, performance of every obtained circuit are analysed with respect to the development timeline.

(c) (e) (d) (f) (g) (h) (i) (j) (k) (k) (l) (l) (n) (m) (o)
Fig. 6: Area and latency of the Montgomery multiplier vs. time for 256-bit operands and 32-bit digits.

The mere results are also analysed with a non-exhaustive state of the art comparison, in order to validate the obtained architecture.

A. Architecture Performances

In this paragraph, we will discuss Figure 6 which presents the performance and the cost of the obtained Montgomery multiplier implementations relatively to a three-weeks timeline. Version controlled code allowed us to plot such data. It should be noted that with the IID method, every obtained circuit is a fully functional ECC coprocessor. However, we focus on the Montgomery multiplier in order to provide the simplest illustration of our development method. All field multipliers support arbitrary digit width, but for ease of clarity only those based on 32-bit digits are represented.

During the first days of development, a functional software application was developed. Automatic test based on the classical cryptography library openssl were also implemented to prevent non functional circuits. The first functional hardware circuit was obtained after this period: (a) in Figures 6 and7. This circuit consumed few resources because parallelism is not identifiable in the first written code. For the same reason, this implementation was inefficient. Quickly (b), useless carry propagations were identified and suppressed: computation time was reduced by 10 times. The multiplier algorithm contains several accumulation and bit shifting operations. The next improvement (c), illustrated in Figure 8, was to merge the bit shifting with the previous operation by modifying the writeback destination. Right immediately, the C code was relaxed to permit parallelism between partial products and a Montgomery specific digit multiplication before shifting.

After some time, a latency increase can be observed (d). The final step of our multiplication is a conditional subtraction, which is always computed to avoid timing attacks. The writeback is thus disabled when this operation is not required. So far, the digits of the multiplier operands are incrementally scanned: in order to reduce routing, operands were hence copied into local rotating buffers (e). However, it appears that the area cost of these rotating buffers is very large. Merging two close loops

(c) (e) (d) (f) (g) (h) (i) (j) (k) (l) (n) (m) (o)
Fig. 7: Area × latency of the Montgomery multiplier vs. time for 256-bit operands and 32-bit digits.

Require: a = (a m-1 , ..., a 2 , a 1 , 0) Ensure: a = a >> digit width = (0, a m-2 , ..., a 1 , a 0 )

1: a 0 = 0 2: for i = 1 to m -1 do Compute a 3:

a i = ... 4: end for 5: for i = 0 to m -2 do Shift a 6:
a i = a i+1 7: end for if i > 0 then 3:

a i-1 = ...

4:

end if 5: end for 6: a m-1 = 0 Fig. 8: Description of the iteration (c). offered a latency reduction by increasing the parallelism, but it also increases the occupied area (f). After two days, it was clear that our rotating buffers do not fit well on FPGA targets, hence they were abandoned. Anticipating this fact was not obvious. Dedicated memory components exist in our FPGA target (Xilinx Virtex 7 family): RAM blocks were therefore introduced in the design (g). This significantly reduced the area occupied by the Montgomery multiplier. In addition, these memories are large enough to store several integers. Operating frequency also improved slightly.

The Montgomery multiplier requires a register for the accumulator, which was hence moved into the BRAMs of its first operand thus simplifying the multiplier usage (h). This reduced the occupied area but was unfavourable for performances. Loop unrolling was useful for easiest parallelism of multiplications without data dependency, but it had a significant cost on area and it indirectly slew down the system frequency due too additional for j = 0 to m -1 do Body without data dep. ... = a j × r b + ... Read a j

5:

end for 6: end for

⇓ 1: for k = -2 to m 2 -1 do
Body without data dep.

2:

if k >= 0 then compute i and j indexes for the next cycle 7: end for Fig. 9: Description of the iteration (k).

routing. After having studied the generated unrolled scheduling, rewriting an internal loop of our multiplier for executing only operations without data dependency was easy (i). This code describes a hardware pipeline, hence loop unrolling for this particular portion of the code was forbidden (j). The pipelined loop is in a second loop: each loop is used to scan one of the operands. For each digit of the second operand, the final writeback of the pipeline was stalled. Figure 9 illustrates the loop merging (k). Less pipeline stalls reduce the pipeline latency, significantly reducing the global latency of our multiplier.

Finally, the Montgomery multiplier requires a conditional final subtraction, which has been deleted by relaxing the constraint on the range of integers (l). RAM blocks are inferred by Xilinx's Vivado from the VHDL generated by Augh whenever a register linked to a memory is found (i.e. an array in C). However, the synthesis tool sometimes chose to absorb this register into an input register of a DSP block instead of using it as output register for a RAM block. After solving this problem, some additional data was moved from LUTRAM to BRAM (m). The latency was reduced thanks to a shorter critical path, and some area was saved by using a more compact digit index coding (n). After a final compression of pipeline stages, the final circuit was built (o).

A lot of radical hardware architecture changes happened during the development: moving data into BRAMs; pipelining the Montgomery multiplier; algorithmic update; deletion of the final subtraction in the multiplier. There has also been changes in several parts of the design out of the multiplier. However, all these modifications were always related to small iterations in the C source code. A lot of small improvements that were done by iterative development on our high-level description would have been very time consuming in classical hardware design.

B. Related Work

In this paragraph, the results of the final architecture for computing the scalar multiplication (ECSM) are compared with state-of-the-art solutions. This architecture contains the final Montgomery multiplier obtained after the iterations previously detailed. It also contains other operations over GF (p), point addition, point doubling, and the top-level scalar multiplication. All these components have been quickly designed using the same methodology.

ECSM computation time is selected for performance comparison. To evaluate the complete circuit cost, two metrics are proposed. The first one (named unified metric) takes into account BRAMs, DSPs, and slices. The number of BRAMs and DSPs is converted in equivalent-slices to allow a fair comparison. Conversion ratios are computed for each FPGA device as the ratio between the number of available slices and the number of available BRAMs and DSPs. The unified metric as seen in Figure 11 is the maximum of occupied equivalent-slices by BRAMs, DSPs, and really occupied slices. This unified metric and the convertion ratio can be calculated with the following formulas:

r D = #slices target #DSPs target r B = #slices target #BRAMs target unif. metric = max(#slices, r D .#DSPs, r B .#BRAMs)
In the case of our target (xc7v585t) the ratios become:

r D 585t = 91050 1260 72.3 r B 585t = 91050 1292 70.5
If we consider the last circuit obtained after the iterations and synthesized for a 256-bit elliptic curve with a digit width of 32 bits (cf. Table I), the unified metric can be calculated as follows:

unif. metric = max(633, r D 585t × 8, r B 585t × 1) = 633

Another metric taking into account DSPs and BRAMs is proposed in [START_REF] Vansteenkiste | Analyzing the divide between FPGA academic and commercial results[END_REF]. Unfortunately, it uses FPGA-specific architectural information that is tedious to gather for the different targets presented here. The second circuit cost metric is classic: it is the number of occupied slices, where DSPs and BRAMs are not taken into account. The former metric is more adapted to the case where all system parts use BRAMs, DSPs and slices. In Figure 11, area is represented by the latter, which is for cases where BRAMs and DSPs are not used by other parts of the system. Thus, they can be considered as carrying no cost and only slices are included in the area cost. In order to compare area fairly, the study is limited to Xilinx Virtex families 5, 6 and 7 because in all these targets each slice has 4 6-input LUTs. For equitable performance comparisons, only implementations supporting 256-bit elliptic curves over the prime field are represented.

In order to obtain interesting area/performance ratios, multiplications over GF (p) must be sequential. [START_REF] Roy | Tile before multiplication: An efficient strategy to optimize DSP multiplier for accelerating prime field ECC for NIST curves[END_REF] exhibited the advantages of using DSP in this context on a Virtex 5 target. The number of occupied slices can be reduced from 9100 to us digit width=32 digit width=54 digit width=18 (61 circuits) our circuit set Fig. 10: ECSM Implementation efficiency vs. area for stateof-the-art circuits and our circuit set using final Montgomery multiplier implementation (*: unified metric is defined in IV-B) [START_REF] Donda | Efficient implementation of a generic coprocessor for elliptic curve cryptography on reconfigurable hardware[END_REF] [25] [START_REF] Bajard | Double level montgomery cox-rower architecture, new bounds[END_REF] [27] [START_REF] Roy | Tile before multiplication: An efficient strategy to optimize DSP multiplier for accelerating prime field ECC for NIST curves[END_REF] [24] [28] [START_REF] Ma | A high-speed elliptic curve cryptographic processor for generic curves over GF(p)[END_REF] us Fig. 11: ECSM Implementation efficiency vs. area (slices) for state-of-the-art circuits and our circuit set using final Montgomery multiplier implementation 4505 by using 16 DSPs. This number of DSPs is equal to 4320 equivalent-slices (16×r D = 16× 51840 192 for a target xc5vlx330). The usage of these specific components is reasonable because it is close to the number of occupied slices. In the same way, computation time is reduced to 0.59 ms for each scalar multiplication. This work is based on a manual tiling of DSPs by taking into account the asymmetry of these components. This tiling step should be repeated if the implemented circuit has to support a larger prime. In addition, the modular reduction is specific to a special prime number [START_REF] Fips | -2. digital signature standard (DSS)[END_REF] and is not compatible with random curves.

To avoid the reconfiguration required to support multiples curves, it is possible to design crypto-processors supporting several curves. For example, [START_REF] Alrimeih | Fast and flexible hardware support for ECC over multiple standard prime fields[END_REF] proposes an implementation supporting 5 NIST curves with different prime widths. It is based on modular operations over a restricted set of prime fields GF (p). The proposed circuit is efficient (scalar multiplication in 0.4 ms for a 256-bit curve on a Virtex 6 target). On the other hand, a lot of resources are required because the design also supports a 521-bit curve: 11200 slices, 289 DSPs, and 128 BRAMs. In addition, the number of DSPs is large comparing to occupied slices. Our proposed area metric in Figure 10 is 39657 equivalent-slices whereas only 11200 slices are occupied. For an optimal resource occupation, this design could be used in conjunction with an application that does not consume DSPs. In order to minimize area, our choice is different because the generated designs support only one curve but an automatic and fast process permits easier curve update.

Another method for increasing the frequency is to use RNS arithmetic [START_REF] Bajard | Double level montgomery cox-rower architecture, new bounds[END_REF]. This can be more complex than classical arithmetic but more efficient due to shorter critical path (0.612 ms for a 256-bit scalar multiplication on a Kintex 7). In order to support a new curve with another prime, however, a new RNS base must be generated.

All our generated designs use less resources than the other circuits. Our most efficient circuit occupies a little less slices than designs proposed by [START_REF] Donda | Efficient implementation of a generic coprocessor for elliptic curve cryptography on reconfigurable hardware[END_REF] on Kintex 7. In this paper, a GF (p) Montgomery Multiplier is used as well. The development was done at the RTL level whereas we used HLS and IID methods. For comparable performances, our circuits consume less area due to the larger exploration that is done by the HLS tool. Althought these circuits use as much slices as our best circuits, the unified metric (in Figure 10) illustrates their intense DSP usage. It should be pointed out that performance and area overheads in their circuits come in part from a protection against physical attacks.

It is possible to increase the Montgomery multiplier throughput. For example, a previous implementation on Virtex 5 target used quotient pipelining to achieve more parallelism [START_REF] Ma | A high-speed elliptic curve cryptographic processor for generic curves over GF(p)[END_REF]. This permits to compute a scalar multiplication in 0.376 ms with 1725 slices and 37 DSPs. If we want to generate more efficient circuits by increasing occupied resources, another iteration in our development may integrate this improvement on the multiplier pipeline. A simpler improvement may also be developed to increase parallelism, i.e. compute at the same time several multiplications. Data dependencies in basic point operations allow this but in order to minimize area this was left as future work.

Another very efficient circuit was proposed by Sasdrich and Günneysu on Zynq 7 target [START_REF] Sasdrich | Efficient elliptic-curve cryptography using curve25519 on reconfigurable devices[END_REF]. Data rolls between two memories through the arithmetic operator. The multiplier is based on a school method with an efficient reduction because the field is based on a pseudo Mersenne prime. This implementation can support only one specific curve (Montgomery curve). For this particular case, the number of prime field multiplications required for computing one point doubling and one point addition is reduced from 26 to 10. It seems difficult to achieve the same performance when targeting generic curves. A future iteration in the development focusing on specific curves and keeping our actual solution as fall-back for generic curves is always possible, but out of the scope of this paper.

C. Generated Set

Our generated circuits support an arbitrary generic elliptic curve contrary to some cited designs [START_REF] Alrimeih | Fast and flexible hardware support for ECC over multiple standard prime fields[END_REF], [START_REF] Sasdrich | Efficient elliptic-curve cryptography using curve25519 on reconfigurable devices[END_REF]. We focused on 256-bit curves for comparison to the state-of-the art, TABLE I: Area and performance of the scalar multiplication over curve P-256 [START_REF] Certicom | Sec 2: Recommended elliptic curve domain parameters[END_REF] for an extract of our generated circuits but other security levels can be chosen. In addition, several crypto-processors with different area/performance ratios can be generated by modifying the digit width used in the large integer partitioning. The target is a Virtex 7 FPGA (XC7V585T) and all results are post place and route. Our circuits in Figures 10 and 11 are generated with digit widths from 4 to 64 on curve P-256 [START_REF] Certicom | Sec 2: Recommended elliptic curve domain parameters[END_REF]. With the same high-level source code, our generated circuits occupy from 417 to 1019 slices (respectively with digit width 4 and 64) and can compute a scalar multiplication on 256-bit from 97.8 ms to 2.4 ms (respectively with digit width 4 and 54). A subset of all the generated circuits is characterized in Table I. The design with best performance uses 54-bit digits and computes the scalar-multiplication in 2.4 ms (651000 cycles at 271 MHz). The 64-bit digits is faster in term of cycles but the critical is longer hence computation time is stretched. Our designs consume less area than other designs thanks to the proposed method, while achieved performances are close to equivalent hand-made designs.

D. Discussion

The present contribution is focusing on FPGA for various reasons. A highly versatile architecture was required given the fast changing trends in elliptic curve standards. Nonetheless, it is possible to apply the presented method to any hardware design that lacks HLS support. Indeed, when a satisfying architecture is found after few iterations, the hardware developer can directly write this kind of architecture in HDL and proceed to the low level improvements that differentiates the output of an HLS tool from a fully hand-made design. With this solution, the developer can avoid potentially costly architecture changes and have final results quicker than with a classical method.

V. C O N C L U S I O N

The reported case study demonstrates the efficiency of the proposed approach, inspired by some principles of Agile development methods. Thanks to HLS, powerful and flexible circuits can be obtained quickly, with various performance/area trade-offs. The proposed flow also enabled the generation of a highly versatile Elliptic Curve Cryptography coprocessor, adapted to any existing and future curve specification.

Future works include the parallelization of some architecture parts. This improvement could allow us to reach state-ofthe-art latencies. Our verification based on a cryptographic library is application specific. Better interoperability by using standardized verification methodology could be beneficial. Moreover, proving equivalence of designs between iterations could further improve verification. Taking into account side channel vulnerabilities (in addition to constant computation time) during the high-level specification is another objective.

R E F E R E N C E S

Fig. 3 :

 3 Fig. 3: Generic design flow and tool chain requirements.

Fig. 4 :

 4 Fig. 4: Flowchart of HLS design methodology.

Fig. 5 :

 5 Fig. 5: Set of tools used to setup design flow.

⇓ 1 :

 1 for i = 0 to m -1 do Compute and shift a 2:

  Require: a = (a m-1 , ..., a 0 ) b = (b m-1 , ..., b 0 ) Ensure: a × b >> m × digit width mod p 1: for i = 0 to m -1 do 2: r b = b i Read b i 3:

3 :

 3 ... = a j × r b + ... r b update: r b = b i Read b i 6:


Jérémie Schutz

Optimal "sporadic" and systematic preventive maintenance policy for leased equipment under various operating conditions

Keywords: Leased equipment, Working conditions, Maintenance policy, Finite horizon

In this paper, the considered system corresponds to a specific equipment proposed for leasing. This equipment is leased to several lessees who use it under various working conditions. The aims of this paper consist in determining optimal maintenance plans which minimize maintenance costs. Two maintenance policies are compared : systematic policy (imperfect preventive maintenance actions are performed, with the same effectiveness factor, after each mission) and "sporadic" policy (imperfect preventive maintenance action can be performed after a mission with its optimal effectiveness factor).

Introduction

In this research area, the first work on the themes of maintenance is usually attributed to Barlow and Hunter [START_REF] Barlow | Optimum Preventive maintenance policies[END_REF]. In this work, these authors were interested in determining an optimal time for replacement equipment thanks to the development of two major strategies known under the name "Age-Based Maintenance" and "Block-Based Maintenance". Subsequently, many authors have determined, for these periods of renewal, maintenance policies in order to minimize maintenance costs. Also, we can refer to the work of Nakagawa [START_REF] Nakagawa | Periodic and Sequential Preventive Maintenance Policies[END_REF] who has made a very significant contribution in this area with periodical and sequential policies. Periodic policies are characterized by preventive maintenance actions performed at constant time intervals. Sequential policies consist in determining the optimal interval between each preventive action. Generally, for these policies, preventive maintenance actions are considered perfect (the equipment is restored to a state "As Good As New"). Corrective actions are considered minimal (the failure rate of equipment remains unchanged "As Bad As Old"). In 1988, Kijima et al. [START_REF] Kijima | Periodical replacement problem without assuming minimal repair[END_REF] proposed a first model for imperfect maintenance. It is characterized by an age reduction of the system.This reduction corresponds to a proportional amount of the time elapsed since the previous preventive maintenance activity. Nakagawa [START_REF] Nakagawa | Sequential Imperfect Preventive Maintenance Policies[END_REF] proposed an alternative approach to modeling imperfect maintenance. It is characterized by an intensity increase in the failure rate after every preventive maintenance. However, after each imperfect action, the failure rate is reduced to zero ; the system can be considered new. All these models assume that the working of the studied system is constant over time. Alas, in reality, an equipment can operate under different operational conditions (e.g. production rate) or environmental conditions (e.g. temperature). Therefore, the degradation of the system depends on these conditions and the equipment can thus degrade faster or more slowly. In the literature, two models are used to represent the degradation variations. Accelerated Life Model (ALM) and Proportional Hazards Model (PHM) influence the reliability function (or hazard function) by adding a risk function [START_REF] Martorell | Age-dependent reliability model considering effects of maintenance and working conditions[END_REF] [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF]. The ALM changes the age of the system while the PHM defines a variation of the failure rate proportional to the working conditions. Thanks to over half a century of research work and considering many important economic crises in recent years, recent work has sought to develop and implement maintenance policies for leased equipment. In recent work, the most notable contributions are those of Jaturonnatee et al. [START_REF] Jaturonnatee | Optimal preventive maintenance of leased equipment with corrective minimal repairs[END_REF], Pongpech and Murthy [START_REF] Pongpech | Optimal periodic preventive maintenance policy for leased equipment[END_REF] and Yeh et al. ( [START_REF] Yeh | Optimal threshold value of failure-rate for leased products with preventive maintenance actions[END_REF], [START_REF] Yeh | Optimal preventive maintenance policy for leased equipment using failure rate reduction[END_REF], [START_REF] Yeh | A study of maintenance policies for second-hand products[END_REF], [START_REF] Yeh | Preventive-maintenance policy for leased products under various maintenance costs[END_REF]). In these works, we find the usual policies such as sequential and periodic ones with imperfect maintenance based on the age reduction or failure rate increase. In 2007, Yeh and Chang offered an innovative policy where preventive maintenance actions are realized when the failure rate reaches a preset threshold.

Based on these observations, the aims of this paper consist in determining optimal maintenance plans which minimize maintenance costs for a leased equipment. This equipment is leased to several lessees with various working conditions. The remainder of this paper is organized as follows. Section 2 gives a studied system description, notations and working assumptions. In section 3, mathematical formulations of the maintenance policies and resolution methods are proposed. A numerical example is presented in the section 4. Finally, section 5 gives the conclusion and the future work.

Problem description

In this paper, the considered system corresponds to a specific equipment proposed for leasing. Lessees are placed in a queue and obtain equipment for the desired period (duration of the mission) as it becomes available. Between two missions, this equipment can undergo imperfect preventive maintenance actions (based on age reduction) performed by the lessor (owner). These actions do not have negligible durations unlike the minimal corrective maintenance activities. The duration to perform preventive actions is a time percentage α of the length of all missions. The minimal corrective maintenance actions are performed during the missions by the lessor, when the system fails. The effect of working conditions, for a mission m, is based on the proportional hazard model (PHM) where they are modeled by a risk function g m [START_REF] Cox | Regression Models and Life-Tables[END_REF]. So, the hazard function is given by g m • λ (t).

Notations

Throughout the paper, the following notations will be used: 3 Mathematical formulation

Minimal maintenance policy

The purpose of this policy consists only in maintaining the system in a state of working to assess the contribution of preventive maintenance policies. This maintenance policy consists solely of minimal corrective maintenance actions. These activities are carried out to overcome the immobilizing failures that arise over time. Maintenance total cost is given by:

Γ (M) = C CM • φ (M) (1)
Usually, when the system operates under constant working conditions during a period [0, T ], the average number of failures corresponds to the cumulative hazard function at the time T . However, in this research, the missions have various operating conditions and the failure rate evolves with the latter. So, the reliability of the system depends also on the working conditions but with minimal maintenance activities, the reliability must be continuous over time. In Schutz et al. [START_REF] Schutz | An integrated strategy for efficient business plan and maintenance plan for systems with a dynamic failure distribution[END_REF], a functional age was defined to ensure a continuous reliability. To illustrate this functional age, let's consider two consecutive missions i and j. The reliability at the end of the mission i (at time ν i +δ i ) under condition z i must be equal to the reliability of the functional age (denoted ν j ) under condition z i . The functional age is given by:

ν m = R -1 [R (ν m-1 + δ m-1 )] zm z m-1 (2)
where R (•) and R -1 (•) respectively denote the reliability function and its inverse. Therefore, from the equation (1), the average number of failures can be expressed by:

φ (M) = dim(M) m=1 z m • νm+δm νm λ (t) dt (3)

Systematic preventive maintenance policy

The first improved maintenance policy consists in planning preventive maintenance activities systematically after the completion of each mission. More precisely, the same preventive maintenance actions must be achieved after the (dim (M) -1) first missions. Indeed, after the last mission, a perfect preventive maintenance action is carried out. For this policy, total maintenance costs are specified by the following equation:

Γ (M, ρ) = C CM • φ (M, ρ) + C P M (ρ) • (dim (M) -1) (4)
where the average number of failures is given by:

φ (M, ρ) = dim(M) m=1 z m • νm,ρ+δm νm,ρ λ (t) dt (5)
From the equation ( 5), the average number of failures may seem identical to the equation (3) (minimal maintenance policy). However, the difference is in functional age expression. As preventive maintenance (system age reduction according to the effectiveness factor ρ) are carried out after the missions, functional age is expressed by:

ν m,ρ = (1 -ρ) • R -1 [R (ν m-1,ρ + δ m-1 )] zm z m-1 if m > 1 0 else (6)
Determination of the optimal effectiveness factor The effectiveness factor "plays" a role in the quality of preventive action implementation. Due to this factor, preventive maintenance is considered imperfect and the system is returned to a state between ABAO (As Bad As Old) and AGAN (As Good As New). Therefore, it seems logical that this factor may also be involved in modeling the costs and lengths of preventive maintenance. For example, the competence of the operators can influence the quality of maintenance activities, but the maintenance action cost also depends on these skills. Although costs may depend on the effectiveness factor, they are not necessarily proportional to the latter. Consequently, the cost of preventive action is based on a fixed one (e.g. parts and products used) and a variable one (qualification level, experience, etc.). Similarly, the length of preventive maintenance is composed of a fixed and a variable duration based on ρ. In this research, the cost and duration are modeled by:

C P M (ρ) = C P M,F + C P M,V (ρ)
(7)

µ P M (ρ) = µ P M,F + µ P M,V (ρ)
where C P M,F and C P M,V correspond respectively to the part of fixed and variable costs. Similarly, the fixed and variable durations are µ P M,F and µ P M,V .

As mentioned in section 2, the total time allocated to preventive actions is a percentage (noted α in the equation ( 9) of cumulative duration of all missions. For this systematic maintenance policy, the effectiveness factor is the same for all preventive actions. Thus, the maximum value of rho is obtained from equation ( 8) and its interval is given by:

ρ ∈          0, min     1, log - C P M,F •(dim(M)-1)-α• dim(M) m=1 δm (dim(M)-1) log (C P M,V )     maximum value of ρ          (9)
Integral limits of the function φ sys (M, ρ) are too complex to analytically determine the ρ factor. The latter will be determined by a numerical resolution.

This maintenance policy, compared to the minimal maintenance policy, can be further improved by determining the optimal value of the effectiveness factor for each preventive action.

"Sporadic" preventive maintenance policy

This preventive maintenance policy is called "sporadic" because after each mission, preventive actions can be (or not) performed with diverse effectiveness factors. Basically, this model expression is very similar to the previous model. The total maintenance cost is given by:

Γ (M, P ) = C CM • φ (M, P ) + dim(M)-1 m=1 C P M (ρ m ) • ρ m ( 10
)
with P the effectiveness factors vector of the first (dim (M) -1) missions and ρ m represents the ceil of the factor ρ m . So, if the effectiveness factor is equal to 0, no preventive maintenance is performed (one of type ABAO has no interest).

The average number of failures, for this policy, is expressed by:

φ (M, P ) = dim(M) m=1 z m • νm,ρ m +δm νm,ρ m λ (t) dt (11)
and the expression of the functional age remains unchanged (cf. equation ( 6))

This model resolution is more difficult. Indeed, there is not one decision variable but... dim (P ) decision variables. Here again, the determination of these effectiveness factors cannot be computed analytically. Given the duration allocated to preventive maintenance actions, effectiveness factors must satisfy the following relationship:

dim(P ) m=1 µ P M (ρ m ) ≤ α • dim(M) m=1 δ m (12)
Unlike the systematic maintenance policy, the ρ factor can take any value in the interval ]0, 1[. The next section presents a numerical example where the objective consists in determining the effectiveness factor of preventive maintenance in the case of systematic and sporadic policies.

Numerical example

Let's consider the following arbitrarily chosen input data to illustrate our model:

-λ (t) = 2.5 300 • t 300 (2.5-1) (Weibull distribution, shape = 2.5, scale = 300)

-C CM = 500 mu (money unit) -

C P M = 300 + (600) ρ mu -µ P M = 2 + (8)
ρ tu (time unit) -α = 0.05 Various missions are presented in Table 1. When maintenance actions performed are only of minimal type, the average number of failures is about 6.66 and the total cost of maintenance activities is 3331.5 mu. With this cost, it will be possible to estimate profit from the two preventive policies described in the previous section.

Systematic maintenance policy

As the total duration of all missions is 628 tu, the duration for preventive maintenance actions is 31,40 tu (α • 628). After each mission, preventive activities can last, at the maximum, 6.28 tu ( α•628 6-1 ). During this period allocated to (imperfect) preventive maintenance, the maximum value of effectiveness factor is 0.69. Based on maintenance costs given above, the total cost of maintenance activities amounted to 2289,7 mu with ρ = 0.58 and an average number of failures about 1.17. This first improved policy generates a profit of 31 %. When the effectiveness factor is less than 0.19, this maintenance policy is less interesting than the minimal policy. This is due to the slight decrease in the number of failures compared to preventive actions cost.

"Sporadic" maintenance policy

As this preventive maintenance policy is more flexible, results are even better. From the numerical resolution, the best choice is to perform preventive maintenance activities after the second mission (with ρ = 0.74) and after the fourth mission (with ρ = 0.78). With these optimal decision variables, total maintenance cost is reduced by 46 % compared to the minimal maintenance policy. The average number of failures is 1.87, the cumulative duration of preventive maintenance actions is 13.72. In this case, extra time (31, 4 -13.72) can be used to add preventive maintenance actions to reduce the number of failures. However, total cost will be higher...

Conclusion and prospects

In this paper, we considered a leased equipment for a defined horizon. During the lease period, the system was used in various working conditions. The latter affected the system, which degraded differently. To maintain a continuous reliability, functional age was introduced to generate equivalences between different operational conditions and working times. Thanks to this functional age, three maintenance policies have been studied. The first maintenance, minimal, is used to assess preventive maintenance policies benefit. Systematic preventive maintenance policy, characterized by the same preventive activities after each mission, give the best economical result. However, this result depends on the ratio between preventive and corrective costs (in the section referring to the numerical example, if ρ < 0.19, this policy should be avoided). The " sporadic" policy allows preventive maintenance action after the chosen missions. It is therefore more flexible but the determination of the decision variables requires much more time. It would be interesting to use a meta-heuristic for solving this problem when the number of missions increases.

Among possible perspectives, it would be interesting to study the case where the maximum duration after each preventive mission depends directly on mission duration. These maintenance policies may also be more realistic considering imperfect corrective actions and no negligible duration for these activities. The lessee can impose criteria such as equipment availability. In this case, preventive actions could be performed during missions. Another perspective would be to consider a set of different possible missions (with time windows). The aim would be to select missions to perform and to determine an optimal maintenance plan.

-

  M: Missions Vector to perform -δ m : Duration of the mission m -λ (t): Hazard function for nominal conditions -Γ (•): Total cost of maintenance policy -Φ (•): Average number of failures -ρ: Effectiveness factor of preventive maintenance action -C CM : Corrective maintenance action cost -C P M (ρ): Preventive maintenance action cost based on effectiveness factor ρ Other notations used in this document are defined below equations.

Fig. 1 .

 1 Fig. 1. Evolution of the total maintenance cost based on the effectiveness factor

Table 1 .

 1 Missions to perform during the lease period

	Mission 1	2	3	4	5	6
	δm	106 107 82 104 108 121
	gm	1.24 1.48 1.41 1.42 1.28 1.12


Yuqi Mo
email: mo.yuqi@insa-lyon.fr

Minh-Tien Do
email: minhtien.do@sigfox.com

Claire Goursaud
email: claire.goursaud@insa-lyon.fr

Jean-Marie Gorce
email: jean-marie.gorce@insa-lyon.fr

Optimization of the predefined number of replications in a Ultra Narrow Band based IoT network

Keywords: Wireless sensor network, UNB transmission, Random-FTMA schemes, replication mechanism

Ultra Narrow Band networks using Random-FTMA scheme are mainly characterized by the randomness both in time and frequency domain. The frequency randomness prevents from allocating orthogonal resources for transmission, and induces uncontrolled interference. By introducing diversity, replication mechanism is a promising candidate to enhance the reliability of such wireless network. Therefore, in this paper, by taking into account the randomness both in time and frequency domain, the replication mechanism has been considered. Considering the outage probability, we theoretically evaluate the system performance and show that there exists an optimal number of transmissions. Finally, we highlight that this number of repetitions can be easily optimized by considering a unique global parameter.

I. INTRODUCTION

IoT (Internet of things) networks are gaining more and more interest as they will provide connectivity for countless number of devices [START_REF] Ishaq | IETF standardization in the field of the internet of things (IoT): a survey[END_REF]- [START_REF] Zanella | Internet of Things for smart cities[END_REF]. When the goal is to collect information (for example, for temperature monitoring, parking availability monitoring, smart metering, ...), the devices will require to transmit only a few amount of data, during a small proportion of time. In this case, the base station should cover as much devices as possible to reduce the cost of infrastructure deployment, and not waste transmission resources. This is the reason why long range technologies such as LoRa [5] and Ultra Narrow Band (UNB) [6] have emerged recently. The main advantage of UNB over LoRa is that it does not rely on access protocol to handle multiple access [START_REF] Goursaud | Dedicated networks for IoT: PHY / MAC state of the art and challenges[END_REF]. Thus, the transmission resources needed for channel reservation are saved.

Nevertheless, UNB suffers from interference that may lead to the loss of some transmitted packets. Indeed, in SigFox UNB based network, devices access to the medium at will. Transmissions are seen as random, both in time and frequency (R-FTMA : Random Frequency and Time Multiple Access), due to the non-negligible jitter impact [START_REF] Do | On the benefits of Random-FDMA schemes in ultra narrow band networks[END_REF]. This induces some simultaneous transmission in the same frequency area, leading to interference. The system performance of such network was evaluated in [START_REF] Do | Interference modelling and analysis of random fdma schemes in ultra narrowband networks[END_REF], and the capacity in term of the devices number was estimated. However, no control is made on the transmission, so when the interference is too high the packet information is definitively lost.

In the literature [START_REF] Mahmood | Reliability in wireless sensor networks: A survey and challenges ahead[END_REF], the reliability of data delivery can be achieved by two ways: retransmission-based mechanism and redundancy-based mechanism. The basic principle of the retransmission schemes is to repeatedly transmit the failed message which cannot be recovered at the receiver. The demand of replication happens until a correct message arrives. To know that the message is correctly received or not between the transmitter and receiver, an acknowledgment mechanism must be applied. However, this kind of protocol requires a downlink transmission, so the reservation of transmission resources. This might further aggravate the inherent congestion of the medium.

Besides, in the redundancy-based mechanism, multiple copies of the same message are transmitted based on Erasures Coding that allows a receiver to recover the information from independent message loss. The popular types of Erasure Coding are Reed-Solomon codes, lowdensity parity-check (LDPC) codes, Fountain Codes, and the simplest one is the replication mechanism.

The performance of retransmission and redundancybased mechanism is theoretically evaluated in [START_REF] Wen | replication or redundancy: Transmission reliability in wireless sensor networks[END_REF] in term of the energy efficiency. Indeed, both retransmission-based mechanism and redundancy-based mechanism are able to improve the data delivery reliability, but also consume a lot of battery-energy of a sensor node which is inherently limited. Besides, this theoretical analysis highlights that Erasure Coding mechanism is out-performing compared to the replication mechanism in term of reliability and energy efficiency, for the case of a low packet loss rate and low number of hops.

In the UNB network using Random-FTMA schemes, the one-hop communication between active node and basestation is essential to reduce the energy consumption. We thus consider the replication redundancy-based mechanism in this paper. We consider only the transmission reliability in terms of the outage probability according to the number of replications. To this aim, the coding, channel and the efficient consumption of battery energy is neglected in our study for the sake of simplicity. Based on the PHY analysis (SINR: Signal to Interference plus Noise Ratio), we define that the message failure occurs when the BER calculated at the base-station is less than the BER threshold (typically 10 -3 or 2 • 10 -4 for BPSK modulation). Each node transmits its message to basestation with a predefined number of occurrences regardless of the energy consumption, and the success state of the previous transmissions. One may expect that the repetition increases the chances of receiving at least a copy correctly. But, it also increases the load, thus, the probability of collisions, leading to a compromise. Such phenomenon has already been studied for classical transmissions [START_REF] Liu | Novel packet Retransmission in ofdma systems using frequency diversity[END_REF]- [START_REF] Sudhaakar | A Novel QoS-Aware MAC Scheme Using Optimal Retransmission for Wireless Networks[END_REF], but, to the best of the authors' knowledge, not for UNB transmissions where the frequency multiple access is specific. In this paper, we thus aim at finding the optimal number of replications for such network.

The rest of this paper is organized as follows. The assumptions and modeling are detailed in Section II. Next, we theoretically evaluate the system performance for the simple replication scheme in Section III. Then, we exploit the obtained formulas to estimate the optimum replication number in Section IV. Finally, Section V gives the conclusion.

II. MODELING AND ASSUMPTION

A. Network Topology

In this study, the network is reduced to a unique basestation (BS) covering a large number of nodes positioned in this coverage. We consider that all the nodes are received at the base station with the same power. This corresponds to the case with a power control loop.

In this study, we consider a point in time where N nodes are active. At a given time, only a few portion of these N nodes are actually transmitting. The sleep or wake-up duty cycle may vary among the different nodes. Each node randomly chooses its transmission time slot.

More specifically to UNB, the choice of carrier frequency is performed independently in each node without any prior knowledge on the occupation of the channel and can be discrete or continuous [START_REF] Do | On the benefits of Random-FDMA schemes in ultra narrow band networks[END_REF]. At the base station, the whole transmission band BW is analyzed. When a transmission is detected, the received signal is filtered at the corresponding frequency before demodulation. In case of simultaneous transmission, the interference depends on the spacing δ f = |f 1 -f 2 | between the desired node frequency f 1 and the interferer's one f 2 :

P I (t) = P (δ f , t) (1)
The interference level for a realistic UNB filter shape is plotted on Fig. 1. We can observe that the interference level is low if the frequency spacing δ f is high. Thus, a unique interferer leads to failure if δ f is small.

The system performance is evaluated through the outage probability (OP). A data transmission is considered successful if the received SINR is above a predefined threshold γ * , with :

SINR = P s (W + P I ) ( 2
)
where W is the filtered additive white Gaussian noise power assumed to be negligible compared to the interference level, which is a reasonable value for UNB systems. We consider γ * = 8 dB (corresponding to the Bit-Error-Rate 2 • 10 -4 for our numerical results as we consider the BPSK modulation. We have determined that, with the

B. Replication Mechanism

Any active node i ∈ [1, • • • , N ] randomly selects a carrier f i and a time slot t i . Failure occurs when 2 active nodes (a desired node x and interfering node y) transmit simultaneously, and the frequency shift between their carriers δ f = |f x -f y | falls in the interval [0; 123] Hz.

To improve the service probability for each transmitting node, the replication mechanism illustrated in Figure 2 is applied. We assume that each message (with transmission duration d) should be transmitted in less than the message lifetime T b . As the system is considered slotted, there are T b d slots available for transmission during the message lifetime. We consider that the message is transmitted exactly n r times during the message lifetime T b , regardless of the success of previous transmissions. n r is thus considered predefined and identical for all nodes. The length of the window allocated for a replication is thus T b nr•d . Each node randomly selects a time slot in this window to repeat its message. The message is successfully transmitted when at least 1 of the n r attempts succeeds. Otherwise, the message is considered lost.

III. THEORETICAL ANALYSIS A. Outage Probability derivation

In this section, we derive a theoretical expression for the outage probability, as a function of the main system parameters : the whole transmission band BW ; the frequency occupancy with respect to the carrier [-b; b] Hz; the wake-up duty cycle of nodes T b ; the time duration of a message d, and the replication number n r . The first four factors influence the channel occupancy of each replication in time and frequency domains.

The goal of the base station is to decode all transmitted messages, from all the users. However, without sake of generality, in this analysis, we consider the transmission status of a unique user, called "desired user". We evaluate its probability to have its message lost, due to other users called "interferers".

We first consider the case where the message is transmitted with a single replication. A node wakes up every T b seconds (one may note that the time unit can also be minute, hour, ...) to send a message who lasts d seconds. Thus, for time-slotted multiple access, with slot duration d, the probability that a given potential interferer transmits during the same slot as the desired user can be expressed by :

P t = d T b (3)
The same principle applies for the frequency domain. In this theoretical analysis, to simplify the expression, we suppose that transmission errors are due to a unique interferer. We neglect the cases where 2 or more users jointly generate an interference level that leads to an error. Therefore, with this assumption, there is an error if an interferer elected a carrier frequency spaced by less than b compared to the desired user's one. Thus, the frequency interfering zone that causes a failure if the frequency gap of non-desired node and the desired node has a width 2b. We can deduce that the frequency collision probability is:

P f = 2b BW (4)
A collision occurs when two nodes send their messages at the same time, and when their frequency spacing is too small. As P f and P t are independent, the collision

probability is P t • P f = 2b × d BW × T b .
As a consequence, when considering N independent potential users, the probability for a given user message being correctly received, without considering replication mechanism, is the probability that none of the N -1 other user creates collision:

P (C = 0 | n r = 1) = 1 - 2b × d BW × T b N -1 (5)
with C the number of collisions, and n r the number of replications.

We now extend this expression to the case with replication mechanism. First, we modify (5) to take into account the fact that the number of transmissions increases. Indeed, the collision probability of a replication with any replication of another user is multiplied by n r and becomes

2b × d BW × T b × n r .
Thus, the probability that a given replication is correctly received can be expressed as:

P (C = 0 | n r ) = 1 - 2b × d BW × T b × n r N -1 (6)
As replications are identical, the message is successfully received if at least one of this replication is correctly received at the base station. We make the assumption that the replications experience independent transmission conditions, and that the success probability of any replication is independent of the collision on the previous ones. Accordingly, the outage probability of one useful message, transmitted with n r replication, is:

OP = 1 -1 - 2b × d BW × T b × n r N -1 nr (7)
One may note that, in practice, the collision probability of a given replication depends on the previous replications state. Indeed, a given replication of an interferer can collide with at most one replication of the desired user. Thus, the collision probability should be updated for each replication, depending on the number of collisions observed during the previous replications. However, as the number of users (thus the number of transmissions) is very high, the collision probability decrease can be neglected.

B. Validation

In this section, we validate by simulation the theoretical OP established in eq.7. We consider the impact of 4 parameters : BW , T b , n r and N . We do not consider b variation as this parameters is constrained by the transmission technology. Furthermore, as the system is slotted, the exact values of d and T b are not of importance, but rather the number of available slots for each message: T b d . Thus, we fix d = 1, and only vary T b .

We have plotted on Figures 345, the comparison of the theoretical OP with the simulated one, as a function of the replication numbers, when varying N , BW , or T b . We can first observe that the curves does not match perfectly, but theory still provides a good estimation of the observed OP, especially for low n r . More importantly, the theory correctly describes the evolution of the OP as a function of n r .

Besides, these figures highlight the fact that the replication mechanism is able to greatly improve the probability of successful transmission of a message. Indeed, the OP can be decreased compared to the one replication case. This is due to the fact that the increased number of replications n r initially permits to increase the chance that at least one packet is successfully received. However, the repetition also multiplies the number of messages sent by the nodes within a fixed time interval, and thus increases the channel utilization and the collision probability. Therefore, for a large number of replication n r , the temporal and frequency resources become overloaded, and the benefit of adding redundancy is disabled by the increase of interference level. Thus the replication process is both increasing the collision probability (thus the message error probability), while also increasing the chance to have at least one replication successfully received (thus decreasing the message error probability). Consequently, we can observe on Figure 3-5, that there is an optimum number of replications, that permits to obtain the highest message throughput.

As the theory evolution is consistent with the simulated one, we use from now on eq.7 to determine the optimum number of replications.

IV. REPLICATION NUMBER OPTIMIZATION

In this section, we focus on the number of replications that permits to get the lowest OP. A. Preliminary estimation of n ropt OP is influenced by many parameters. For sake of simplification, we call λ the time-frequency collision factor, with:

λ = 2b × d BW × T b (8)
Hence OP simplifies to:

OP = 1 -(1 -λ × n r ) N -1 nr (9)
We derive eq.9 with respect to n r to find out n ropt .

(

) 10
dOP dn r = (N -1)λn r (1 -λn r ) N -2 1 -(1 -λn r ) N -1 + ln(1 -1 -λn r ) N -1 × 1 -(1 -λn r ) N -1 nr
As OP is convex, n ropt is obtained when dOP dn r = 0. We do not consider the case where

(1 -(1 -λn r ) N -1 ) nr
is equal to 0, as this would lead to a null OP. So the expression of n ropt can be deduced from the first term of the multiplication. However, n ropt can not be solved analytically. Nonetheless, its value can be numerically evaluated, depending on the values of λ and N .

We have plotted on Figure 6 the evolution of n ropt when considering different time-frequency factor λ and number of active users N . This figure shows that when the time-frequency factor λ increases, the optimal number of replication n ropt decreases. This means that when collision probability is already important, the replication mechanism make transmission even more difficult. On the contrary, for low λ, the replication mechanism permits to take advantage of the under-used channel to improve the transmission quality. Similarly, for a given λ, the number of replications should be decreased as the number of active users increases.

Thus, from this theoretical analysis, we can confirm that the optimal number of replication decreases when the network is overloaded, which could be caused by : too many active users at the same time; users sending message too often; or a small available bandwidth.

B. Parameters simplification

We can further our analysis by evaluating OP when considering a constant repartition of users in time or frequency. To do so, we have computed eq.7 for special cases of N BW = cste and N T b = cste, and plotted the results on Figure 7 and Figure 8.

As shown in Figure 7 and Figure 8, we can observe that when the parameters BW and T b evolve proportionally to N , the behavior of the OP according to the number of replications n r is unchanged. Therefore, it is possible to simplify our study by reducing the number of parameters from 3 (i.e., BW , T b and N ) to 2 (i.e., N BW = cste and N T b = cste). But, by jointly exploiting the above properties, we can further simplify to a single parameter : From the previously tested sets, we can extract the value of the optimal number of replications n ropt for different node density Ω. It may be noted that obtained couples (n ropt , Ω) correspond to optimal values of different OP. The evolution of the optimal number of replications n ropt according to user density Ω is plotted in Figure 9. We can observe that, n ropt decreases when the node density gets high (this is due to the fact that the increase of the network load reduces the available time/frequency resources for the replication process). This result is very interesting and promising because its urges the potential ability of configuring the network by relying on a single global parameter: node density over the temporal and frequency resources Ω.

V. CONCLUSION

In this paper, we have considered a UNB based network using Random-FTMA scheme. We have proposed and studied a repetition mechanism to improve the reliability of the transmissions. The outage probability is theoreti-cally derived and evaluated under different scenarios with parameters such as the number of replications n r , the number of nodes N , the bandwidth length BW and the message lifetime T b . We have shown that the outage probability can be evaluated by considering, in addition to n r a single criteria : the node density over the temporal and frequency resources Ω. Besides, we have highlighted that, for each node density Ω, there is an optimal number of replications n ropt , which we have estimated. With such results, the number of replications can be easily predefined and adapted to the network load. Finally, one may note that this study can be furthered by estimating the minimal number of replications needed to obtain a targeted OP.

Figure 1 .

 1 Figure 1. Interference behavior P I as a function of the frequency spacing δ f between desired user and interferer.

Figure 2 .

 2 Figure 2. Illustration of replication nr = 2 for duration of message d, period of message T b

Figure 3 .

 3 Figure 3. OP vs number of replications nr with different number of active users N , BW = 12 kHz, T b = 75s.

Figure 4 .

 4 Figure 4. OP vs number of replications nr with different bandwidth BW , T b = 75s, N = 1000.

Figure 5 .

 5 Figure 5. OP vs number of replications nr with different message period T b , BW = 12 kHz, N = 1000.

Figure 6 .Figure 7 .

 67 Figure 6. Optimal replication nr opt vs time-frequency factor λ and active user numbers N

Figure 8 .

 8 Figure 8. OP vs number of replications nr with BW = 192 kHz, N T b = 100.

Figure 9 .

 9 Figure 9. Node density Ω vs optimal number of replications nr opt .

ACKNOWLEDGMENT

The authors would like to thank Christophe Fourtet (Sigfox, Toulouse) for his helpful inputs regarding the UNB parameters during this study.