Antoine Perquin
email: antoine.perquin@irisa.fr

Gwénolé Lecorvé
email: gwenole.lecorve@irisa.fr

Damien Lolive
email: damien.lolive@irisa.fr

Laurent Amsaleg
email: laurent.amsaleg@irisa.fr

Phone-Level Embeddings for Unit Selection Speech Synthesis

HAL is

Introduction

Unit selection speech synthesis concatenates pre-existing segments of recorded speech, producing high-quality, natural sounding, oral renderings of sentences [START_REF] Hunt | Unit selection in a concatenative speech synthesis system using a large speech database[END_REF]. This process optimises a target cost function selecting the units that best match the linguistic descriptions of the phonemes to synthesize. Quality results from the involvement of linguistic experts who carefully design that function. These methods, however, suffer from concatenation errors and cannot generalize outside the pre-recorded units. Furthermore, they necessitate extremely costly human expertise, which might not exist when targeting or adapting to other domains or languages.

Instead, this paper proposes an embedding-based method allowing to rely on euclidean distances between units when optimizing their selection. Not only is this cheap compared to human expertise, but embeddings also facilitate domain adaptation. This paper presents four models, with two of them resulting in phonelevel embeddings for unit selection. Three of them rely on a feed-forward Deep Neural Network (DNN) whereas the last one uses a Long Short Term Memory layer (LSTM). Experiments using a large French speech corpus show that the use of the embeddings outperforms expert-based unit selection.

A few remarks are in order. One, temporal information can be integrated in different ways. We compare its early or later integration in DNNs, and study the use of temporal dependencies brought by an LSTM layer. Two, in an attempt to understand if mitigating the effects of the curse of dimensionality on the embeddings is beneficial, three models use varying layers sizes at their core. Three, as far as we know, no objective metrics for assessing the quality of such embeddings exist. As a proxy, we evaluate their ability to predict appropriate acoustic features. Four, in contrast to state of the art contributions that work for English, we work with French. This is the opportunity to assess how the embedding-based approach performs on French.

The remainder of this paper is as follows. Section 2 presents related work. Section 3 discusses the handling of time and dimensionality in the embeddings. Finally, Section 4 compares our approach to an expert unit selection system.

Related Work

The Statistical Parametric Speech Synthesis (SPSS) approach uses vocoders to synthesize speech from acoustic parameters predicted by a so-called acoustic model based on linguistic features [START_REF] Black | Statistical parametric speech synthesis[END_REF]. Recent models are now based on DNNs [START_REF] Ze | Statistical parametric speech synthesis using deep neural networks[END_REF]. Beside pros and cons of this approach (good flexibility but low quality), the resulting prediction models are interestingly independent of any linguistic expertise.

This approach gave rise to different extensions. On the one hand, the Wavenet model proposes to integrate vocoding in the prediction model, i.e., the raw waveforms are directly predicted [START_REF] Van Den Oord | Wavenet: A generative model for raw audio[END_REF]. End-to-end approaches such as [START_REF] Wang | Tacotron: Towards end-to-end speech synthesis[END_REF] go further by predicting these waveforms directly from text, that is linguistic features are removed.

On the other hand, prediction models have also been used to transform linguistic features into other representations, based on which a new target cost for unit selection can be defined. In turn, the new representation space can be handled more easily than the symbolic space of the linguistic features. In [START_REF] Yan | Rich-context unit selection (RUS) approach to high quality TTS[END_REF], linguistic features are converted into acoustic ones based on hidden Markov models. The Kullback-Leibler divergence is then used as a target cost. Alternatively, [START_REF] Merritt | Deep neural networkguided unit selection synthesis[END_REF][START_REF] Wan | Googles next-generation real-time unit-selection synthesizer using sequence-to-sequence lstm-based autoencoders[END_REF] have proposed intermediate representations extracted from hidden layers of DNNs, also called embeddings.

In [START_REF] Merritt | Deep neural networkguided unit selection synthesis[END_REF], each phone is divided into 4 sections for which the mean and variance of the embedding of each frame is computed. The target cost function is the sum of the Kullback-Leibler divergence on the mean and variance of the 4 sections. This showed that there is no loss in the quality of speech when using embeddings instead of acoustic features. The subdivisions, however, still involve some human expertise.

In [START_REF] Wan | Googles next-generation real-time unit-selection synthesizer using sequence-to-sequence lstm-based autoencoders[END_REF], a phone level embedding is obtained with a multi-modal model composed of an acoustic encoder and a linguistic encoder. The target cost is the euclidean Fig. 1. Architecture of the 4 models. Layer sizes are reported. denotes the vector of linguistic features, d the phone duration, ti the time position in the phone for the i-th frame, ai the corresponding acoustic features, and N is the total number of frames in the phone.

distance in the embedding space. This work is very similar to ours, and was actually conducted at the same time as [START_REF] Perquin | Big deep voice: indexation de données massives de parole grâce à des réseaux de neurones profonds[END_REF]. Our work focuses on the effect of various time integration and information compression schemes, rather than those of multi-modality.

Handling of time and dimensionality

A key assumption in our work is that the quality of an embedding is correlated to the quality of the model from which it is derived, in our case acoustic models. In this section, different variants of acoustic models are studied, among which those proposed for unit selection. These models are objectively and perceptually evaluated to assess the behaviour of embeddings with respect to information compression and different strategies for time integration. This section first presents the models under study before describing the experimental dataset, then the objective and perceptual tests.

Models

Four models are studied, see Figure 1. The first model is a standard acoustic model (DNN-S), a simple feed-forward DNN, comparable to the one proposed in [START_REF] Wu | Improving trajectory modelling for dnn-based speech synthesis by using stacked bottleneck features and minimum generation error training[END_REF]. For a given phoneme, the model predicts acoustic features a i of the i-th frame based on the linguistic feature vector . Those linguistic features provide information about the phoneme, e.g., its identity, the one of its close neighbours, its position in the syllable/word/utterance it belongs to, etc. The timing information is encoded as two numerical features: the phone duration d in seconds and the relative position t i of the frame i inside the phone. This timing information is useful to take into account the dynamics of acoustic features when realising a phone. For the frame i, the output of the middle layer can be seen as the frame embedding e i .

The second model, DNN-B, has a similar architecture as DNN-S but introduces a bottleneck layer. This bottleneck is obtained by gradually decreasing then increasing the size of the hidden layers. e i is then a compression of the linguistic features . As a side effect, compressing the embedding space avoids the curse of dimensionality and allows for tractable similarity measures.

The limitation of both models is that the resulting embeddings would only represent the frame currently predicted, whereas phone-level embeddings are needed for unit selection. To solve this problem, we propose to postpone making use of the timing information (d and t i ) until after the embedding layer, as in [START_REF] Wan | Googles next-generation real-time unit-selection synthesizer using sequence-to-sequence lstm-based autoencoders[END_REF]. Applying this principle on DNN-B gives the model DNN-BP, P for Phone-level. For a given linguistic feature vector , we obtain a phone-level embedding e.

Finally, to attempt to model the timing dependency across frames, we propose to replace the layer after the embedding layer with an LSTM in DNN-BP to obtain the model LSTM-BP. While the previous models are trained on independent frames, LSTM-BP is trained on the full frame sequence of the considered phone. This method decreases the shuffling possibilities over the training set but could lead to better predictions, and thus better embeddings.

The implementation details of the different models are as follows:

-DNN-S: All hidden layers rely on the hyperbolic tangent (tanh) activation function. For all models, the output layer is with a linear activation.

Dataset and experimental setup

Our models were trained on a corpus corresponding to the reading of a French audio-book by a professional French speaker resulting in approximately 11 hours of speech for a total of 3300 utterances (approximately 390 000 phonemes). Speech is expressive (narration, acted dialogues), and sentences are complex (long sentences, formal register) due to the style of the audiobook's author (Marcel Proust). 105 utterances were held out for the listening test, while the rest was shuffled at the frame level (or phone level for model LSTM-BP) and then split into three sets : 90% for the training set and 5% each for the validation set and test set. About 110 linguistic features are considered for each phone. Categorical attributes represent information about quinphones, syllables, articulatory features, and part of speech for the current, previous and following words. They are encoded in one-hot. 34 other features are numerical, such as the position of the phone inside the word or the utterance. After encoding, the overall linguistic vector is of size 341. The linguistic features and the timing information were normalised to the range [0.01, 0.99]. Each linguistic feature was manually extracted, without automatic annotation.

The acoustic features, extracted using the WORLD vocoder [START_REF] Morise | WORLD: a vocoder-based high-quality speech synthesis system for real-time applications[END_REF], consist of a 60 dimension Mel-Generalized Cepstral coefficients (MGC) vector, a 5 dimension band-aperiodicity (BAP) vector and the fundamental frequency F 0 . Those features were extracted every 5 ms. The F 0 coefficient was linearly interpolated on unvoiced parts, a boolean attribute keeps track whether the frame was voiced or not and the logarithm was applied to F 0 . Finally, the deltas and delta-deltas were computed for MGC, BAP and F 0 . In total, the acoustic vector is of size 199. The acoustic features were centered and normalized to unit variance.

The implementation was done using Keras with TensorFlow. Training was done on a GTX 1080 Ti, over 100 epochs using RMSPROP with the mean square error as a loss function. The model weights with the best performance on the validation set were saved. Those models were trained using the true duration values.

Objective Evaluation

The 4 models were evaluated in an acoustic modelling perspective. Table 1 reports the quality of the predicted acoustic features according to the following measures:

-MCD : Mel-Cepstral Distortion on MGC coefficients.

-BAP : a distortion measure on BAPs.

-V/UV : Voiced/unvoiced error rate.

-RMSE(F0) : Root mean squared error on F 0 .

Those measures are computed between the acoustic features predicted by the DNNs and the reference acoustic features. The results from a state of the art acoustic model are also directly reported from [START_REF] Wu | Merlin: An open source neural network speech synthesis system[END_REF]. They were not trained on the same data, nor even the same language. They are presented for the sake of a sanity check.

First, our measures are higher than those from [START_REF] Wu | Merlin: An open source neural network speech synthesis system[END_REF], especially regarding F 0 and voicing error. We believe the high error on those two aspects is due to the high expressiveness of our speech corpus. Taking this into account, these results can be considered as acceptable. Second, by comparing model DNN-B and DNN-BP, we can see that the displacement of the timing information raises the error for all measures by only a small margin. Then, while one would expect DNN-S to lead to the best results, it appears that its performance is a bit worse than for DNN-B and DNN-BP. Possible explanations are the larger number of parameters for model DNN-S (more difficult to reach the global optimum for weights) or the larger number of layers in the DNN-B and DNN-BP. At least, the results show that the bottleneck does not hurt. On the contrary, LSTM-BP leads to the worst results. This is surprising since adding extra information about the previous frames should not degrade the prediction. In our opinion, these results mainly come from the fact that training data are organised as sequences of frames, thus reducing the diversity of observations during training, whereas the other models are trained on shuffled frames. Despite those bad results, we chose to keep LSTM-BP to evaluate the effect of an acoustic model quality on its embeddings and on the resulting synthesized speech.

Perceptual Evaluation

The previous objective measures are not perfect estimators of the overall speech quality. For example, a significant improvement in MCD does not necessarily translate to a perceptual improvement of the synthesized speech. Thus we want to observe the perceptual impact of embeddings on SPSS.

For each model, synthetic speech has been generated based on the predicted acoustic features using an SPSS approach. In practice, the WORLD vocoder is used. Long utterances were split into breath groups of 4-5 seconds. The timing information (d and t i ) was derived from a duration model (DNN with 6 hidden layers and tanh activation) predicting the number of frames N of a phone based on its linguistic feature vector . The mean absolute error of this model is 3.7 frames. d is derived by multiplying N by the sampling rate, t i is simply i N . A listening test was conducted with 21 French native speakers. They were asked to rate between 0 and 10 the overall quality of speech utterances. Synthetic speech coming from the 4 models was presented along with natural speech (NAT) and utterances vocoded based on the reference acoustic features (VOC). Each listener was given 10 in-domain utterances and 10 out-of-domain utterances. For in-domain evaluation, we used the 105 sentences left-out from the dataset, while for out of domain evaluation we used 100 phonetically balanced French sentences. Thus, every utterance was rated by at least two different listeners.

In accordance with the objective results, the listening test shows that our models do not perform well in SPSS mode, as can be seen on listeners reacted well to perfect SPSS (analysis-synthesis), giving a 7.2 mean score to the system VOC, they gave a mean score of around 3 to our first 3 models, and agreed that system LSTM-BP's productions were incomprehensible with a mean score of almost 0. This is not surprising since this system had a higher MCD than other systems, which is measured on a log-scale. There is no real statistically significant (p-value greater than 0.05) difference between systems DNN-S, DNN-B and DNN-BP which proves that the displacement of the timing information does not cause any perceptual loss in quality. Surprisingly, the natural speech received a couple of marks below 7. The corresponding utterances were perceived as having unnatural prosody because they were cut to be shortened.

Comparison with Expert Unit Selection

In this section, the proposed unit selection method is compared to a system where the target cost is defined based on expert knowledge.

Unit Selection Engine

Once phone-level embeddings have been extracted, we use them to guide the unit selection process. Before synthesis, for each phone in the database, the corresponding embedding is computed and stored. At synthesis time, a preselection reduces the set of candidate units to those corresponding to the same phoneme as the target, in order to reduce the computation time. target phone in the embedding space among the pre-selected units. Finally, the lattice of these nearest neighbours for each phone is decoded to find the sequence of units minimizing the sum of the target and concatenation cost. The target cost is defined as the euclidean distance in the embedding space between the candidate and target phone. For the expert system, the target cost was originally defined as a weighted sum of linguistic features and has since been improved over the years. The join cost for all systems is the same as in [START_REF] Lolive | The irisa text-to-speech system for the Blizzard challenge[END_REF], defined as a sum of euclidean distances on acoustic features between following candidate units.

Perceptive Evaluation

Since there is no proposed measure in the literature to evaluate the quality of an embedding in relation to speech synthesis, we directly address the subjective evaluation of the embeddings with listening tests. The results for in-domain and out-of-domain utterances can be found on Figure 3 and Figure 4 respectively.

For in-domain utterances, the 3 proposed models were all awarded really high grades (multiple times rated with a 10/10 score for each system) but also really low notes (up to 0 for system DNN-BP). However, on average, the expert system received a mean score of 5.4, system DNN-BP received 6.8 and system LSTM-BP received 6.6. While we cannot statistically distinguish the two systems with embeddings (p-value greater than 0.05), both are a statistically significant improvement over the expert system. Interestingly, even if the system LSTM-BP performed really poorly for SPSS, it receives good results for the unit selection paradigm. For out-of-context utterances, natural speech from the same voice as the one used in the database was unavailable. The expert system received a mean score of 5.9, system DNN-BP received 7.0 and system LSTM-BP received 6.4. This time, the difference between the two systems using embeddings is significant. However, they both perform without significant loss over in-domain synthesis. Surprisingly, the expert system appears to have been better rated on out-of-domain utterances. A probable explanation is the absence of natural speech as a higher baseline during the MUSHRA test. Still, the embedding systems remain significantly better than the expert one proving that the embeddings are general enough to be used in other domains.

EXP

Conclusion

In this paper, we proposed two models to extract phone-level embeddings in the context of DNN-driven unit selection. The models were compared to other DNNs on the task of acoustic modeling. The experiments highlighted that late integration of time and information compression (bottleneck) do not impact the quality of feature prediction, even though the use of LSTM did not seem conclusive yet. Then, experiments on unit selection showed that quality of synthesized speech based on embeddings resulting from our models is perceptually preferred over an approach with an expertly defined target cost. They also demonstrated that the proposed embeddings are general enough to be used in multiple domains.

Besides these results, our study also highlights the fact that the link between the acoustic feature prediction and the quality of unit embeddings is not clear.

For instance, the LSTM-BP model led to good unit selection results whereas it was very bad for SPSS. It would thus be interesting to further study how embeddings could be evaluated with objective measures, e.g. by analyzing the topological properties of the embedding spaces. Then, extensions of the proposed models should be tested to produce better embeddings, and to better understand dependencies across different types of information. Multi-modality, as in [START_REF] Wan | Googles next-generation real-time unit-selection synthesizer using sequence-to-sequence lstm-based autoencoders[END_REF], is a first direction. The integration of duration into the embeddings is another. Finally, our objective in the long term is to deal with large heterogeneous speech databases (different speakers, different languages, etc.). Apart from the previously raised questions, this would require to study the compliance of the embeddings with fast database searches, especially in regard to approximate nearest neighbours techniques.

Figure 2 Fig. 2 .

 22 Fig. 2. Results of the listening test for the models in statistic parametric mode for in-domain utterances.

Fig. 4 .

 4 Fig. 4. Results of the listening test for the models working in unit selection mode for out-of domain utterances.

Table 1 .

 1 Objective evaluation of the predictions of our models.

		MCD BAP V/UV RMSE(F0)
		(dB) (dB) (%)	(Hz)
	DNN-S	5.22 0.48 17.2	18.3
	DNN-B	5.06 0.35 12.6	17.9
	DNN-BP	5.09 0.36 13.7	18.2
	LSTM-BP	5.80 0.49 19.7	19.5
	DNN (reported from [11])	4.54 0.36 11.38	9.57

  Results of the listening test for the models working in unit selection mode for in-domain utterances.

		10
		8
	Rating	4 6
		2
		0
		NAT	EXP	DNN-BP	LSTM-BP
	Fig. 3.
					The number of
	candidate units is reduced further by searching the 25 nearest neighbours of the

Acknowledgments

This study has been realized under the ANR (French National Research Agency) project SynPaFlex ANR-15-CE23-0015.


Hugues De Lassus Saint-Geniès

Nicolas Brunie

Guillaume Revy

Meta-implementation of vectorized logarithm function in binary floating-point arithmetic

Keywords: logarithm function, vector micro-architectures, faithful rounding, automated design

Besides scalar instructions, modern microarchitectures also provide support for vector instructions. They enable to treat packed inputs (typically 4 or 8) in a single instruction. The challenge is now to write vector programs to support mathematical functions like sin, cos, exp, log, • • • which efficiently exploit those vector instructions. This article focuses on the design of vectorized implementation of log(x) function, and more particularly on its automation for different formats and micro-architectures. First it rewrites a classic range reduction in a branchless fashion so as to use at best recent micro-architecture features, like rcp (reciprocal) instruction, and to treat all inputs in the same flow. Second it details rigorously how to achieve "faithfully rounded" implementations. Third it shows how to automate this implementation process using the MetaLibm framework, on SSE/AVX and AVX2 supporting micro-architectures. Finally we illustrate that this process enables to achieve high throughput implementations for the binary32 and binary64 formats in a fully automated way.

I. INTRODUCTION

Nowadays, most modern micro-architectures embrace new hardware support for vector instructions in floating-point arithmetic, besides scalar instructions. They enable to treat packed inputs (typically 4 or 8) in a single instruction, thus potentially increasing the performance of floating-point programs. The challenge is now to use these instructions properly to write efficient vectorized implementations of mathematical functions. These are widely used in scientific programs manipulating floating-point computations, for which performance improvement is crucial (like in physics, for example).

In this article, we deal with the vectorized implementation of the logarithm function, log(x), in floating-point arithmetic, with a particular focus on its automation through the MetaLibm framework [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF], [START_REF] Brunie | Code generators for mathematical functions[END_REF]. It enables to describe the implementation of a function using a meta-language and to generate C codes optimized for different micro-architectures. For performance purpose, we target "faithfully rounded" implementations, that is, with 1-ulp accuracy. Even if, since 2008, the IEEE 754 standard recommends the correct rounding for scalar elementary functions like log(x), no output accuracy is required nor recommended for functions with vector arguments [START_REF]IEEE standard for floating-point arithmetic[END_REF], and "faithful rounding" is often sufficient in highperformance computing context, for example in physics [START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF].

Several vectorized implementations of mathematical functions have already been proposed. The commercial Intel SVML library 1 provides intrinsic support for vectorized mathematical functions, but it is not freely available. Other, noncommercial implementations have been proposed. In 2015, Intel contributed the GNU Libmvec to glibc 2.22. 2 AMD also contributed the proprietary ACML_MV (no longer available) and the AMD LibM. 3 All these libraries support a set of vectorized mathematical functions. GNU Libmvec claims reasonable testing to pass 4-ulp maximum relative error, while ACML_MV was announced returning results accurate to at most 1 ulp, but with unpredictable behavior for subnormal inputs. Unfortunately, GNU Libmvec is written in assembly code, which reduces portability. The VDT is another open source effort [START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF]: in terms of accuracy, log(x) function is announced with at most 2 bits different compared to AMD LibM. Two recent projects are Yeppp 4 and SLEEF, 5 the latter providing 3-ulp vectorized function versions. Finally in 2016, a new SIMD vector library is proposed [START_REF] Lauter | A new open-source SIMD vector libm fully implemented with high-level scalar C[END_REF]. Written in high-level C language, log(x) is accurate to at most 1.0825 ulp. Note that the accuracy of these implementations is not guaranteed since, to our knowledge, it is obtained by non-exhaustive testing.

We provide here high throughput "faithfully rounded" implementations of log(x), in the binary32 and binary64 formats, and with a guarantee on the output accuracy. We propose also a meta-implementation of this function in the MetaLibm framework, enabling to write different optimized binary32 versions. This work being parametrized by the precision and code generation being based on MetaLibm backends, this can easily be extended to any other precision. Our functions do not handle special inputs, producing unpredictable outputs in these cases, and the subnormal inputs are considered in the main flow. This results in branchless programs, which avoids an expensive fallback treatment for these inputs (like most of the vector libraries presented above do) at a small overhead.

Therefore the main contributions of the article are:

• A branchless range reduction, designed to use at best some special features available on SSE/AVX and AVX2 instruction sets, • A rigorous description of how to achieve "faithful rounding" in a guaranteed fashion,

• The meta-implementation of this range reduction in the MetaLibm framework, parametrized by the precision, and enabling to generate different function versions, optimized for various instruction sets and vector sizes. The article is organized as follows: Section II details the range reduction we use to achieve required accuracy, while Section III gives some guidelines on how to build accurate enough programs. Then Section IV provides some implementation details. In Section V, we present the work done in the MetaLibm framework. And Section VI gives some numerical results, before a conclusion in Section VII.

II. ALGORITHM FOR log(x) WITH FAITHFUL ROUNDING

We consider here that the input x is neither a special input (NaN, ±∞, ±0, or x ≤ 0), nor the particular input x = 1. All the floating-point numbers manipulated, as well as the floatingpoint operations carried out, are in precision p ≥ 2.

A. Range reduction

Let x = 1 be a positive floating-point number as defined in [START_REF]IEEE standard for floating-point arithmetic[END_REF] with the sign s = 0:

x = m • 2 e , (1)
where m ∈ [1, 2 -2 1-p ] and e ∈ {e min , • • • , e max }. Many techniques have already been designed for the implementation of log(x). They are generally based on table lookup and/or polynomial evaluations [START_REF] Tang | Table-driven implementation of the logarithm function in IEEE floating-point arithmetic[END_REF]- [START_REF] Revy | Automated design of floating-point logarithm functions on integer processors[END_REF]. Here let r i denote an approximation of 1/m. Then using log(x) = e • log(2) + log(m), the log(x) function is computed as follows:

log(x) = (e + τ ) • log(2) -log(2 τ • r i ) + log(1 + u), where u = r i • m -1. Here τ = [m > √ 2
] is introduced to avoid the catastrophic cancellation that may occur while computing e • log(2) + log(m) when e = -1 and m ≈ 2 as in [START_REF] Revy | Automated design of floating-point logarithm functions on integer processors[END_REF]. Note that this rewriting enables to avoid code branches, and to take advantage of the rcp instruction, returning an approximation of the reciprocal and available on recent vector micro-architectures.

Furthermore when x is a subnormal floating-point number, a usual way consists in rescaling x in the normal range. Hence let x be a floating-point number defined as:

x = x • 2 λ with λ ∈ {0, • • • , p -1}, ( 2
)
where λ is the number of leading zeros in the binary representation of m. It follows that x = 2 e • m ≥ 2 emin and

log(x) = (e + τ -λ) • log(2) -log(2 τ • r i ) + log(1 + u), (3)
with m ∈ [1, 2 -2 1-p ],
represented exactly with only one precision-p floating-point number. The quantity r i is thus an approximation of 1/m obtained with the rcp instruction truncated on i bits, and u = r i • m -1. In addition the value -log(2 τ • r i ) is retrieved from a lookup table using the i > 0 most significant fraction bits of the significand of r i . Note the way subnormal inputs are handled here. It is different from the one usually seen in vector function implementations, for which calling a fallback routine for these inputs is often preferred. However, treating subnormal inputs in the general flow has no great impact on performance as discussed in Section VI-D.

B. How to determine r i , u, τ , and -log(2 τ • r i )?

Let us now detail the computation of r i , u, τ , and -log(2 τ •r i ).

Determination of r i . Let assume an rcp instruction available, that returns an approximation of 1/m with a relative error bounded by κ:

rcp(m ) = 1/m • (1 + rcp ), with | rcp | ≤ κ. (4)
The quantity r i is then computed by truncating rcp(m ) on i > 0 fraction bits. In order to improve and to re-center error enclosure, we first add one half-ulp on i bits. Thus we have

r i = rcp(m ) • (1 + ri ), with | ri | ≤ 2 -i-1 , since m , r i > 0. It follows that r i = 1/m •(1+ rcp )•(1+ ri ).
Hence:

r i = 1/m • (1 + ) with = rcp + ri + rcp • ri , (5)
where

| | ≤ κ + 2 -i-1 + κ • 2 -i-1 . (6)
This assumption on rcp is not restrictive at all. This kind of instructions is available on a many vector micro-architectures, like on Intel's ones since SSE and considered in this article, or on ARM's ones through the VRECPE (Vector Reciprocal Estimate) instruction. 6 If it is not the case, it can be emulated with full division or lookup table, as discussed in Section V. For example, using Intel intrinsics guide, 7 we know that SSE/AVX instruction sets provide a binary32 rcp instruction with a maximum relative error less than 1.5•2 -12 , that is, with

| rcp | ≤ 1.5 • 2 -12
, while on more recent micro-architectures, AVX-512 instruction set provides also a binary64 rcp instruction with | rcp | ≤ 2 -14 . On SSE/AVX-supporting microarchitectures, the binary32 instruction can also be used to approximate 1/m for a binary64 input: m must first be casted to the binary32 format, thus slightly increasing the approximation error:

| rcp | ≤ 1.5 • 2 -12 + 2 -24 + 1.5 • 2 -36 .
Determination of u. Since u = r i •m -1, we deduce from (5) that the quantity u = , that is, u is the error occurring when approximating 1/m by r i .

Property 1: Let x = 1 be a positive floating-point number as in [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF]. The quantity u in (3) can be represented exactly with only one precision-p floating-point number as long as the parameter i is chosen so that:

i ≤ e max -2 and κ • (2 i + 2 -1 ) < 2 -1 . (7)
Proof: Let M and R be the integer significands of m and r i , respectively. Since r i is obtained after truncating rcp(m ) on i > 0 fraction bits, its significand can be stored on at most i + 1 bits. Hence

m = M • 2 1-p and r i = R • 2 δ-i , Table I SUITABLE VALUES FOR i FOR DIFFERENT COMBINATIONS (eMAX, κ). emax κ i ∈ • • • 127 1.5 • 2 -12 [1, 10] 1023 2 -14 [1, 12] 1023 1.5 • 2 -12 + 2 -24 + 1.5 • 2 -36 [1, 10]
with δ ∈ {-1, 0}. Using u = r i • m -1, we rewrite u as:

u = R • M • 2 δ-i+1-p -1 = 2 δ-i+1-p • R • M -2 p-1+i-δ .
Since p ≥ 2 and i > 0, then if i ≤ e max -2, we deduce that

δ -i + 1 -p ∈ [e min + 1 -p, 1 -p]
is the exponent of a floating-point number. Moreover the value U defined as:

U = R • M -2 p+i-1-δ
is an integer. Using [START_REF] Tang | Table-driven implementation of the logarithm function in IEEE floating-point arithmetic[END_REF], since u = , we deduce that:

|U | ≤ κ • 2 p • (2 i + 2 -1 ) + 2 p-1 . Therefore if i is chosen so that κ • (2 i + 2 -1 ) < 2 -1 , then |U | < 2 p
, and u is a precision-p floating-point number, which concludes the proof.

From Property 1, we know that for certain values of i, the quantity u fits exactly in only one floating-point number in precision p. These values are given in Table I, for different combinations (e max , κ). As a consequence, in order to implement the sequence computing u, a first choice may be to use an fma, returning exactly u in a single instruction, and available on recent micro-architectures. If it is not available, an alternative is to use multi-word arithmetic, like double-double arithmetic in [START_REF] Lauter | Arrondi correct de fonctions mathématiques -fonctions univariées et bivariées, certification et automatisation[END_REF], [START_REF] De Dinechin | Fast and correctly rounded logarithms in double-precision[END_REF]. This is discussed in Section VI.

Determination of τ and -log(2 τ •r i ). The value -log(2 τ •r i ) is an approximation of -log(2 τ /m ), and it is retrieved from a lookup table, indexed only by the i most significant fraction bits of r i , and not by τ . It follows that this table is built such that for the entries corresponding to τ = 0, it contains approximations of -log(1/m ), while for the others (τ = 1), it contains approximations of -log(2/m ).

As a consequence, each cell in the table must represent m values for which τ is the same. In this sense, τ cannot be decided using m , but it must be decided using the i most significant fraction bits of r i , as well, denoted by idx below. And since r i is a floating-point number, τ is computed as:

τ = 1 if idx ≤ (2/ √ 2 -1) • 2 i , 0 otherwise. ( 8
)
This way, when m ≈ 1, we may have idx = 0 and τ = 1, while τ should be 0. Hence in this particular case, the value of τ must be explicitly set to 0. Note that τ in (8) may not correspond exactly to the definition of τ in Section II-A: this is not an issue since we have a certain latitude for the value involved in the test. Usually √ 2 is chosen to re-center the fraction around 1, but other values around √ 2 do also work.

C. How to ensure faithful rounding?

Let r be an approximation of log(x), with x = 1 a positive floating-point number. For accuracy purpose, the value r is represented as the unevaluated sum of two floating-point numbers r h and r in precision p, that is

r = r h + r with |r | < 1 2 ulp (r h ) ≤ 2 -p • |r h |. ( 9
)
Our goal is to compute r so that r h is a "faithful rounding" of log(x), that is, either RD (log(x)) or RU (log(x)). To do this, following [START_REF] Muller | On the definition of ulp(x)[END_REF], a sufficient condition is:

|r h -log(x)| < ulp (log(x)) . (10)
Using the triangular inequality, we have:

|r h -log(x)| ≤ |(r h + r ) -log(x)| + |r |.
Hence from ( 9) and assuming ulp (log(x)) = ulp (r h ), we deduce that if

|(r h + r ) -log(x)| ≤ 1 2 ulp (r h ) , (11)
then [START_REF]Fast hardware-based algorithms for elementary function computations using rectangular multipliers[END_REF]

holds. When log(x) ∈ ]2 k-1 • (2 -2 -p ), 2 k [, with k ∈ Z, we have ulp (log(x)) = 1
2 ulp (r h ): these cases are ignored at implementation time, and the accuracy of the result for these particular inputs is checked a posteriori by exhaustive testing. Now let us distinguish three cases:

• Case 1: When e + τ -λ = 0, and x ∈ X with • Case 3: Otherwise, e + τ -λ = 0: thus we have

X = 1 + κ 1 + 2 -i-1 -2 1-p , 1 -κ 1 -2 -i-2 , we show that log(2 τ • r i ) = 0,
| log(x)| > log( √ 2)
. Hence without loss of generality, we deduce that

| log(x)| > 2 -2 , then |r h | ≥ 2 -2 and ulp (r h ) ≥ 2 -p-1 . It follows that if |r -log(x)| ≤ θ with θ =      2 -2p in Case 1, ulp (min) /2 in Case 2, 2 -p-2 in Case 3, (12)
then [START_REF]Fast evaluation of the elementary functions in single precision[END_REF] holds, and r h is a "faithful rounding" of log(x). For example, for (p, i) = (24, 8), we have: min ≈ 6.1 • 10 -4 and ulp (min) = 2 -34 .

Distinguishing these three cases will help us in certifying the error entailed by the evaluation of the program in finite precision as explained in Section IV-B.

III. GUIDELINES TO BUILD ACCURATE ENOUGH PROGRAM

This section gives some guidelines to build a program that returns a "faithful rounding" of log(x). It is based on the evaluation of a particular polynomial. Then, this section presents how to constraint some coefficients of this polynomial, so that the special input 1 can also be handled in the general flow.

A. Program to handle the general flow

Let x = 1 be a floating-point number. Our goal is now to build a program that computes r = r h +r as in [START_REF] Schulte | Exact rounding of certain elementary functions[END_REF]. Following [START_REF] Revy | Automated design of floating-point logarithm functions on integer processors[END_REF], we use a 3-step process:

• We consider log(x) as the exact result of the function F defined as in (3):

F (x) = (e + τ -λ) • log(2) -log(2 τ • r i ) + log(1 + u).
• Since F cannot be evaluated directly, we approximate the function F by another function P :

P (x) = (e + τ -λ) • L 2 + log_tbl(i) + a(u).
Here L 2 approximates log(2) so as |L 2 -log(2)| ≤ θ 1 , and log_tbl(i) returns -log(2 τ • r i ) with an error no greater than θ 2 . And finally a(u) is a polynomial approximant of the function log(1 + u) over the interval I defined in ( 6) for , and with the approximation error defined as:

max u∈I |log(1 + u) -a(u)| ≤ θ 3 .
• We finally evaluate P by a finite-precision evaluation program P, that computes r h + r and returns r h . Using the triangular inequality, we have:

|r -log(x)| ≤ |F (x) -P (x)| + |P (x) -P(x)| ≤ |e + τ -λ| • θ 1 + θ 2 + θ 3 + θ 4 ,
where θ 4 is a bound on the evaluation error of P. Assuming the three cases of Section II-C, in order to ensure that the condition in (12) holds, we must build a program so that

|e + τ -λ| • θ 1 + θ 2 + θ 3 + θ 4 ≤ θ. (13)
This is detailed in Section IV-B.

B. How to handle the special input 1 in the general flow?

Let us now consider the case x = 1. Hence we have x = 1, and m = 1. And from (4), we know that

rcp(m ) ∈ [1 -κ, 1 + κ].
Recall that the quantity r i equals rcp(m ) plus one halfulp on i bits, then truncated on i fraction bits. Therefore if the parameter i is chosen such that

2 -i-2 ≥ κ and 2 -i-1 ≤ 2 -i -2 1-p -κ, (14)
then, we can deduce that r i = 1. In this case, the computation of log(x) is reduced to the evaluation of log(1 + u). For our implementations, this appears for all the values i in Table I, but for i = 10 when (p, κ) = (24, 1.5

• 2 -12 ) and (p, κ) = (53, 1.5 • 2 -12 + 2 -24 + 1.5 • 2 -36 ).
As a consequence, since u is the error entailed by the approximation of 1/m by r i , we deduce also that u = 0. It follows that to handle the special input 1 in the general flow, i must be chosen to satisfy [START_REF] Lauter | Arrondi correct de fonctions mathématiques -fonctions univariées et bivariées, certification et automatisation[END_REF], and the polynomial approximant a(u) in Section III-A must be built so as the first coefficient (a 0 ) is zero, to ensure that a(u) = 0 in this case.

IV. SOME IMPLEMENTATION DETAILS

This section gives some details on our implementation process, for the example of the binary32 format. After some elements on the way used to write branchless programs, it presents how to certify the accuracy of the output programs.

A. Branchless program

For performance purpose, a key point consists in writing branchless programs. First τ is a boolean introduced to avoid the catastrophic cancellation that may occur while reconstructing the result. It is actually an integer ∈ {0, 1} determined in a branch-free fashion using comparison instructions.

Second, we must compute x = x • 2 λ , defined in (2). Once λ is known, this could be performed with a floatingpoint multiplication. However, as observed in [START_REF] De Lassus Saint-Geniès | Performances de schémas d'évaluation polynomiale sur architectures vectorielles[END_REF], multiplying two floating-point numbers may lead to a huge overhead, when one of both operands is a subnormal number which may be the case. An alternative is to work on the bit string encoding x using the standard binary interchange format encoding [3, § 3.4]. Let the following bit string encode a binary32 subnormal floating-point number x: 0 00000000 000000000 λ=9 101010000000000.

Computing x • 2 λ just consists in shifting this bit string by λ + 1 bits to the left. With λ = 9, we obtain: 0 00000001 01010000000000 0000000000 λ+1=10 which encodes x ≥ 2 emin . This sequence of operations relies on arithmetic and logical instructions, but it requires a routine returning λ if x is a subnormal number, and 0 otherwise. The instruction set we target does not embrace any instruction enabling to count the number of leading zeros in a bit string. Various techniques can be implemented in software for computing this, but either they are costly or they use branches [START_REF]Hacker's Delight[END_REF]. In our implementation, we use the following piece of code, presented here for the binary32 arithmetic, but that can be adapted to any other formats. It works on an integer X, that represents the bit string of a floating-point x. 1) If x is a normal number, its exponent field (X >> 23) is non-zero, and mask = 0. Hence the routine returns 0. 2) Otherwise if x is a subnormal number:

1 typedef union { uint32_t i; float f; } cfloat32; 2 3 uint32_t lambda_or_zero(uint32_t X) { 4 cfloat32 Z; 5 Z.i = X | 0x3f800000; // exponent
x = 2 -126 × 0. 000 • • • 000 λ 1 • • • .
Then at Line 5, we have:

Z = 2 0 × 1. 000 • • • 000 λ 1 • • • .
At Line 7, the subtraction is exact due to Sterbenz lemma. Thus as long as 1 -p ≥ e min , we have Z = 2 -λ × 1. • • • .The exponent field of x being zero, we deduce that mask = 2 32 -1 = 1111 • • • 1111, and value = -λ. Hence the routine returns λ. This is a well-known technique, but to our knowledge, this has never been published.

B. Certified evaluation program

Recall that in this section, we target the binary32 arithmetic, that is, with p = 24. For accuracy purpose, we must build an evaluation program, so that the bound in ( 13) is satisfied. And as mentioned in Section II-C, we must distinguish three cases. Indeed, the error bound θ = 2 -48 , which is the tightest in (12), could not have been proven on the whole input interval.

From now we start by determining the polynomial degree and coefficients. To avoid branches for deciding in which cases the inputs fall, we choose to implement the log(x) function with a single degree-d polynomial. Since θ 1 , θ 2 , θ 3 ≥ 0, we know from ( 12) and ( 13) that L 2 and a(u) must be built so as θ 1 , θ 3 ≤ 2 -48 . To reach such θ 1 bound, L 2 must be represented as the unevaluated sum of two binary32 floatingpoint numbers. Remark also that the lookup table index size i influences the approximation interval, and thus the degree d: once the couple (i, d) is chosen, the polynomial approximant together with the certified error bounds θ 1 , θ 2 and θ 3 are computed using the Sollya tool. 8 This way, we obtain θ 1 ≈ 2 -53.33 . In addition Table II shows the polynomial degree d, the error bounds θ 2 and θ 3 , as well as the memory occupancy in bytes for lookup table and polynomial coefficients, for different table index sizes i. In the rest of the example, we choose to implement log(x) using a table indexed by 8 bits and a degree-4 polynomial, which represents the best trade-off in terms of performance as shown in Section VI.

The remaining part is to compute the program evaluation error bound θ 4 , and to check if the bound (13) holds. This error corresponds mainly to the error due to the evaluation of the polynomial approximant a(u) in finite precision arithmetic, combined with the error entailed by the reconstruction. Note that the way used to evaluate a(u) may be determined by using the software tool CGPE, 9 which enables to build several schemes to evaluate a given polynomial. Then bounding the evaluation error is carried out using the Gappa tool, that uses interval arithmetic and rewriting rules in order to bound error entailed by the evaluation of a program in finite precision. 10 Notice that our implementation works with doublesingle arithmetic, that is, with numbers represented as the unevaluated sum of two binary32 floating-point numbers. Thus we first need to adapt the error bounds computed in [22, § 4] for double-double addition and multiplication to our context, and to pass them to Gappa: we found a relative error bound of 2 -44 and 2 -45 for double-single addition and multiplication, respectively. Using them, we thus find that the absolute evaluation error of the polynomial is not greater than ≈ 2 -51.75 , which is actually less than 2 -48 . Then according to the case, we obtain:

θ 4 ≤      2 -50.75 in Case 1, 2 -45.49 in Case 2, 2 -36.21 in Case 3,
and in all cases, the condition in ( 13) holds. 11V. TOWARD AUTOMATED IMPLEMENTATION This section presents some insights on how to automate the implementation of this function. It starts with a description of the MetaLibm framework, before some details on its extension to our context.

A. The MetaLibm framework

The MetaLibm framework enables to describe the implementation of a function using a meta-language in a Python syntax. More particularly it consists in describing a mathematical function implementation as a set of high-level basic blocks. Then the MetaLibm backend translates each of these basic blocks in a specific instruction set, thus automatically outputting programs optimized for various floating-point formats and micro-architectures [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF], [START_REF] Brunie | Code generators for mathematical functions[END_REF]. For example, the following piece of code describes the multiplication x • x in this metalanguage, where x is the input variable of our program. As other examples, when fma or rcp instructions are not available, MetaLibm may replace them by a call to fma/fmaf libm functions or a loookup table, respectively. So far this framework has been actively developed as a fast and efficient code generation tool for two kinds of final users. The first kind is experienced software developers who aim at implementing elementary mathematical functions libraries optimized for different architectures. The second kind is nonexpert numerical software users, interested in customized mathematical code generation for their applications. Among the main features, it enables:

• To develop multi-architecture-specific meta-codes, and to support transparently different accuracies (improved, normal, degraded or customized), producing automatically different code versions, • To factor common code across software and hardware architectures, for various IO precisions and microarchitectures, reducing thus the development time, • To vectorize meta-code according to different microarchitectures, to fit at best the underlying architecture. The latter mainly consists in removing code branches, and inserting tests and data selections, instead. This technique implies speculative execution, and thus increases the code size, and eventually the evaluation latency. However for throughput purpose, this may be quite efficient as shown on Cephes library [START_REF] Piparo | The VDT Mathematical Library[END_REF]. Our programs are already written in a branchless fashion: thus this vectorization step will not have any impact on the produced optimized code. But this remark reinforces our choice to write branchless codes, including the handling of subnormal inputs in the general flow.

In this work, an effort has been made to enhance the backend of MetaLibm, with all the needed instructions unavailable so far, but also to provide support to build efficient polynomial evaluation schemes, as detailed in Section V-B below.

B. The case of polynomial evaluation schemes

As shown in Section II-A, our implementations of log(x) function rely on the evaluation of a polynomial, whose degree may vary from 4 up to 9 (Table II). Various schemes may be used to evaluate a given degree-d polynomial. In order to achieve high throughputs, we chose to use CGPE (Code Generation for Polynomial Evaluation) to generate an efficient polynomial evaluation scheme on the targeted architecture. Given the polynomial coefficients and a bound on the evaluation error, it enables to automatically write and certify programs to evaluate this polynomial. At first developed for VLIW integer processors to provide fixed-point computation abilities to this kind of hardware, it has recently been extended to handle floating-point computations. But above all, one of its main features is its capability in computing polynomial evaluation schemes, exposing more or less instruction-level parallelism. To do this, it is based on exhaustive and heuristic algorithms that look for low latency polynomial evaluation schemes on abstract customizable architectures.

CGPE is written in C++ and is available as a command-line tool. Since MetaLibm is written in Python, a first step was to develop Python bindings for CGPE (called PythonCGPE). This package provides a non-exhaustive interface to CGPE features in Python that gives means to automatically generate efficient polynomial schemes given adder and multiplier latencies. It uses heuristics for polynomials of degree greater than 6 to avoid blocking MetaLibm whenever the search space grows too large. Note that even with these heuristics, the search space may get huge when the degree increases. Hence as soon as d ≥ 12, we choose to skip the schemes computation step, and to rely only on the classic Estrin's rule. This is a well-known limit of the CGPE software tool, and solving this issue is out of the scope of this article.

VI. NUMERICAL RESULTS

We have written various versions of log(x), and evaluate their performance. This section compares our implementations to existing solutions, and it studies the impact of three design parameters: table size, fma/rcp use, and subnormal handling.

A. Experimental protocol

We measured the reciprocal throughput of handmade and generated routines, and evaluate them against the Libmvec routines. All benchmarks were run on an Intel R Xeon R Processor E5-2650v4, which features the AVX2 instruction set extensions. The operating system was Cen-tOS Linux release 7.4.1708, running a Linux kernel version 3.10.0-693.2.2.el7.x86_64. All source code was compiled using GCC 7.2.0, linked against glibc 2.26, both compiled from source. Shared compiler flags included -O3 -mtune=native. For Libmvec benchmarks, we also had to enable -ftree-loop-vectorize -ffast-math and link against GNU libm. The former option is used to vectorize loops, while the latter enables mathematical simplifications, flushes subnormal arguments and results to zero, among other optimizations, and is required to activate the Libmvec. For all benchmarks, to enable AVX2 code generation we simply used the -march=core-avx2 flag, while we used -march=corei7 to restrict GCC to emitting SSE4 code.

To evaluate these performances, we used either automated micro-benchmarks provided by MetaLibm or custom microbenchmarks for the handmade version and Libmvec routines. These custom micro-benchmarks tried to reproduce plausible workloads by using pseudo-random inputs, but also rare workloads by using only subnormal inputs. The benchmarks were designed to warm up the L1 cache before the main measurement, so that cache-miss penalties could be kept to a minimum. Also, we chose to take the minimum measured reciprocal throughput for each micro-benchmark. We claim that this represents the best the CPU can actually compute when there is the least noise perturbing the measures. This can be further justified by the fact that we benchmarked separately pseudo-random input vectors versus subnormals-only vectors, which may yield lower performance. We also benchmarked constants-only vectors, which might have yielded better throughputs for vectorized memory gather.

B. Performance of the log(x) function

In this section, we measured the reciprocal throughput of the generated routines, and compared them to SSE4 and AVX2 reference implementations coded by hand and to SSE4 and AVX2 versions of Libmvec logf/log routines. In this experiment, our implementations use fma, rcp, and tables indexed by i = 8 bits as in Section IV-B, and we use MetaLibm and custom micro-benchmarks, for generated and handmade codes, respectively. To gain in confidence, the binary32 versions have been exhaustively verified and compared with the values returned by the MPFR library. 12 Performance results are given in Table III, in cycle per element (CPE). We observe that using our approach, we achieve a binary32 implementation with a reciprocal throughput of ≈ 18.8 CPE for vector size = 4 and ≈ 10.5 CPE for vector size = 8. In binary64, the measured performances are ≈ 80.2 CPE for vector size = 2 and ≈ 42.2 CPE for vector size = 4, that is, by a factor of ≈ 4 compared to binary32 versions. We observe similar performances for binary32 routines automatically generated with MetaLibm. Obviously, we are far away from the Libmvec implementations, whose reciprocal throughputs vary from 0.3 up to 6.7 in binary32 according to the vector sizes and the input ranges. This might be due to the fact that binary32 Libmvec implementation relies on the evaluation of a small degree polynomial, done using single precision arithmetic only, while in our cases polynomials are evaluated using double-single arithmetic. A direct consequence is the accuracy of the output: indeed our implementations are correct to 1 ulp, while Libmvec provides functions correct to at worse 4 ulp. Furthermore Libmvec has a scalar fallback for special inputs such as NaNs, infinities or zero. Since subnormal numbers are flushed to zero with -ffast-math, the Libmvec routines are much less efficient when dealing with a subnormals-only vector, by a factor of ≈ 5 in binary32 and ranging from 6.8 up to 11.7 in binary64. Conversely as our scheme unifies normal and subnormal handling, performance is not affected by a subnormals-only vector. (The same remarks hold for the binary64 format.) 12 See http://www.mpfr.org/ and [START_REF] Fousse | MPFR: A Multiple-precision Binary Floating-point Library with Correct Rounding[END_REF]. An interesting observation is that the multiple-indices vectorized accesses to tabulated values were not penalizing, although data were accessed at non-contiguous locations. This may be explained by the fact that all tables are smalltypically less than 2 kB for binary32 precision and less than 4 kB for binary64 -therefore fitting easily in the 32-kB L1d cache of the tested processor.

C. Impact of table size and architectural features

The range reduction presented in Section II-A heavily relies on some particular architectural features, like fma and rcp instructions. In this section, we modified our metaimplementation, to produce binary32 routines relying on multiword arithmetic and table lookup, instead, respectively, to observe the impact on performance of using these features. Table IV reports the performance of these routines in CPE, for different table index sizes, from 3 up to 9 bits, and for different vector sizes (typically 4 and 8),

As expected, the performance of the routines increases when the index size and consequently the table size increase. Using this table we may conclude that, for performance purpose, a table indexed by 8 bits seems to be a good choice, with a thoughput of 19.0 CPE.

Note that fma and rcp instructions have a great impact on performance of our generated implementations. Indeed by replacing fma with multi-word arithmetic, we observe a slowdown factor between 1.29 and 2.61 for vector size = 4, and between 1.31 and 2.4 for vector size = 8. Regarding rcp, using table lookup, instead, leads to a slowdown of a factor between 1.41 and 2.7 for vector size = 4, and between 1.39 and 2.53 for vector size = 8. In the latter case, Property 1 cannot be verified: u must be represented as the unevaluated sum of two floating-point numbers, leading to huge overhead in the polynomial evaluation. This reinforces the interest and efficiency of our approach on recent micro-architectures, supporting both fma and rcp instructions.

D. Impact of subnormal handling

In this section, we measured the cost of supporting subnormals in the main flow. Thus we modified our meta-implementation, to generate binary32 routines with and without subnormal handling. To unplug subnormal handling, it suffices to remove λ from all the computations. Table V shows the performance of these routines in CPE, for a parameter i ranging from 3 up to 9, for different vector sizes (typically 4 and 8).

An interesting observation is that treating subnormal inputs in the main flow has no great impact on the performance of the function. Indeed, for vector size = 4, the overhead remains no greater than 2.3 CPE, expect for i = 5 (22.2 CPE), and it is always less than 2 CPE but for i = 4, which is acceptable. And for vector size = 8, this overhead is always less than 1 CPE, but for i = 3 (1.6 CPE). Remark that for i ∈ {8, 9}, this technique degrade routine performance, compared to the original. This reinforces our choice to avoid branches and fallback routines to treat subnormals separately.

VII. CONCLUSION AND PERSPECTIVES

This article addresses the implementation of log(x) functions on vector micro-architectures, with a particular focus on its automation through the MetaLibm framework. This enables to achieve high throughput implementations with 1-ulp accuracy, optimized for AVX/SSE and AVX2 micro-architectures, with a relatively small overhead due to handling subnormals in the main flow (less than 2 CPE).

In Table III, "n/a" indicates that MetaLibm were not able to generate these functions. This is actually due to the lack of certain meta-instructions in the framework backend. Ongoing work focuses on integrating these requirements. In a near future, this will enable generation of binary64 implementations, as well. More generally, the IEEE 754 standard defines two other formats, namely, the binary16 and binary128. These techniques being parametrized by the precision, as a consequence, using these to write implementations for these formats would be reachable and of interest. This would require efforts to eventually emulate underlying required missing arithmetic (SIMD support for binary128, for example), and to adapt the MetaLibm backend in consequence.

Research direction is twofold: First we could extend our efforts to the implementation of certain elementary functions.

Following [START_REF] Revy | Automated design of floating-point logarithm functions on integer processors[END_REF], a direct extension would be the optimized implementations of log 2 (x) and log 10 (x). In addition, we could concentrate on the design of efficient exponential or trigonometric functions, which are also widely used in HPC. Second we could adapt our work to produce programs with different accuracy levels, like in [START_REF] Lauter | A new open-source SIMD vector libm fully implemented with high-level scalar C[END_REF] where 8 ulp is guaranteed.

  and log(x) = log(1 + u). In this case, since | log(x)| > 2 -p , we have |r h | ≥ 2 -p and ulp (r h ) ≥ 2 -2p+1 . (See [13, Prop. 3] for details.) • Case 2: When e + τ -λ = 0 and x / ∈ X : | log(x)| > min with min = min(| log(inf(X ))|, | log(sup(X ))|). Then |r h | ≥ min and ulp (r h ) ≥ ulp (min) .

  y = Multiplication(x, x, tag="y", precision = ML_Binary32)Then, if scalar architectures are targeted, MetaLibm simply generates the sequence y = x * x, while if x86 architectures with AVX2 are targeted, it provides the following piece of C code, where vec_x is the vector equivalent of x.carg = GET_VEC_FIELD_ADDR(vec_x); tmp = _mm256_load_ps(carg); y = _mm256_mul_ps(tmp, tmp);

Table II IMPLEMENTATION

 II PARAMETERS FOR VARIOUS VALUES i AND p = 24, AND TABLE AND POLYNOMIAL COEFFICIENT SIZES (# BYTES).

	i	d	θ 2	θ 3	table	coefficient
	3	9	2 -52.69	2 -52.01	64	68
	4	7	2 -52.29	2 -49.62	128	52
	5	6	2 -52.10	2 -49.82	256	44
	6	5	2 -51.36	2 -48.88	512	36
	7	5	2 -51.20	2 -54.50	1024	36
	8	4	2 -51.04	2 -48.98	2048	28
	9	4	2 -51.03	2 -52.93	4096	28
						of the
						// floating-point number 1	6
						Z.f -= 1.f;	7
						uint32_t mask = 0xffffffff + ((X >> 23) != 0);	8
						uint32_t value = ((Z.i >> 23) -127) & mask;	9
						return -value;	10
						}	11

Table III MEASURED

 III PERFORMANCES OF OUR IMPLEMENTATIONS COMPARED TO THE GNU LIBMVEC, IN CYCLES PER ELEMENT (CPE).

	Version	ISA	-ffast-math	Accuracy	Workload	Reciprocal throughput (CPE) Reciprocal throughput (CPE)
			enabled			in binary32	in binary64
	Libmvec 128 bits	SSE4	yes	4 ulp	Pseudo-random	1.2	5.9
	Libmvec 256 bits	AVX2	yes	4 ulp	Pseudo-random	0.3	1.3
	Libmvec 128 bits	SSE4	yes	4 ulp	Subnormals	6.7	40.5
	Libmvec 256 bits	AVX2	yes	4 ulp	Subnormals	2.6	15.2
	Handmade 128 bits	AVX2	no	1 ulp	Any	18.8	80.2
	Handmade 256 bits	AVX2	no	1 ulp	Any	10.5	42.2
	ML-generated 128 bits	AVX2	no	1 ulp	Any	19.0	n/a
	ML-generated 256 bits	AVX2	no	1 ulp	Any	11.0	n/a

Table V MEASURED

 V PERFORMANCES (CPE) ON AVX2 OF OUR GENERATED IMPLEMENTATIONS, IN BINARY32 ARITHMETIC AND FOR VECTORS OF SIZE 4 (V4) AND 8 (V8), WITH OR WITHOUT SUBNORMAL HANDLING.

	i	3	4	5	6	7	8	9
	v4 / no-sub	56.3	47.3	21.2	19.0	19.0	17.8 17.8
	v8 / no-sub	28.7	24.7	12.9	11.5	11.5	16.6 16.7

See https://software.intel.com/en-us/node/524289.

See https://sourceware.org/glibc/wiki/libmvec.

See https://developer.amd.com/amd-cpu-libraries/amd-math-library-libm/.

See http://www.yeppp.info/.

See http://sleef.org.

See ARM Compiler armasm Reference Guide.

See https://software.intel.com/sites/landingpage/IntrinsicsGuide/.

See http://sollya.gforge.inria.fr/ and[START_REF] Chevillard | Sollya: An environment for the development of numerical codes[END_REF].

See http://cgpe.gforge.inria.fr/ and[START_REF] Mouilleron | Automatic Generation of Fast and Certified Code for Polynomial Evaluation[END_REF].

See http://gappa.gforge.inria.fr/ and[START_REF] Melquiond | De l'arithmétique d'intervalles à la certification de programmes[END_REF].

The Sollya and Gappa scripts are available upon request.


Valerie Pommier-Budinger
email: valerie.budinger@isae.fr

Olivier Cherrier

Baffle Silencer with Tunable Resonators for Adaptive Control of Variable Tonal Noise

Keywords: Baffle silencer, tunable resonators, neck geometry, partitioned cavities, acoustic impedance, adaptive noise control

The article presents a baffle silencer with tunable resonators consisting of two superimposed and identically perforated plates associated with a partitioned cavity made of thermoplastic resin. One plate is fixed while the other is movable. Displacement of the mobile plate changes the internal shapes of the resonator necks and shifts the resonance frequency of the system to lower values. The contributions of this paper are firstly, the modeling of a panel with tunable resonators made of necks with a variable geometry and a partitioned cavity in resin and, secondly, the use of the model to elaborate a control strategy to attenuate variable tonal noise. All of the theoretical studies are validated by experimental measurements. Final results show the efficiency of the silencer in attenuating a tonal noise that varies between 2000 and 2800 Hz.

INTRODUCTION

The article focuses on silencers with acoustic resonators made of perforated sheets bonded onto honeycomb cavities for noise reduction. The advantage of these resonators is that they can be used in harsh conditions. Their main drawback is that they are selective in frequency. They are thus efficient for tonal noise reduction. In the case of variable tonal noise, the efficiency is limited if the resonators are not tunable. Many studies have been performed on adaptive resonators and have led to patents: [START_REF] Kotera | Resonator type silencer[END_REF], [START_REF] Mclean | Variably tuned Helmholtz resonator with linear response controller[END_REF], [START_REF] Cheng | Tunable acoustic system[END_REF], [START_REF] Paschereit | Apparatus for damping acoustic vibrations in a combustor[END_REF], [START_REF] Stuart | Variable resonator[END_REF], [START_REF] Kotsun | Helmholtz resonator[END_REF], [START_REF] Kudernatsch | Exhaust gas system with Helmholtz resonator[END_REF], [START_REF] Ciray | Exhaust processor with variable tuning system[END_REF]. The reading of these patents enables identifying two types of solutions to obtain adaptive resonators: either by the control of the volume of the cavities, or by the control of the neck (perforated holes) geometry.

In order to control the volumes of the cavities, [START_REF] Konishi | Tunable acoustic absorbing system using a deep hole array[END_REF] proposed tunable acoustic resonators whose depth is controlled by C-shape micro actuators. De [START_REF] Bedout | Adaptive-passive noise control with self-tuning helmholtz resonators[END_REF] used a rotating internal radial wall inside the resonator cavity to vary the volumes of the cavities. [START_REF] Kostek | Hybrid noise control in ducts[END_REF] studied the robust feedback control of such systems. [START_REF] Kobayashi | Newly-developed adaptive noise absorption control technology for high speed fan noise reduction[END_REF] successfully implemented resonators with tunable cavities in a turbofan.

With respect to the control of the neck geometry and the study of its influence on the resonator acoustic behavior, [START_REF] Birnbach | Increasing of sound absorption proprieties of acoustical panels based on Helmholtz resonators of double plates with orifice inlet neck by geometrical configuration variations[END_REF] studied a resonator with two perforated plates and an inlet air gap. They demonstrated that for a very small distance between the plates from 0.05 mm to 0.1 mm, the resonators have high acoustic impedance and low absorption and that, by increasing the distance, the impedance decreases and the absorption becomes maximum. [START_REF] Chanaud | Effects of geometry on the resonance frequency of Helmholtz resonators[END_REF] studied the radiation impedance for geometries of non-circular orifices. [START_REF] Tang | On Helmholtz resonators with tapered neck[END_REF] studied a resonator with a tapered neck. Results show that the resonance frequency increases with the tapering length and that absorption increases with the tapered neck slope. Esteve et al. (2004) proposed to tune the resonator by using an accordion neck to vary the length of the neck or by using an iris diaphragm to vary the resonator opening.

As changing the volumes of cavities requires heavy machinery that may be unsuitable for some applications (e.g. aerospace applications), this article proposes to study adaptive resonators whose neck geometry can be controlled. An easy way to produce such resonators is to associate two superimposed and identically perforated plates with air cavities. One plate is fixed and bonded to the cavities and the other plate is movable. This concept enables changing the internal shapes of the holes of the perforated layers and the resonance frequency of the system. These tunable resonators have already been studied by [START_REF] Powell | NASA's subsonic jet transport noise reduction research[END_REF]. They showed that the resonance frequency shifts to lower values but they failed to develop a model. A circular silencer with adaptive resonators based on the same principle (two inner cylinders, a fixed one and a movable one, associated with a cavity) has been developed by [START_REF] Nagaya | Silencer consisting of a two stage Helmholtz resonator with auto-tuning control[END_REF] in order to reduce a blower noise with two varying peaks. He proved the efficiency of the system but did not develop a model for the silencer. In [START_REF] Cherrier | Panel of resonators with variable resonance frequency for noise control[END_REF], we developed a model of such tunable resonators with two perforated plates.

We now have built a baffle silencer with panels of tunable resonators made with two plates and this article proposes the model and the control of this baffle silencer.

Part 2 of the article develops the geometry of the tunable acoustic resonators and describes the baffle silencer. Part 3 deals with the modeling of a silencer panel made of an array of tunable resonators. The acoustic impedance at normal incidence of a tunable resonator panel is computed, as well as the acoustic absorption coefficient. All computation results are validated by experimental results in an impedance tube. Finally, in part 4, the performance of the silencer is evaluated and the control of the silencer is demonstrated. The capability of the silencer to track and attenuate a tonal noise that varies according to a frequency ramp is measured.

BAFFLE SILENCER WITH TUNABLE RESONATORS

Description of tunable resonators

Helmholtz resonators are characterized by a high absorption coefficient but are selective in frequency: such resonators are efficient around a resonance frequency that depends on the size, shape and porosity of the necks and on the thickness of cavities. In order to vary the resonance frequency, we propose to use tunable resonators whose necks are generated by two perforated plates, side by side, one fixed and one mobile (Figure 1). The mobile plate can translate a distance that can vary from zero to the hole perforation diameter (2r). These resonators have been studied in [START_REF] Cherrier | Panel of resonators with variable resonance frequency for noise control[END_REF].

Description of the baffle silencer

The silencer is made of three baffles with tunable resonators on each side (Figure 2 and Figure 3).

There are so six noise absorbent panels. The mobile parts of the three baffles are actuated by a single electromechanical linear actuator (Part B). The link between the actuator and the mobile parts of baffles is rigid, so the displacement of the actuator is the same as that of the baffles. A loudspeaker (Part A) is connected to one end of the silencer to generate noise to test the silencer. estimation of the resonance frequencies is accurate but the thicknesses and the amplitudes of the resonance peaks are either much lower or greater than the measured ones. The differences can be explained by two facts:

1 -Compared to the resonators studied in [START_REF] Cherrier | Panel of resonators with variable resonance frequency for noise control[END_REF], the cavity of the silencer panel is partitioned in several small cavities and, because of the manufacturing method, the thickness of the cavities wall is not negligible (1mm).

2 -The silencer cavities are made of ABS resin, and therefore the acoustic impedance of the bottom of the cavity cannot be regarded as infinite, as opposed to the bottom of a cavity in metal.

Moreover, the surface of the cavities is rough because of the manufacturing process (3D printing) and there may be traces of glue in the cavities due to the bonding of the fixed perforated plate onto the cavities. This increases the acoustic resistance of the silencer panel.

Up-dating the model

The model proposed in [START_REF] Cherrier | Panel of resonators with variable resonance frequency for noise control[END_REF] can be improved to take into account the repartition of the resonators, the wall thickness and the influence of the cavities material.

To model the repartition of the resonators and the wall thickness, it is necessary to define another open area ratio as proposed in [START_REF] Allard | Propagation of sound in porous media: modeling sound absorbing materials[END_REF]. It corresponds to the open area ratio of the plate:

(1)

To model the losses in the ABS material and on its surface, an acoustic normalized resistance at the interface of the cavities and of the perforated plate 1 (Figure 1) must be added to the model. It has been measured on a sample in ABS material in an impedance tube and its mean value on the frequency range of the model is:

11 loss R  (2)
Finally, the normalized impedance at normal incidence of the panel with tunable resonators in plastic resin can be expressed as:

panel panel panel Z R jX  (3) with normalized acoustic resistance R panel : 12 ll 1 22 2 1 2 1 2 1 4 EE hole E E E panel loss s s s plate E E hole S PS R R R R R r S S S r                                            (4)
and normalized acoustic reactance X panel :

12 1 cot( ) l l 22 hole EE panel cavity i E e plate E S X kL k k k S                             (5)
that are both functions of the geometry and of mobile plate translation D.

In equation ( 5), interior mass end correction i  equals to    were defined in Cherrier et al. (2012) as the perimeter, the surface and the length of an elliptic piston created by the overlapping of the two perforated plates and representing an additional mass of air in motion near the orifice lips of the tunable resonator.

Model validation

The new model is validated by comparing simulation results with experimental results obtained in an impedance tube with a sound level of 110dB. Figure 6 to Figure 8 show respectively the absorption coefficient, the normalized reactance and the normalized resistance for different displacements D in mm of the mobile plate. Figure 6 shows that the amplitudes of the resonance peaks and their thicknesses are now closer to the measurements results. Figure 7 shows that the estimation of the resonance frequency of the panel (frequency for which one reactance is null) is very accurate.

Table 1 gives the percentage error between computed and measured resonance frequencies. The estimation of the resistive effects that come from the viscous interaction between acoustic flow and structure (Figure 8) is less accurate. The behavior is not well-characterized in low frequencies. Modeling these phenomena is more difficult.

CONTROL OF THE BAFFLE SILENCER

Silencer attenuation

Before tackling the issue of control, the attenuation of the silencer is measured on the frequency range 1800-2800Hz with a step of 100Hz. An input tonal noise that consists of a single sine wave at fixed frequency is generated by the loud speaker and then the acoustic sound power of the source is measured with and without the silencer. The acoustic attenuation is the difference between these two acoustic sound powers.

For each frequency, the mobile plates of all panels are moved of the same displacement in order to get the best attenuation. The experimental results are given in Figure 9 and correspond to the legend "silencer in active mode".

The acoustic attenuation is also measured when there is no displacement of the mobile plates. The silencer is then said to be in "passive mode".

For frequencies around 2800 Hz, the silencer has the same behavior in active and passive mode since the displacement of the mobile plates is null. The maximum gain is obtained at 2400 Hz: the attenuation is increased by 20 dB with the silencer in active mode compared with the silencer in passive mode. Below 1800 Hz, the attenuation produced by the silencer is weak or null, even in active mode, since the maximum displacement of the mobile plates is reached. The electromechanical actuator that enables the displacement of the mobile plates is controlled in position with a proportional-integral controller in order to ensure good accuracy. The displacement of the actuator rod is supposed to be equal to that of the mobile plates since the connection between the mobile plates and the actuator is rigid.

The mobile plates displacement required to best attenuate a tonal noise is computed from a FFT analysis of the measured noise. The FFT enables the determination of the frequency of the most disturbing tonal noise (tonal noise with the highest magnitude). The displacement required to attenuate the tonal noise is then computed by using a relation between the displacement of the mobile plates and the resonance frequency of the absorbent panels in the silencer. This relation has been extracted from the panel model of section 3.2. Indeed, for frequencies between 1800 and 3000Hz, the displacement for which the attenuation is the highest has been computed and a polynomial expression of the silencer acoustic behavior has been established. The polynomial that gives displacement D (in meter) versus frequency f (in Hz) is:

( ) (6)

Experimental results

Some experimental tests are performed to validate the proposed adaptive noise control system. The silencer is placed in the middle of a room of large volume where reflection from the walls has been reduced with absorbent materials. The distance between the noise source and the floor is 1.2m. The distance between the microphone and the silencer output is 1m. The measurements are performed in the output axis of the silencer.

The silencer is first used to reduce a pure variable tonal noise consisting of a sine wave whose frequency varies as a ramp between 2800Hz and 2000 Hz in 10s. Figure 12 In passive mode, as the displacement of the mobile plates is null, the silencer is configured to attenuate mainly the frequencies around 2800Hz. So the output noise is minimal at the beginning of data from model polynomial the test that starts at 2800 Hz and increases as the frequency of the input tonal noise decreases. In active mode, the displacement of the mobile plates increases as the frequency decreases. The output noise is also increasing but much less than in passive mode. For 2400Hz, the gain between the silencer in passive and active mode is of 20dB, which is consistent with the results of Figure 9. The silencer is then tested where the noise consists of a stepped sine wave whose frequency varies between 2400Hz and 2800 Hz in 8s. Figure 13(a) represents the frequency step of the sine wave. the good position when the frequency of the sine wave varies from 2800Hz to 2400Hz and 0.9s when the frequency of the sine wave varies from 2400Hz to 2800Hz.

CONCLUSION

This article presents a baffle silencer with tunable resonators for the adaptive noise control of varying tonal noise. The tunable resonators consist of two superimposed and identically perforated plates associated with cavities. One plate is mobile and when it is moving, the necks of the resonators are restricted. The resonance frequency then shifts to lower values. This model is used to determine the frequency range at which the silencer will be efficient. It is also used to elaborate a strategy to control the silencer. Results show the efficiency of the silencer to attenuate a tonal noise that varies between 2000 and 2800 Hz.

The proposed adaptive noise control system is easy to implement and different from active noise control systems based on the interference with secondary noise sources. It can be used for machines with variable tonal noises such as fans of aircraft engines and rotating machines controlled at variable speed.

In reference to the model of the tunable resonators baffle, as some of the foreseen applications (e.g. fans controlled at variable speed) are exposed to flow, it would be relevant in future work to take into account the flow over the perforated interface. It could also be useful to build a model that includes high acoustic excitation levels.

With respect to the silencer, improvements can be made to the system by choosing a smaller actuator with a better dynamics to get a more compact solution and to improve the response time of the silencer. The control strategy can also be improved by regularly measuring the behavior of the silencer and by calculating from these measurements the polynomial used in the control law. This update of the polynomial would make the silencer less sensitive to system changes during its use.



  cavity : open area ratio of the cavity  plate : open area ratio of the plate

Figure 1 .

 1 Figure 1. Principle of the tunable resonator with two perforated plates and a single cavity

Figure 2 .Figure 3 .o

 23 Figure 2. Scheme of the silencer

Figure 4 .

 4 Figure 4. Cavities in honeycomb

Figure 5 .

 5 Figure 5. Absorption coefficient at normal incidence for displacements D of the mobile plate (D=0mm, D=0.78mm, D=1.38mm): computed with the model of Cherrier et al. (2012) (---) and measured ( _____ )

  E , S E and E

Figure 6 .Figure 7 .Figure 8 .

 678 Figure 6. Absorption coefficient at normal incidence for displacements D of the mobile plate Curves computed with the new model: D=0mm ( _____ ) , D=0.78mm (----), D=1.38mm (----)

FigureFigure 10 .

 10 Figure 9. Silencer attenuation

Figure 11 .

 11 Figure 11. Acoustic behavior of the silencer

  (a) represents the frequency sweep of the sine wave. Figure 12(b) shows the sound level in dB at the output of the silencer when the silencer is in passive mode (no displacement of the mobile plates) and in active mode. The displacement of the mobile plates measured when the silencer is in active mode is shown in Figure 12(c).

Figure 12 .

 12 Figure 12. Experimental results of the baffle silencer with tunable resonators with a swept sine wave

Figure 13 (

 13 Figure 13(b) gives the sound level in dB at the output of the silencer in active mode and Figure 13(c) the displacement of the mobile plates. This test shows that the control is efficient and allows evaluating the dynamics of the silencer under these conditions. It takes 1.6s for the plates to reach

Figure 13 .

 13 Figure 13. Experimental results of the baffle silencer with tunable resonators with a stepped sine wave

  a tunable resonator panel is proposed and validated by experimental measurements.

Table 1 .

 1 Comparison of computed and measured resonance frequencies