Florent De Dinechin
email: florent.de.dinechin@ens-lyon.fr

Bogdan Pasca

Octavian Cret
email: octavian.cret@cs.utcluj.ro

Radu Tudoran

AN FPGA-SPECIFIC APPROACH TO FLOATING-POINT ACCUMULATION AND SUM-OF-PRODUCTS LIP RESEARCH REPORT RR2008-22

This article studies two common situations where the flexibility of FPGAs allows one to design application-specific floating-point operators which are more efficient and more accurate than those offered by processors and GPUs. First, for applications involving the addition of a large number of floating-point values, an ad-hoc accumulator is proposed. By tailoring its parameters to the numerical requirements of the application, it can be made arbitrarily accurate, at an area cost comparable for most applications to that of a standard floating-point adder, and at a higher frequency. The second example is the sum-of-product operation, which is the building block of matrix computations. A novel architecture is proposed that feeds the previous accumulator out of a floating-point multiplier without its rounding logic, again improving both area and accuracy. These architectures are implemented within the FloPoCo generator, freely available under the GPL.

INTRODUCTION

Most general-purpose processors have included floatingpoint (FP) units since the late 80s, following the IEEE-754 standard. The feasibility of FP on FPGA was studied long before it became a practical possibility [START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machine[END_REF][START_REF] Li | Implementation of single precision floating point square root on FPGAs[END_REF][START_REF] Ligon | A re-evaluation of the practicality of floating-point operations on FPGAs[END_REF]. As soon as the sizes of FPGAs made it possible, many libraries of floating-point operators were published (see [START_REF] Belanović | A library of parameterized floating-point modules and their use[END_REF][START_REF] Lee | Parameterisable floating-point operators on FPGAs[END_REF][START_REF] Lienhart | Using floatingpoint arithmetic on FPGAs to accelerate scientific Nbody simulations[END_REF][START_REF] Roesler | Novel optimizations for hardware floating-point units in a modern FPGA architecture[END_REF] among other). FPGAs could soon provide more FP computing power than a processor in single precision [START_REF] Lienhart | Using floatingpoint arithmetic on FPGAs to accelerate scientific Nbody simulations[END_REF][START_REF] Roesler | Novel optimizations for hardware floating-point units in a modern FPGA architecture[END_REF], then in double-precision [START_REF] Zhuo | Scalable and modular algorithms for floating-point matrix multiplication on FP-GAs[END_REF][START_REF] Dou | 64-bit floating-point FPGA matrix multiplication[END_REF][START_REF] Delorimier | Floating-point sparse matrix-vector multiply for FPGAs[END_REF]. Here single precision (SP) is the standard 32-bit format consisting of a sign bit, 8 bits of exponent and 23 bits of significand (or mantissa), while double-precision (DP) is the standard 64-bit format with 11 bits of exponent and 52 significand bits. Since then, FPGAs have increasingly been used to accelerate scientific, financial and other FP-based computations. This acceleration is essentially due to massive parallelism: Basic FP operators in an FPGA are typically slower than their processor counterparts by one order of magnitude.

Most of the aforementionned applications are very close, from the arithmetic point of view, to their software implementations. They use the same basic operators, although the internal architecture of the operators may be highly optimised for FPGAs [START_REF] Liu | An iterative division algorithm for FPGAs[END_REF]. Most published FP libraries are fully parameterisable in significand length and exponent length, but applications that exploit this flexibility are rare [START_REF] Roesler | Novel optimizations for hardware floating-point units in a modern FPGA architecture[END_REF][START_REF] Zhuo | Scalable and modular algorithms for floating-point matrix multiplication on FP-GAs[END_REF].

The FloPoCo project 1 studies how the flexibility of the FPGA target can be better exploited in the floating-point realm. In particular, it looks for operators which are radically different from those present in microprocessors. In the present article, such operators are presented for the ubiquitous operation of floating-point accumulation, and applied to sums of products.

All the results in this article are obtained for Virtex4, speedgrade -12, using ISE9.1, and should be reproducible using Xilinx WebPack and FloPoCo, both available at no cost. Very similar results have been obtained for Altera Stratix II. FloPoCo produces portable VHDL, and is designed with the goal of automatic fine-tuning the architectural parameters to the target hardware and frequency.

FLOATING-POINT ACCUMULATION

Summing many independent terms is a very common operation. Scalar product, matrix-vector and matrix-matrix products are defined as sums of products. Numerical integration usually consists in adding many elementary contributions. Monte-Carlo simulations also involve sums of many independent terms. Many other applications involve accumulations of floating-point numbers, and some related work will be surveyed in section 6.

If the number of summands is small and constant, one may build trees of adders, but to accomodate the general case, it is necessary to design an iterative accumulator, illustrated by Figure 1.

It is a common situation that the error due to the computation of one summand is independent of the other summands and of the sum, while the error due to the summation grows with the number of terms to sum. This happens in integration and sum of products, for instance. In this case, it makes sense to have more accuracy in the accumulation than in the summands.

A first idea is to use a standard FP adder, possibly with a larger significand than the summands. The problem is that FP adders have long latencies: typically l = 3 cycles in a processor, up to tens of cycles in an FPGA (see Table 1). This is explained by the complexity of their architecture, illustrated on Figure 2.

This long latency means that an accumulator based on an FP adder will either add one number every l cycle, or compute l independent sub-sums which then have to be added together somehow. This will add to the complexity and cost of the application, unless at least l accumulation can be interleaved, which is the case of large matrix operations [START_REF] Zhuo | High performance linear algebra operations on reconfigurable systems[END_REF][START_REF] Bodnar | Floating-point accumulation circuit for matrix applications[END_REF].

In addition, an accumulator built out of a floating-point adder is inefficient, because the significand of the accumulator has to be shifted, sometimes twice (first to align both operands and then to normalise the result, see Figure 2). These shifts are in the critical path of the loop of Figure 1.

In this paper, we suggest to build an accumulator of floating-point numbers which is tailored to the numerics of each application in order to ensure that 1/ its significand never needs to be shifted, 2/ it never overflows and 3/ it eventually provides a result that is as accurate as the application requires. We also show that it can be clocked to any frequency that the FPGA supports. We show that, for many application, the determination of operator parameters ensuring the required accuracy is easy, and that the area can be much smaller for a better overall accuracy. Finally, we combine the proposed accumulator with a modified, errorless FP multiplier to obtain an accurate application-specific dot-product operator.

A FAST AND ACCURATE ACCUMULATOR

This section presents the architecture of the proposed accumulator. Section 3.2 will discuss the determination of its many parameters in an application-specific way.

Overall architecture

The proposed accumulator architecture, depicted on Figure 3, removes all the shifts from the critical path of the loop by keeping the current sum as a large fixed-point number. Figure 4 illustrates the accumulation of several floatingpoint numbers (represented by their significands shifted by their exponent) into such an accumulator.

There is still a loop, but it is now a fixed-point addition for which current FPGAs are highly efficient. Specifically, the loop involves only the most local routing, and the dedicated carry logic of current FPGAs provides good performance up to 64-bits. For instance, a Virtex4 with speed grade -12 runs such an accumulator at more than 220MHz, while consuming only 64 CLBs. Section 3.3 will show how to reach even larger frequencies and/or accumulator sizes.

The shifters now only concern the summand (see Figure 3), and, being combinatorial, can be pipelined as deep as required by the target frequency.

As seen on Figure 3, the accumulator stores a two's complement number while the summands use a sign/magnitude representation, and thus need to be converted to two's complement. This can be performed without carry propagation: if the input is negative, it is first complemented (fully in parallel), then a 1 is added as carry in to the accumulator. All this is out of the loop's critical path, too.

1 1 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1 0 1 0 1 0 0 0 0 0 0 1 1 1 1 1 w A = MSB A -LSB A + 1 Accumulator w F + 1 LSB A = -12 MaxMSB X = 8 MSB A = 16

Parameterisation of the accumulator

Let us now introduce, with the help of Figure 4, the parameters of this architecture.

• w E and w F are the exponent and significand size of the summands

• The weight MSB A of the most-significant bit (MSB) of the accumulator has to be larger than the maximal expected result, which will ensure that no overflow ever occurs.

• The weight LSB A of its least-significant bit can be arbitrarily low, and will determine the final accuracy as Section 4 will show.

• For simplicity we note w A = MSB A -LSB A the width of the accumulator.

• MaxMSB X is the maximum expected weight of the MSB of a summand. MaxMSB X may be equal to MSB A , but very often one is able to tell that each summand is much smaller in magnitude than the final sum. In this case, providing MaxMSB X < MSB A will save hardware in the input shifter.

The main claim of the present article is the following: For most applications accelerated using an FPGA, values of MaxMSB X , MSB A and LSB A can be determined a priori, using a rough error analysis or software profiling, that will lead to an accumulator smaller and more accurate than the one based on an FP adder. This claim will be justified in section 4.

This claim sums up the essence of the advantage of FP-GAs over the fixed FP units available in processors, GPUs or dedicated floating-point accelerators: we advocate an accumulator specifically tailored for the application to be accelerated, something that would not be possible or economical in an FPU.

Fast accumulator design using partial carry-save

If the dedicated carry logic of the FPGA is not enough to reach the target frequency, a partial carry-save representation allows to reach any arbitrary frequency supported by the FPGA. As illustrated by Figure 5, the idea is to cut the large carry propagation into smaller chunks of k bits (k = 4 on the figure), simply by inserting (MSB A -LSB A )/k registers. The critical path is now that of a k-bit addition, and the value of k can therefore be chosen to match the target frequency. This is a classical technique which was in particular Of course a drawback of the partial carry-save accumulator is that it holds its value in a non-standard redundant format. To convert to standard notation, there are two options. One is to dedicate (MSB A -LSB A )/k cycles at the end of the accumulation to add enough zeroes into the accumulator to allow for carry propagation to terminate. This comes at no hardware cost. The other option, if the running value of the accumulator is needed, is to perform this carry propagation in a pipelined way before the normalisation -this is the dashed box on Figure 3. The important fact is again that this carry propagation is outside of the critical loop.

Post-normalisation unit, or not

Figure 3 also shows the FloPoCo LongAcc2FP component, that performs the conversion of the long accumulator result to floating-point. However, let us first remark that this component is much less useful than the accumulator itself. Indeed, some applications need the running sum at each cycle, but most applications need only the final sum and have no need for all these intermediate sums. Let us take a few examples.

In [START_REF] Cret | FPGA-based acceleration of the computations involved in transcranial magnetic stimulation[END_REF], the FPGA computes a very large integrationseveral hours-and only the final result is relevant. In such Fig. 5. Accumulator with 4-bit partial carry-save. The boxes are full adders, bold dashes are 1-bit registers, and the dots show the critical path. cases, of course, it makes no sense to dedicate hardware to the conversion of the accumulator back to floating-point. FPGA resources will be better exploited at speeding up the computation as much as possible. FloPoCo provides a small helper program to perform this conversion in software.

Another common case is that a single post-normalisation unit may be shared by several accumulators. For instance, a dot product of vectors of size N accumulates N numbers before needing to convert the result back to floatingpoint. Therefore, in matrix operations, one pipelined LongAcc2FP may be shared between N dot product operators [START_REF] Zhuo | High performance linear algebra operations on reconfigurable systems[END_REF].

For these reasons, it makes sense to provide LongAcc2FP as a separate component, as on Figure 3. Note that the same holds for an accumulator based on an FP adder of latency l (that actually computes l intermediate subsums). If only the final sum is needed, it may be computed in software at no extra hardware cost. However, if the running sum is needed at each cycle, it will take l -1 additional FP adders (and some registers) to add together the l sub-sums. And finally, such a post-normalisation unit can be shared as in the long accumulator case for matrix operations [START_REF] Zhuo | High performance linear algebra operations on reconfigurable systems[END_REF][START_REF] Bodnar | Floating-point accumulation circuit for matrix applications[END_REF].

To sidestep this discussion, in the following, we choose to compare the performance only for the accumulators themselves. We nevertheless provide for completeness some synthesis results for the post-normalisation unit in Table 2, but the reader should have in mind that comparable extra costs will plague an accumulator based on an FP adder.

Back to LongAcc2FP, it mostly consists in leadingzero counting and shifting, followed by conversion from 2's complement to sign/magnitude, and rounding. If the accumulator holds a partial carry-save value, the carries need to be propagated -this simply requires w A /k pipeline levels, each consisting of one k-bit adder and w A /k -1 registers of k bits. Again, all this may be performed at each cycle and pipelined arbitrarily. For an IEEE-754-like correctlyrounded conversion, one must in addition compute a sticky bit out of the discarded bits, but it is yet unclear if users are ready to pay this cost.

Performance results

Table 1 illustrates the performance of the proposed accumulator compared to one built using a floating-point adder from the Xilinx CoreGen tool.

For each summand size, we build accumulators of twice the size of the input significand (MSB A = w E , LSB A = -w E ) for two configurations: a small one where MaxMSB X = 1, and a larger one where MaxMSB X = MSB A = w E . As section 4 will show, a typical accumulator will be between these two configurations. Again, these results are for illustration only: An accumulator should be built in an application-specific way.

Table 2 provides preliminary results for the current implementation of the LongAcc2FP post-normalisation unit.

APPLICATION-SPECIFIC TUNING OF ACCUMULATOR PARAMETERS

Let us now justify the claim, made in 3.2, that the many parameters of our accumulator are easy to determine on a per-application basis.

A tale of performance, and also accuracy

First note that a designer has to provide a value for MSB A and MaxMSB X , but these values do not have to be accurate. For instance, adding 10 bits of safety margin to MSB A has no impact on the latency and very little impact on area. Now from the application point of view, 10 bits mean 3 orders of magnitude: it is huge. A designer in charge of implementing a given computation on FPGA is expected to understand it well enough to bound the expected result with a margin of 3 orders of magnitude. An actual example is detailed below in 4.2. As another example, in Monte-Carlo simulations of financial markets, each accumulation computes an estimate of the value of a share. No share will go beyond, say, $100,000 before something happens that makes the simulation invalid anyway.

It may be more difficult to evaluate MaxMSB X . In doubt, MaxMSB X = MSB A will do, but the accumulator will be larger than needed, and application knowledge should help reduce it. For instance, implementers of Monte-Carlo simulation will exploit the fact that probabilities are [START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machine[END_REF] 142 slices, 6 cycles @ 386 MHz (8,23) SP 304 slices, 7 cycles @ 366 MHz [START_REF] Li | Implementation of single precision floating point square root on FPGAs[END_REF]37) 546 slices, 10 cycles @ 325 MHz (11,52) DP 863 slices, 12 cycles @ 306 MHz Table 2. Synthesis results for a LongAcc2FP compatible with Table 1, rounding an accumulator of size 2w F to an FP number of size w F . smaller than 1. Another option is profiling: a typical instance of the problem may be run in software, instrumented to output the max and min of the absolute values of summands. Again, the trust in such an approach comes from the possibility of adding 20 bits of margin for safety.

(w E , w F ) LongAcc2FP, 2w F → w F (7,
In some cases, the application will dictate MaxMSB X but not MSB A . In this case, one has to consider the number n of terms to add. Again, one will usually be able to provide an upper bound, be it the extreme case of 1 year running at 500MHz, or 2 53 cycles. In a worst-case scenario on such simulation times, this suggests the relationship MSB A = MaxMSB X + 53 to avoid overflows. For comparison, 53 is the precision of a DP number, so the cost of this worst case scenario is simply a doubling of the accumulator itself, but not of the input shifter which shifts up to MaxMSB X only. It will cost just slightly more than 53 LUTs in the accumulator (although much more in the post-normalisation unit if one is needed).

The last parameter, LSB A , allows a designer to manage the tradeoff between precision and performance. First, remark that if a summand has its LSB higher than LSB A (case of the 5 topmost summands on Figure 4), it is added exactly, entailing no rounding error. Therefore, the proposed accumulator will be exact (entailing no roundoff error whatsoever) if the accumulator size is large enough so that its LSB is smaller than those of all the inputs. Conversely, if a summand has an LSB smaller than LSB A (case of the bottommost summand on Figure 4), adding it to the accumulator entails a rounding error of at most 2 LSB A -1 . In the worst case, when adding n numbers, this error will be multiplied by n and invalidate the log 2 n lower bits of the accumulator -note that this worst-case situation is antagonist to the one leading to the worst-case overflow. A designer may lower LSB A to absorb such errors, an example is given in next section. Again, a practical maximum is an increase of 53 bits for 1 year of computation at 500MHz. However, in a global sense, it doesn't make much sense to consider this worst case situation, because when a summand is added exactly, it is usually the result of some rounding, so it carries an error of the order of its LSB, which it adds to the accumulator (by not adding bits lower than its LSB). These errors, which are not due to the accumulator, will typically dwarf the rounding errors due to the accumulator. They can be reduced by increasing w F , but this is of course outside of the scope of this article.

All considered, it is expected that no accumulator will need to be designed larger than 100 bits.

Finally, the proposed accumulator has sticky output bits for overflows in the summands and in the accumulator, so the validity of the result can be checked a posteriori.

A case study

In the inductance computation of [START_REF] Cret | FPGA-based acceleration of the computations involved in transcranial magnetic stimulation[END_REF], physical expertise tells that the sum will be less than 10 5 (using arbitrary units due to factoring out some physical constants), while profiling showed that the absolute value of a summand was always between 10 -2 and 2.

Converting to bit weight, and adding two orders of magnitude (or 7 bits) for safety in all directions, this defines MSB A = log 2 (10 2 × 10 5 ) = 24, MaxMSB X = 8 and LSB A = -w F -15 where w F is the significand width of the summands. For w F = 23 (SP), we conclude that an accumulator stretching from LSB A = -23 -15 = -38 (least significant bit) to MSB A = 24 (most significant bit) will be able to absorb all the additions without any rounding error: no summand will add bits lower than 2 -38 , and the accumulator is large enough to ensure it never overflows. The accumulator size is therefore w A = 24 + 38 + 1 = 63 bits.

Remark that only LSB A depends on w F , since the other parameters (MSB A and MaxMSB X ) are related to physical quantities, regardless of the precision used to simulate them. This illustrates that LSB A is the parameter that allows one to manage the accuracy/area tradeoff for an accumulator.

Accuracy measurements

Table 3 compares for accuracy and performance the proposed accumulator to one built using Xilinx CoreGen. To evaluate the accuracies, we computed the exact sum using multiple-precision software on a small run (20,000,000 summands), and the accuracy of the different accumulators was computed with respect to this exact sum. The proposed accumulator is both smaller, faster and more accurate than the ones based on FP adders. This table also shows that for production runs, which are 1000 times larger, a single-precision FP accumulator will not offer sufficient accuracy.

Table 4 provides other examples of the final relative accuracy, with respect to the exact sum, obtained by using an FP adder, and using the proposed accumulator with twice as large a significand. In the first column, we are adding n numbers uniformly distributed in [0,1]. The sum is expected to be roughly equal to n/2, which explains that the result becomes very inaccurate for n = 1, 000, 000: as soon as the sum gets larger than 2 17 , any new summand in [0,1] is simply shifted out and counted for zero. This problem can be anticipated by using a larger significand, or a larger MSBA Table 4. Accuracy of accumulation of FP [START_REF] He | Group-alignment based accurate floating-point summation on FPGAs[END_REF][START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machine[END_REF] numbers, using an FP [START_REF] He | Group-alignment based accurate floating-point summation on FPGAs[END_REF][START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machine[END_REF] adder, compared to using the proposed accumulator with 32 bits (M SBA = 20, LSBA = -11).

in the accumulator as we do. In the second column, numbers are uniformly distributed in [-1,1]. The sum grows as well (it is a random walk) but much more slowly. As we have taken a fairly small accumulator (LSB A = -11), for the first sums floating-point addition is more accurate: while the sum is smaller than 1, its LSB is smaller than -16. However, as more numbers are added, the sum grows. More and more of the bits of a summand are shifted out in the FP adder, but kept in the long accumulator, which becomes more accurate. Note that by adding only 5 bits to it (LSB A = -16 instead of -11), the relative error becomes smaller than 10 -10 in all cases depicted in Table 4: Again, LSB A is the parameter allowing to manage the accuracy/area tradeoff.

We have discussed in this section only the error of the long fixed-point accumulator itself (the upper part of Fig. 3). If its result is to be rounded to an FP [START_REF] He | Group-alignment based accurate floating-point summation on FPGAs[END_REF][START_REF] Shirazi | Quantitative analysis of floating point arithmetic on FPGA based custom computing machine[END_REF]) number using the post-normalisation unit of Figure 3, there will be a relative rounding error of at most 2 -17 ≈ 0.76 • 10 -5 . Comparing this value with the relative errors given in Table 4, one concludes that the proposed accumulator, with the given parameters, always leads to a result accurate to the two last bits of an FP(7,16) number.

ACCURATE SUM-OF-PRODUCTS

We now extend the previous accumulator to a highly accurate sum-of-product operator. The idea is simply to accumulate the exact results of all the multiplications. To this purpose, instead of standard multipliers, we use exact multipliers which return all the bits of the exact product: for 1 + w F -bit input significand, they return an FP number with a 2 + 2w F -bit significand. Such multipliers incur no rounding error, and are actually cheaper to build than the standard (w E , w F ) ones. Indeed, the latter also have to compute 2w F + 2 bits of the result, and in addition have to round it. In the exact FP multiplier, results do not need to be rounded, and do not even need to be normalised, as they will be immediately sent to the fixed-point accumulator. There is an additional cost, however, in the accumulator, whose input shifter is twice as large.

This idea was advocated by Kulisch [START_REF] Kulisch | Advanced Arithmetic for the Digital Computer: Design of Arithmetic Units[END_REF] for inclusion in microprocessors, but a generic DP version requires a 4288 bits accumulator, which manufacturers always considered too costly to implement. On an FPGA, one may design an application-specific version with an accumulator of 100-200 bits only. This is being implemented in FloPoCo, and Table 5 provides preliminary synthesis results for single and double precision input numbers, with the same 2w F accumulators as in Table 1. Some work is still needed to bring the FloPoCo multiplier generator to CoreGen level.

COMPARISON WITH RELATED WORK

Luo and Martonosi [START_REF] Luo | Accelerating pipelined integer and floating-point accumulations in configurable hardware with delayed addition techniques[END_REF] have described an architecture for the accumulation of SP numbers that uses two 64-bit fixedpoint adders. It first shifts the input data according to the 5 lower bits of the exponent, then sends it to one of the fixedpoint accumulators depending on the higher exponent bits. If these differ too much, either the incoming data or the current accumulator is discarded completely, just as in an FP adder. This design is much more complex than ours, and the critical path of the accumulator loop includes one 64-bit adder and a 3-2 compressor. The main problem with this approach, however is that it is a fixed design that will not scale beyond single-precision. Another one is that the detection of accumulator overflow may stall the operator, leading to a variable-latency design. The authors suggest a workaround that imposes a limit on the number of summands to add.

He et al [START_REF] He | Group-alignment based accurate floating-point summation on FPGAs[END_REF] have suggested group-alignment based floating-point accumulation. The idea is to first buffer the inputs into groups of size m (with m = 16 in the paper), while keeping trace of the largest exponent of a group. Then all the numbers in a group are aligned with this largestmagnitude number, and added in a fixed-point accumulator larger than the significands (30 bits for SP numbers). All this can be run at very high frequency. Then, these partial sums are fed to a final stage of FP accumulation that may run at 1/m the frequency of the first stage, and may therefore use a standard unpipelined FP adder. The authors also show that this accumulator is more accurate than one based on a standard FP adder. Again, this is a very complex design (for SP, 443 slices without the last stage, 716 with it). Besides, the frequency of the group accumulator will not scale well to higher precisions without resorting to techniques similar to our partial carry save.

Zhuo and Prasanna [START_REF] Zhuo | High performance linear algebra operations on reconfigurable systems[END_REF], then Bodnar et all [START_REF] Bodnar | Floating-point accumulation circuit for matrix applications[END_REF] have described high-throughput matrix operations using standard FP adders scheduled optimally thanks to additional buffers or memories. Performance-wise, this approach should be comparable to ours. The advantage of our approach is its better genericity and its finer control of the accuracyperformance tradeoff.

CONCLUSION AND FUTURE WORK

The accumulator design presented in this article perfectly illustrates the philosophy of the FloPoCo project: Floatingpoint on FPGA should make the best use of the flexibility of the FPGA target, not re-implement operators available in processors. By doing so, one comes up with a fairly simple design which can be tailored to be arbitrarily faster and arbitrarily more accurate than a naive floating-point approach, without requiring more resources. This is done in an application-specific manner, and requires the designer to provide bounds on the orders of magnitudes of the values accumulated. We have shown that these bounds can be taken lazily. In return, the designer gets not only improved performance, but also a provably accurate accumulation process. We believe that this return is worth the effort, especially considering the overall time needed to implement a full floating-point application on an FPGA.

The challenge is now to integrate such advanced operators in the emerging C-to-FPGA compilers. A profilingbased approach might be enough to set up the proposed accumulator. In parallel, FloPoCo should be extended with tools that help a designer set up a complete floating-point datapath, managing synchronisation issues but also accuracy ones. Such tools are still at the drawing board.

Fig. 1 .

 1 Fig. 1. Iterative accumulator

Fig. 2 .

 2 Fig. 2. A typical floating-point adder (w E and w F are respectively the exponent and significand sizes)

Fig. 3 .

 3 Fig.3.The proposed accumulator (top) and postnormalisation unit (bottom). Only the registers on the accumulator itself are shown. The rest of the design is combinatorial and can be pipelined arbitrarily.

Fig. 4 .

 4 Fig. 4. Accumulation of floating-point numbers into a large fixed-point accumulator

sign and exception handling

				X	Y
			wE + wF + 3	wE + wF + 3
			swap/difference
			X	Y	EX -EY
		wE + wF + 3	wE + wF + 3	wE
	EX wE	MX wF + 1	MY wF + 1		+/-		MY wF + 1 shift
		wF + 3	M Z		EX	MX wF + 1	/	wF + 4 M Y
		LZC			wE	M Z	wF + 4
		k log (wF + 3)	FZ shift/round			wF + 1 FZ round
	EZ wE + 1		wF + 1				far path
	close path		close/far
				wE + wF + 2
					final normalization
							wE + wF + 3
						Z

Table 1 .

 1 summand (w E , w F )CoreGen FP adder (w E , w F ) 2w F accumulator, MaxMSB X = 1 2w F accumulator, MaxMSB X = MSB A Compared synthesis results for an accumulator based on FP adder, versus proposed accumulator with MSB A = w E , LSB A = -w E , all targetted for 400MHz on a VirtexIV. suggested by Hossam, Fahmy and Flynn[START_REF] Fahmy | The case for a redundant format in floating point arithmetic[END_REF] for use as an internal representation in processor FPUs. For k = 1 one obtains a standard carry-save representation, but larger values of k are prefered as they take advantage of dedicated carry logic while reducing the register overhead. The FloPoCo implementation computes k out of the target frequency. For illustration, k = 30 allows to reach 400MHz on VirtexIV. The additional hardware cost is just the few additional registers -1/4 more in our figure, and 1/30 more for 400MHz accumulation.

	(7,16)	304 slices + 1 DSP, 12 cycles @ 359 MHz	129 slices, 8 cycles @ 472 MHz	176 slices, 9 cycles @ 484 MHz
	(8,23) SP	317 slices + 4 DSP, 16 cycles @ 450 MHz	165 slices, 8 cycles @ 434 MHz	229 slices, 9 cycles @ 434 MHz
	(10,37)	631 slices + 1 DSP, 14 cycles @ 457 MHz	295 slices, 10 cycles @ 428 MHz	399 slices, 11 cycles @ 428 MHz
	(11,52) DP	771 slices + 3 DSP, 15 cycles @ 366 MHz	375 slices, 11 cycles @ 414 MHz	516 slices, 12 cycles @ 416 MHz

Table 3 .

 3 Compared performance and accuracy of different accumulators for SP summands from[START_REF] Cret | FPGA-based acceleration of the computations involved in transcranial magnetic stimulation[END_REF].

		accuracy	area	latency
	SP FP adder acc	1.2 • 10 -3	317 sl, 4 DSP 16 @ 450 MHz
	DP FP adder acc 2.8 • 10 -15	771 sl, 3 DSP 15 @ 366 MHz
	proposed acc	2.0 • 10 -16	247 sl	10 @ 454 MHz

  CoreGen, SP ×, SP + 484 sl + 8 DSP, 26 cycles @ 366 MHz ours, SP ×, DP acc 319 sl + 4 DSP, 13 cycles @ 363 MHz CoreGen, SP ×, DP + 973 sl + 7 DSP, 26 cycles @ 366 MHz CoreGen, DP ×, DP + 1241 sl + 19 DSP, 37 cycles @ 366 MHz ours, DP ×, 105-bit acc 1441 sl + 9 DSP, 23 cycles @279 MHz

Table 5 .

 5 Preliminary results for Sum-Of-Product operators

www.ens-lyon.fr/LIP/Arenaire/Ware/FloPoCo/

This work was partly supported by the XtremeData university programme, the ANR EVAFlo project and the Egide Brâncus ¸i programme 14914RL.


Olivier Déforges
email: odeforge@insa-rennes.fr

Marie Babel
email: mbabel@insa-rennes.fr

Laurent Bédat
email: lbedat@insa-rennes.fr

Véronique Coat
email: coat@insa-rennes.fr

Scalable Lossless and Lossy Image Coding

published or not. The documents may come

INTRODUCTION

Despite many drawbacks and limitations, JPEG is still the most commonly-used compression format in the world. JPEG2000 overcomes this old technique, particularly at low bit rates, but at the expense of a significant complexity overhead. Therefore, the JPEG normalization group has recently proposed a call for proposals on JPEG-AIC (Advanced Image Coding) in order to look for new solutions for still image coding techniques [START_REF]Jpeg-aic : scope and evaluation[END_REF]. Its requirements reflect the earlier ideas of Amir Said [START_REF] Said | Reversible image compression via multiresolution representation and predictive coding[END_REF] for a good image coder: compression efficiency, scalability, good quality at low bit rates, flexibility and adaptability, rate and quality control, algorithm unicity (with/without losses), reduced complexity, error robustness (for instance in a wireless transmission context) and region of interest decoding at decoder level. A. Said and the JPEG committee should add additional functionalities such as image processing at region level, both at the coder or the decoder. The LAR (Locally Adaptive Resolution) tries to address all these features. In [START_REF] Déforges | Color LAR codec: a color image representation and compression scheme based on local resolution adjustment and self-extraction region representation[END_REF], we proposed an original scheme able to perform efficient lossy compression, enabling an unusual hierarchical region representation (without any shape description). Then, in [START_REF] Babel | Interleaved S+P Pyramidal Decomposition with Refined Prediction Model[END_REF], we presented an extension of a more efficient scalable multi-resolution solution in terms of both lossy and lossless compression. New improvements have been shown by the introduction of the Reversible Walsh Hadamard Transform (RWHaT) in [START_REF] Deforges | The RWHaT+P for an improved lossless multiresolution coding[END_REF]. This paper addresses further advances though the concept of inter-coefficient classification, focusing on the lossless compression context.

Section 2 gives an overview of the LAR coding method and the Reversible Walsh Hadamard Transform. These two concepts lead to the introduction of the RWHaT pyramid which can be improved when adding a prediction step. Both of these features are presented in section 3. Significant further improvements in terms of lossless compression efficiency are introduced in section 4. Finally, section 5 concludes the paper.

THE LAR AND RWHAT METHODS

LAR overview

The basic concept of the LAR method is that local resolution should be adapted to suit local activity. Also assuming that an image consists of global information and local texture, we firstly proposed a two-layer, content-based codec, both relying on a quadtree partition [START_REF] Déforges | Color LAR codec: a color image representation and compression scheme based on local resolution adjustment and self-extraction region representation[END_REF]. The first layer, called the FLAT LAR, encodes the global information at block level representation. The additional second layer enables texture compression within blocks using variable block size DCT/Hadamard transforms. Therefore, the method provides natural SNR scalability. The block sizes are estimated through a local morphological gradient. The direct consequence is that the smallest blocks are located round the edges whereas large blocks map homogeneous areas. This being so, the main feature of the FLAT coder consists of preserving contours while smoothing homogeneous parts of the image. This characteristic is also exploited to get a free hierarchical region representation: from the low bit-rate image compressed by the FLAT LAR, both coder and decoder can perform a segmentation process by iteratively merging blocks into regions. A direct application is then Region Of Interest (ROI) enhancement, by first selecting regions and enabling second layer coding only for the relevant blocks.

The RWhaT Transform

In order to enable spatial scalability and lossless compression, we studied different pyramidal representations. Lossless compression means that the construction and decomposition steps lead to a reversible process, resulting in the exact recovery of each element of the original source data. Moreover, the quadtree partition requires dyadic decomposition. Some integer wavelet kernels fit these requirements. For instance, Lux [START_REF] Lux | A novel set of closed orthogonal functions for picture coding[END_REF] proposed a modified version of the 1D Walsh Hadamard Transform (WHT), which has been popularized by Said [START_REF] Said | Reversible image compression via multiresolution representation and predictive coding[END_REF] and known as the "S" transform or "Haar integer wavelet transform" [START_REF] Komatsu | Lossless 2D Discrete Walsh-Hadamard Transform[END_REF]. The S transform is currently known as one of the best integer wavelet bases for reversible compression at the present time [START_REF] Abhayaratne | Lossless and Nearly Lossless Digital Video Coding[END_REF]. We introduced the Interleaved S+P technique [START_REF] Babel | Interleaved S+P Pyramidal Decomposition with Refined Prediction Model[END_REF] as an improved version of S+P by increasing space prediction. Another solution for a dyadic reversible decomposition has been proposed with the RWHaT transform for 2 × 2 blocks [START_REF] Deforges | The RWHaT+P for an improved lossless multiresolution coding[END_REF]. The conventional 2D W HT2×2 basis is defined as follows:

W HT2×2 Matrix: = 1 √ 2 1 1 1 -1 . (1)
Let U2×2 be the input block with:

U2×2 = u0 u1 u2 u3 . (2)
The block transformed Z2×2 by the 2D WHT is defined by:

Z2×2 = W HT2×2(U2×2) = W2×2 U2×2 W2×2 = z0 z1 z2 z3 = 1 2 u0 + u1 + u2 + u3 u0 + u1 -u2 -u3 u0 -u1 + u2 -u3 u0 -u1 -u2 + u3 .
(3) Clearly, if the sum of ui values is odd, then normalization by 2 produces an integer value. On the other hand, if the sum is even, the normalization induces a real value, leading to an ambiguous rounding operation. The solution consists of controlling the rounding step so that the inverse transform can distinguish between odd or even sum cases: for an even sum case, the rounding operations are applied such that the sum of the transformed coefficient is also even. A reversible transform then becomes possible, first by performing the inverse rounding operation to recover actual real transformed values, and secondly by carrying out the inverse 2D WHT transform [START_REF] Deforges | The RWHaT+P for an improved lossless multiresolution coding[END_REF].

RWHAT + P PYRAMID

Notations: I(i, j) denotes the pixel in an image I with the coordinates (i, j), I(b N (i, j)) the block b N (i, j) in I including the set of pixels {I(N.i, N.j), . . . I(N.i + -1, N.j + N -1)}.

The RWHaT pyramid

Let {Y l } Lmax l=0 be the multiresolution representation of an image I of size Nx × Ny, where Lmax is the top of the pyramid and l = 0 the full resolution image. As for the conventional W HT2×2 case, we iteratively construct the pyramid gathering four blocks to form a mean block at the upper level:

l = 0, Y0(i, j) = I(i, j); l > 0, Y l (i, j) = 1 4 1 k=0 1 m=0 Y l-1 (2x + k, 2y + m) (4) with 0 ≤ i ≤ N l
x , 0 ≤ j ≤ N l y , where N l x = Nx/l and N l y = Ny/l.

The top-down decomposition of the pyramid consists of encoding the WHT transformed block Z l (b 2 (i, j)) of each input block Y l (b 2 (i, j)). From ( 3) and ( 4), it is clear that the DC component of each block is unambiguously reconstructed from the upper level plus an additional bit called z 0 l . This bit is separately encoded from the other coefficients. Z l+1 denotes the WHT transformed block of Y l (i, j) and Żl (i, j) the corresponding block with z 0 l as DC value. The reconstruction from the previous level and current WHT transform is given by:

Ỹl (b 2 (i, j)) = EXP AN D[Y l+1 (i, j)) + Ỹl (b 2 (i, j)] with Ỹl (b 2 (i, j)) = W HT -1 2×2 Żl b 2 (i, j) . (5)
The EXP AN D function only duplicates a node value from the tree into its four sons. At this step, we get a common pyramidal representation and encoding process based on W HT2×2 trans-form, but with the exception of a possible lossless decomposition ( Ỹ = Y in this case). The top level has been encoded for the two methods by a simple DPCM. The common pyramid representation involving a dyadic decomposition is generally associated with a multilevel quadtree partition QP [2 Lmax ...2 l ] . Level l of the pyramid also specifies the finest resolution. More generally, we consider here a global quadtree partition of the image QP [Nmax...N min ] according to allowed block sizes. The parameter N l ∈ [Nmax . . . Nmin] gives the upper limit of block sizes to be decomposed at level l of the pyramid. For instance, a global partition QP [64...2] leads to the encoding of the representation from sizes 64 to 2, while N0 = 8 means that blocks of sizes 8, 4 and 2 will be decomposed at level 0. Finally, the decomposition parameter Lmin specifies the last level to be encoded: for all levels lower than Lmin, the value of all nodes in the pyramid are expanded only. The two pyramid scans presuppose a conditional decomposition. Assuming that A l (b 2 (i, j)) is the activity (gradient value) of block Y l (b 2 (i, j)), the conditional decomposition during the first pass is defined by:

If A l (b 2 (i, j)) < T h grad then encode Y l (b 2 (i, j)), else EXP AN D [Y l-1 (b 2 (i, j))] (6)
with typical T h grad = 30. Figure 1 shows the global coding based on this model called C1 coder, and restricted to the first pass. Blocks that are not decomposed during the first pass are processed during the second one.

The use of quadtree decomposition enables the approach to act as an "objective context modeling", separating the error prediction laws characterized by high entropy for the first pass, and low entropy for the second pass.

Prediction in the RWHaT pyramid

A prediction step can precede the transform operation to increase compression. Interpolation can also give a suitable means of providing a full resolution image from a decoded level l > 0. Even if functions have different objectives, we proposed in [START_REF] Deforges | The RWHaT+P for an improved lossless multiresolution coding[END_REF] a unified estimation solution. The main advantage is that, under particular conditions, the image can be progressively reconstructed from a previously interpolated level which also acts as prediction step. Yl (b 2 (i, j)) and Ỹl (b 2 (i, j)) denote the estimated and reconstructed blocks of Y l (b 2 (i, j)) respectively. The estimation

Estim Estim Estim ⊗ + - for C 3 Only ⊕ WHT2x2 -1 WHT 2x2 Coder Entrop Expansion . Exp(Y ) l+1 [ .. l N [Nmax Ỹ l Coder C 2 & C 3 Y l+1 ~Exp(Y ) l+1 l ≥ Lmin l< Lmin Ỹ l [N l ...N min ] Exp(Ỹ l+1 [N l ...N min ] )  l [N l ...N min ] Ỹ l [N max ...N l [ E l [N l ...N min ] Y l [N l ...N min ]
Fig. 2. Pyramidal coder with prediction step: C2 and C3 process consists of rebuilding unknown values around their block mean in a 2D context, according to a linear function. Two versions have been considered: one exploiting only inter-level information (C2), the other one using both inter and intra-level information (C3). Figure 2 shows new coding schemes including the estimation stage. The C2 coder uses only inter-level relationships and can be useful for progressive reconstruction. The C3 coder also makes use of reconstructed values at the current level and, obviously, leads to better compression.

Both SNR and spatial scalabilities are displayed in figure 3 for a lossless compression process. Six levels of decomposition are considered here, implying eleven successive streams ( 1 + 2 × 5). At low bit rates, visual distortions are essentially due to a blurring effect and are less disturbing than blocks or ringing artefacts. At the end of the first pass, reconstructed images present perfect accuracy upon strong contours while homogeneous areas appear as smooth regions.

INTER-COEFFICIENTS CLASSIFICATION

The a priori estimation of entropy laws for error coefficients is a crucial step for the efficiency of a lossless coder. In CALIC [START_REF] Wu | A Context-based, Adaptive, Lossless/Nearly-Lossless Coding Scheme for Continuous-Tone Images (CALIC)[END_REF], which is still one of the best available reversible compression techniques, Wu has introduced the concept of "context modeling" to characterize coefficients from the activity observed in their neighborhood. This method has been widely used in lossless compression context.

Let zi (0 ≤ i ≤ 3) be the transformed coefficients of the estimated block values. δz i = zi -zi denotes the associated errors. As previously mentioned, the exact reconstruction of the DC value only requires one bit, and is considered separately. Consequently, only z1, z2 and z3 have to be encoded. The context modeling approach suggests the estimation of local energy ∆i from the estimated value and previous locally-decoded data. To avoid context dilution problems, ∆i is generally quantized into L levels. The a priori estimation of δz i then consists of trying to maximize the conditional probability function p(δz i |∆i). This means that the various quantization levels have to be adapted to each image. In a lossless context and with full resolution processing, the number of encoded values is sufficient for an in-line estimation. However, in a multiresolution scheme, the top levels induce only a few values and they are not sufficient to enable such an approach. Another difficulty occurs when considering a lossy coding context, with quantized values. The problem maximizes the conditional probability p(δ ẑi |∆i) (with ẑi the quantized value of zi) and no more p(δz i |∆i). δ ẑi is, of course, dependent on the quantization factor qz, and if we consider a linear quantization, we can write

p(δ ẑi |qz) = p(qzδ ẑi ≤ δz i < (qz + 1)δ ẑi ). (7)
p(δ ẑi |qz) is an increasing function of qz. As qz and ∆i are independent, p(δ ẑi |(∆i, qz)) has to be maximized. When ∆i is quantized on L levels, and qz on Qn levels, they define L Qn contexts, making the estimation process much more complex.

To avoid the context quantization problem, we have introduced the concept of inter-coefficient classification. First, the energy function is defined according to the method features by considering:

1. the dynamics of the estimated value zi, 2. the difference δz 0 between the actual mean value of the block, z0, and the estimated value, z0.

The energy function ∆i associated with δz i is defined as follows:

∆i = w0 |δz 0 | + w1 |zi| , 1 ≤ i ≤ 3. (8)
The weights w0 and w1 have been empirically fixed at 0.5 and 1.

Inter-coefficient classification then consists of dividing the coefficients into equiprobable classes. Fixing the maximal number of classes at N b class, the definition of a class CLASSn is given by: (9) Each class constitutes a specific substream sent to the entropy coding layer. δz 3 coefficients present lower dynamics than δz 1/2 ones and are then scanned independently.

CLASSn = δz i |T∆ n-1 > ∆i ≥ T∆ n , 1 ≤ n ≤ N b class with T∆ 0 = ∞,
The C2 codec involves only inter-layer prediction, and thereby is slightly less efficient than C3. Nevertheless, it provides interesting SNR scalability when used with the classification coefficients method. In particular, as predicted coefficients do not depend on reconstructed data at the same level, values can be independently reconstructed. Thus, only coefficients belonging to CLASS0 for instance can be decoded when providing the class thresholds. This particular class contains coefficients with a higher dynamic, but it is also the most significant in terms of prediction error correction.

The first order entropy for the proposed method (C3 mode) is compared in Table 1 with the state of the art: S+P (scalable) and CALIC [START_REF] Wu | A Context-based, Adaptive, Lossless/Nearly-Lossless Coding Scheme for Continuous-Tone Images (CALIC)[END_REF] (non scalable) methods. Both the quadtree partition (Qd) and the inter-coefficient classification (Cl) techniques allow for major improvements. Combining the two techniques leads to an even better solution, significantly outperforming state of the art. This means that the two approaches act differently for data decorrelation and can be complementary.

CONCLUSION

The LAR codec is based on the concept of a two-layer encoding process. The first one is used to compress global information, while the second one adds local texture. One scalable and reversible solution has been proposed with the RWHaT+P pyramid. The two layers of the LAR codec are then implemented as two pyramidal decomposition passes, according to a quadtree partition. In terms of scalability, this process provides spatial and semantic SNR enhancement. In terms of compression, it acts as an objective coefficient classification estimated from the local activity of blocks. We propose here a complementary classification with a finer level of granularity, by considering activity measurements at both block and pixel levels. The inter-coefficient classification introduced in this paper avoids the common problems of context quantization inherent to such approaches, and significantly increases lossless compression efficiency. At the present time, this method also outperforms previous LAR performances for lossy coding. Future work will concentrate on the definition of a suitable entropy layer coder as the present system involves only a simple arithmetic one. Two other major new features have been also presented, leading to additional decorrelation, and further significant improvements i.e. content-based pyramidal decomposition and intercoefficient classification. The global scalable coding scheme surpasses both S+P and CALIC. It also provides locally adaptive multiresolution representations for the rebuilding of good quality images at low bit rates.

Fig. 1 .

 1 Fig. 1. Simple pyramidal coder: C1

  T∆ N b class = 0 and such as ∀n ∈ {0 . . . N b class} , card(CLASSn) N b class = const.

a)Fig. 3 .

 3 Fig. 3. Scalable lossless coding on "Barbara2" with partition QP [64...2]

Table 1 .

 1 2] Comparison of the proposed approaches with state-of-the-art methods S, CALIC and S+P . First-order entropy (bit/pixels).

	Entropy (bpp)


Min Won Lee

Laurent Larger

Jean-Pierre Goedgebuer

Transmission System Using Chaotic Delays Between Lightwaves

Keywords: Chaos, coherence modulation, encryption, unbalanced electrooptic Mach-Zehnder

We report a chaos-based transmission system that uses optical path differences between wavegroups produced by a short coherence source. Chaos synchronization is demonstrated in delay-differential equations involving optical delays larger than the coherence length of the source. The particular feature is all-optical chaos subtraction for message extraction and a very high masking efficiency.

I N RECENT years, secure communication systems based on chaotic dynamics have known an increasingly interest in electronics [START_REF] Pecora | Synchronization in chaotic systems[END_REF], [START_REF] Cuomo | Circuit implementation of synchronized chaos with applications to communications[END_REF], optics [START_REF] Vanwiggeren | Optical communication with chaotic waveforms[END_REF], [START_REF] Fujino | Experimental synchronization of chaotic oscillations in external-cavity semiconductor lasers[END_REF] and optoelectronics [START_REF] Larger | Optical encryption system using hyperchaos generated by an optoelectronic wavelength oscillator[END_REF], [START_REF] Uchida | Chaotic on off keying for secure communication[END_REF]. Most investigations have focused on the problem of chaos synchronization [START_REF] Peng | Synchronizing hyperchaos with a scalar transmitted signal[END_REF] with low-and high-dimensional chaotic systems, and recently on the problem of anticipating chaos [START_REF] Voss | Anticipating chaotic synchronization[END_REF]. Also recently, all-optical chaos synchronization and signal encryption/decryption with signal bandwidths up to the gigahertz range have been demonstrated experimentally from chaotic lasers with external feedback [START_REF] Fischer | Synchronization of chaotic semiconductor laser dynamics on subnanosecond time scales and its potential for chaos communication[END_REF], [START_REF] Tang | Message encoding and decoding at 2.5 Gbits/s through synchronization of chaotic pulsing semiconductor laser[END_REF]. All of those investigations have been sustained by the need of exploring mechanisms by which coupled simple nonlinear elements can entrain one another, one of the last example in this area being the investigation of delay-coupled semiconductor lasers [START_REF] Sivaprakasam | Message encoding and decoding using chaotic external-cavity diode lasers[END_REF]. The basic idea behind the use of chaos to transmit messages is to mask a message within a chaotic carrier before transmission, which is a form of spreadspectrum communication, and makes the presence of the signal difficult to detect in either the time or frequency domain, although breaking methods have been recently proposed [START_REF] Geddes | Extraction of signals from chaotic laser data[END_REF], [START_REF] Hegger | Identifying and modeling delay feedback systems[END_REF]. The masking efficiency is typically 15 dB in most of the present-day chaotic systems.

Unlike optical chaotic encryption systems reported to date, where chaos synchronization has been observed in the chaotic behavior of laser wavelength and laser intensity, we report in this paper the demonstration of chaos synchronization based on an unusual dynamical variable, the optical path-difference (OPD) between wavegroups of light produced by a short coherence source in coherence-modulated systems. Recent studies of chaos in coherence modulation have shown specific routes to chaos, including forward and reverse Hopf bifurcations, and their fundamental relevance [START_REF] Larger | Chaos in coherence modulation: bifurcations of an oscillator generating optical delay fluctuations[END_REF]. We present the synchronization process M. W. Lee is with the University of Wales, Bangor, Gwynedd LL57 2UW, Wales, U.K.

L. Larger and J.-P. Goedgebuer are with the University of Franche-Comté, 25030 Bensançon, France. between two such chaos generators and report numerical and experimental results about its dependence on the chaotic keys.

The encryption system (Fig. 1) consists of a transmitter (subscript 1) and a receiver (subscript 2), which are built symmetrically. The message to be transmitted is . The transmitter consists of a short coherence source [a superluminescent diode (SLD)], and an integrated unbalanced Mach-Zehnder modulator driven by a feedback loop. The feedback loop is formed by a second unbalanced MZ , a photodetector , a delay line (time delay s), an electronic amplifier (gain dB), and a low-pass filter with a time response s (which stands for the time constant of the setup and corresponds to a cut-off frequency ). The continuous wave (CW) SLD operates at m and with a coherence length m. Due to its low temporal-coherence degree, its emission may be regarded as a random sequence of wavegroups of coherence length and whose correlation function is related to the temporal coherence degree of the source [START_REF] Born | Principles of Optics[END_REF]Ch. 10]. The interferometers and have identical static OPDs, m greater than the coherence length. The optical beam at the output of is therefore formed by random pairs of wavegroups that do not produce any detectable intensity modulation since their extension is much smaller than the OPD separating them (for more details about coherence modulation, we refer the reader to [START_REF] Larger | Chaos in coherence modulation: bifurcations of an oscillator generating optical delay fluctuations[END_REF]). Electrooptically modulating results in OPD variations around that are linearly related to the modulation voltage through , where V is the half-wave voltage of . The OPD at output is , with .

A part of this coherence-modulated beam is sent to the receiver, while the other part is launched in . then converts nonlinearly the optical delay into a detectable output intensity modulation given by [START_REF] Larger | Chaos in coherence modulation: bifurcations of an oscillator generating optical delay fluctuations[END_REF] (1) where nW is the average optical power at output. After detection by , the signal is time delayed by through the electronic delay line, filtered by the low-pass filter and serves as the feedback signal applied to . The message is added to the feedback signal [START_REF] Udaltsov | Communicating with optical hyperchaos: information encryption and decryption in delayed nonlinear feedback systems[END_REF] with its amplitude (300 mV) at least ten times smaller than (8.6 V). The amplitude of the feedback voltage is adjusted to be twice higher than . Under these conditions, can be expressed as the addition of the OPD induced by the feedback signal and the OPD induced by with . Then, the transmitter equation, which rules the fluctuations of the OPD of is a delayed differential equation (DDE) that can be shown to be [START_REF] Larger | Chaos in coherence modulation: bifurcations of an oscillator generating optical delay fluctuations[END_REF] (2)

In ( 2), is a nonlinear function derived from (1) and expressed as [START_REF] Vanwiggeren | Optical communication with chaotic waveforms[END_REF] where is the bifurcation parameter and . Finally, the message is encrypted at the transmitter into a chaotic coherence-modulated beam formed by wavepackets separated by a chaotically modulated OPD larger than the coherence length of the source. Under these conditions, we do not have any detectable intensity modulation in the transmitter beam.

We have already described the complicated dynamics and routes to chaos corresponding to [START_REF] Cuomo | Circuit implementation of synchronized chaos with applications to communications[END_REF], showing that forward and reverse Hopf bifurcations can be obtained [START_REF] Larger | Chaos in coherence modulation: bifurcations of an oscillator generating optical delay fluctuations[END_REF]. For the parameters used in the transmitter, the solutions of (2) are hyperchaotic, with the Lyapunov dimension evaluated from Kaplan-Yorke's conjecture [START_REF] Farmer | Chaotic attractors of an infinite-dimensional dynamical system[END_REF], [START_REF] Kaplan | Functional Differential Equations and Approximation of Fixed Points[END_REF] found to be about 76.

The receiver is constructed from the same components as those in the transmitter. A part of the transmitted beam is sent into the loop of the receiver. This loop is formed by , which features a static OPD and a coherence demodulation function , a delay line , an amplifier (gain ) and a low-pass filter the solutions of (2), when there is no message ( , ) at the transmitter, meaning that we have chaos replication in that case.

The other part of the transmitted coherence-modulated beam is launched in . operates as a coherence demodulator that features a chaotically modulated OPD with . That branch is aimed at optically carrying out a subtraction between the chaotic fluctuations in the transmitted beam and those produced by , to retrieve the original message . Noting that and is small, the output intensity of is given by (5) The original message is then directly retrieved at the receiver output. Since matching the transmitter and receiver parameters is one key condition, it is important to know its impact on the recovery of the message. For that, we first evaluated the sensitivity of the synchronization process from numerical simulations using ( 2) and ( 4), and introducing voluntarily mismatches between the different parameters at the receiver. The computations were done with no message and evaluating the power spectrum of the synchronization error signal (which is also chaotic), normalized to the rms level of the chaotic carrier , as a function of the relative parameter mismatch expressed as a percent. Results of the simulations are plotted in Fig. 2 for mismatches ranging from 2% to 2% for and for . In Fig. 2, the magnitude of the decoding noise (expressed in decibels) is encoded into gray scales, where dark gray corresponds to a very small synchronization error (a good message recovery), whereas white indicates a very bad synchronization, and hence a bad message decoding. Typically, following the usual standards in telecommunications, a good decoding quality is obtained when the noise level is smaller than 40 dB that corresponds here to 20 dB for the signal-to-noise ratio and 20 dB for the chaos-to-signal ratio. A contour for those 40 dB is indicated in the figure. Computation shows that synchronization is very sensitive to a time delay mismatch [Fig. 2 (b)]. In that case, the decoding noise can reach 0 dB when the mismatch is a few percent, meaning that the message is completely altered by noise. In contrast, the synchronization process is much less sensitive to a -mismatch, as shown in Fig. 2 (a). Similar computations can be done for parameters and , showing that an intermediate sensitivity is obtained, compared to the one presented for parameters and .

To confirm those numerical simulations, we then investigated experimentally the synchronization process by varying the value of the time delay mismatch, all the other parameters being closely matched one with each other. The results are shown in Fig. 3. Fig. 3(a) shows, for a 0.01% -mismatch, the Lissajou representation of the driving signals of and , which are linearly related to and at the transmitter and receiver, respectively. In spite of their multiple positive Lyapunov exponents [START_REF] Peng | Synchronizing hyperchaos with a scalar transmitted signal[END_REF], these signals show a very good synchronization, and the signal-to-noise ratio was then measured to be 36.5 dB from this figure. Fig. 3(b) corresponds to only a 1% -mismatch, resulting in a large synchronization error and non correlated chaotic signals. The fact that the most critical parameter to consider is the time delay in the feedback loop has already been observed with other chaotic schemes [START_REF] Cuenot | Chaos shift keying with an optoelectronic encryption system using chaos in wavelength[END_REF].

In a third step, we adjusted the parameters to meet the conditions of optimal matching (all the parameters matched to less than 1%), and we introduced a message in the feedback loop of the transmitter. The encryption and decryption results are illustrated in Fig. 4. The residual intensity modulation was measured to be smaller than 1% of the average intensity. This residual intensity modulation is due to the low mutual coherence degree of the wavegroups of light at output, and could be made even smaller by using an OPD in greater than 260 m. We also show the spectrum of the residual intensity modulation of the coherence-modulated beam [Fig. 4 (c)], which provides information that an eavesdropper can detect when tapping the transmission fiber. No detectable peak at the 1-kHz message frequency can be seen in the spectrum, down to 70 dB, which is the masking efficiency (the masking efficiency in chaotic transmission systems is expressed as the message-to-chaos ratio in the transmitted signal). The message is therefore strongly hidden due to the chaotic coherence modulation encryption. Finally, the spectrum of the decoded message from the receiver is shown in Fig. 4 (d), in which the signal-to-noise ratio was measured to be 33 dB. In the case of a realistic broadband information message, we observed similar signal-to-noise ratio performances, which can be explained by the fact that the decoding noise does not depend significantly on the frequency inside the chaos generator bandwidth (see Fig. 2).

To summarize, we have reported for the first time the proofs of principle of an optical cryptographic system based on a chaotic coherence modulation. OPD fluctuations are used to mask an information message. The chaotic waveform exhibits a high Lyapunov dimension ( 76). Such a hyperchaos was successfully synchronized in a transmitter and receiver, and the impact of mismatches has been investigated numerically and experimentally. Encryption and decryption of a message was obtained experimentally with a 33-dB signal-to-noise ratio. Direct recovery of the message at the receiver was performed through a subtraction process performed optically (which seems to be specific to the coherence demodulation scheme), unlike other schemes in which an electronic subtractor is generally used. The system features an extremely high masking efficiency ( 70 dB), a value which is the best performance reported with chaotic systems so far. The system could operate in the gigahertz range with high-speed components [START_REF] Fischer | Synchronization of chaotic semiconductor laser dynamics on subnanosecond time scales and its potential for chaos communication[END_REF]. Another attractive feature of the system concerns the multiplexing capabilities of coherence modulation systems [START_REF] Hauden | Demonstration of a single source bidirectional fiber link using polarization insensitive LiNbO integrated coherence modulator[END_REF], in which each channel can be attached to a different OPD. Finally, the system can be potentially dedicated to optical telecommunications, since it makes use of fiber-based components and integrated Mach-Zehnder modulators.

Fig. 1 .

 1 Fig. 1. Experimental setup. SLD: super luminescent diode. MZ: unbalanced electrooptic Mach-Zehnder. PD: photodiode. TD: time delay. Amp: amplifier. LPF: low-pass filter.

Fig. 2 .

 2 Fig. 2. Spectrum of the synchronization error for various parameter mismatches. The contour line indicates a -40 dB noise level. (a) Response time mismatch. (b) Delay T mismatch.

Fig. 3 .

 3 Fig. 3. Experimental synchronization with a mismatch between T and T . (a) Receiver chaos versus emitter chaos for a 0.01% T -mismatch. (b) Same representation for a 1% T -mismatch.

  Fig. 4(a) shows the original and decoded message in the case of a 1-kHz sine waveform. The driving voltage of thus obtained (lower trace), and the intensity of the coherence-modulated beam carrying the message in the transmission line (upper trace), are shown in Fig. 4(b).

Fig. 4 .

 4 Fig. 4. Experimental encoding and decoding results. (a) Original and decoded 1-kHz sine-waveform. (b) Transmitted intensity and corresponding OPD fluctuations. (c) FFT of the transmitted intensity. (d) FFT of the recovered message.