Marc Daumas

Lionel Lacassagne

Daniel Menard

Éric Petit

Alexis Werey

Amine Najahi

Bernard Goossens

Christophe Negre

David Parello

Jean-Marc Robert

Matthieu Martel

Nasrine Damouche

Philippe Langlois

de l'équipe DALI, actuels ou passés, que j'ai pu côtoyer ces dernières années

Keywords:

porté la lourde charge de relire, corriger et valider ce manuscrit dans le cours temps qui leur était imparti. Leurs retours m'ont permis d'améliorer la qualité de ce manuscrit. Un grand merci également aux examinateurs de cette thèse

Extrait du guide à l'usage du candidat au Doctorat en vue de la soutenance de thèse (version du 19 février 2018) de l'école doctorale n°305 E 2 Énergie Environnement de l'Université de Perpignan Via Domitia : "La langue des thèses et mémoires dans les établissements publics et privés d'enseignement est le français, sauf exceptions justifiées par les nécessités de l'enseignement des langues et cultures régionales ou étrangères ou lorsque les enseignants sont des professeurs associés ou invités étrangers. La maîtrise de la langue française fait partie des objectifs fondamentaux de l'enseignement. [. . . ] Lorsque cette langue n'est pas le français, la rédaction est complétée par un résumé substantiel en langue française.

ED 305 : une synthèse d'une vingtaine de pages en français, intégrée au manuscrit, mettant en évidence les apports principaux du travail de thèse." introduction De plus en plus de découvertes scientifiques sont faites grâce à la puissance de calcul qu'offrent des super-ordinateurs comme la grille de calcul du CERN [START_REF] Aad | Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC[END_REF][START_REF] Chatrchyan | Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC[END_REF]. 1 Par exemple, on peut considérer que la découverte du boson de Higgs en 2012 doit sa rapidité, au moins en partie, au traitement informatisé des 13 petaoctets (13 × 10 15 o) de données brutes que produisent les expériences ATLAS et CMS chaque année [START_REF] Brumfiel | High-energy physics: Down the petabyte highway[END_REF]. Ces données nécessitent en effet des calculs complexes de reconstruction des phénomènes physiques avant de pouvoir être exploitées. Cependant, le temps passé par la grille à évaluer des fonctions élémentaires telles que l'exponentielle ou le logarithme peut représenter jusqu'à un quart du temps de calcul total [START_REF] Hauth | Development and Evaluation of Vectorised and Multi-Core Event Reconstruction Algorithms within the CMS Software Framework[END_REF][START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF].

En effet, beaucoup de processeurs généralistes ne disposent pas d'instructions matérielles permettant d'approximer efficacement l'ensemble des fonctions élémentaires recommandées par la norme IEEE 754 pour l'arithmétique flottante [START_REF][END_REF]. Il est alors nécessaire d'utiliser des bibliothèques de fonctions mathématiques (couramment appelées "libm") qui diffèrent souvent en termes de précision ou de performance.

Le premier défi que pose cette thèse est la conception d'implémentations performantes de fonctions élémentaires en arithmétique flottante, et en particulier leur optimisation pour les architectures matérielles proposant des instructions vectorielles SIMD. Mais trouver le bon équilibre entre performances et précision des calculs est un 1 Conseil Européen pour la Recherche Nucléaire. https://home.cern/about vii viii résumé substantiel en français travail de longue haleine. C'est pourquoi un second défi est l'automatisation des procédés d'écriture de bibliothèques mathématiques (aussi appelées libms). Cela permettrait aux utilisateurs de créer des implémentations de fonctions élémentaires adaptées à leurs besoins et optimisées pour leurs architectures matérielles.

En ce sens a été lancé le projet ANR MetaLibm1 en 2014. Ce projet a pour but de développer des outils de génération de code pour les fonctions mathématiques ainsi que pour les filtes numériques [START_REF] Brunie | Code Generators for Mathematical Functions[END_REF]. Il implique les acteurs suivants : Pequan/LIP6 (Sorbonne Université), Inria/SOCRATE (INSA Lyon), le CERN, et DALI/LIRMM (Université de Perpignan Via Domitia). Cette thèse s'inscrit dans le cadre de ce projet. Nos travaux portent sur la génération de codes hautes performances sous contraintes architecturales pour l'implémentation de fonctions mathématiques.

arithmétique flottante

Le Chapitre 1 de cette thèse rappelle l'état de l'art sur lequel se basent nos travaux :

• l'arithmétique flottante et plus particulièrement sa standardisation avec la norme IEEE 754-2008 ;

• l'évaluation des fonctions élémentaires, de façon générale mais avec un accent sur les implémentations logicielles ;

• et les outils actuels de génération pour automatiser en partie ou en totalité l'implémentation de fonctions élémentaires.

Arithmétique flottante IEEE 754

L'arithmétique des ordinateurs est fondamentalement limitée par la quantité de mémoire disponible. Ainsi, seul des ensembles finis de nombres peuvent être représentés par nos machines. Pour approximer les nombres à virgule, de nombreuses représentations existent, mais une des plus communément utilisées est celle définie par la norme IEEE 754-2008 [START_REF][END_REF]. Ce standard définit plusieurs formats, binaires ou décimaux, ainsi que l'arithmétique dite « flottante » pour les formats dédiés au calcul. Cette thèse se concentre sur les formats arithmétiques binaires comme binary32 ou binary64, décrits dans la Table 1.1 (p.10). L'arithmétique de ces formats est différente de l'arithmétique sur les nombres réels, car chaque nombre, opérande ou résultat, a une précision finie. Ainsi, des erreurs d'arrondi peuvent avoir lieu dès lors qu'un opérande ou un résultat n'est pas représentable dans son format de destination. Pour les formats utilisés dans cette thèse, le standard IEEE 754-2008 requiert que les opérations arithmétiques dites « de base » que sont l'addition, la soustraction, multiplication, la division, ix mais aussi la racine carrée et l'opération fusionnée « multiplication puis addition » (FMA), soient correctement arrondies, c'est à dire que le résultat final doit être identique à l'arrondi du résultat infiniment précis.

Cependant, pour d'autres opérations, et notamment des fonctions élémentaires comme le sinus, l'exponentielle ou le logarithme, cette propriété est seulement recommandée. En effet, contrairement aux opérations de base, ces fonctions sont généralement plus difficiles à arrondir correctement, car beaucoup produisent des résultats transcendants.

Un nombre transcendant n'a pas de suite périodique de chiffres dans sa partie fractionnaire, donc il ne peut être ni un nombre flottant de précision finie, ni exactement au milieu de deux nombres flottants de précision finie. Cependant, il peut être extrêmement proche d'une de ces valeurs. Cela pose un dilemme quant à la précision supplémentaire nécessaire pour arrondir correctement ce type de nombres, appelé dilemme du fabricant de tables (DFT), qui peut être beaucoup plus grande que la précision de destination.

Évaluation de fonctions élémentaires

L'évaluation des fonctions élémentaires repose en général sur des compromis entre performance et précision, à cause du DFT. Certaines méthodes visent l'arrondi correct, et doivent donc prévoir des algorithmes plus précis (et donc plus coûteux) pour les cas difficiles à arrondir. Mais il est possible de séparer l'évaluation en deux étapes pour gagner en performance : une première étape dite rapide fournit un résultat correctement arrondi la plupart du temps, mais si l'arrondi peut ne pas être correct, alors une deuxième étape dite précise est exécutée et renvoie le résultat attendu. À l'inverse, d'autres méthodes privilégient la performance à la précision en garantissant des bornes d'erreur plus souples, par exemple bornées par 1 ulp (on parle alors d'arrondi fidèle).

Dans cette thèse on s'intéresse à l'évaluation logicielle de fonctions élémentaires se décomposant en trois grandes étapes : une réduction d'argument, une évaluation polynomiale et une reconstruction. La réduction d'argument consiste à utiliser des propriétés mathématiques comme la périodicité, la parité, ou des identités remarquables afin de réduire le domaine d'évaluation à un domaine restreint. Par exemple, le cosinus étant 2π-périodique, on peut restreindre son évaluation à l'intervalle [0, 2π]. D'autres propriétés permettent de réduire encore cet intervalle à [0, π/4]. Sur un intervalle raisonnablement petit, il est devient intéressant d'approximer la fonction par un polynôme car le matériel est en général efficace pour ce genre de calculs. Cependant, si le degré nécessaire est trop grand pour atteindre la précision voulue, il est souvent possible de le diminuer à l'aide de tables de correspondances. Il s'agit de valeurs précalculées et stockées en mémoire qui permettent de réduire encore un peu plus l'intervalle d'évaluation.

Pour un polynôme quelconque, il y a en général de nombreuses manières de l'évaluer. On appelle ces méthodes des schémas d'évalua-x résumé substantiel en français tion polynomiale. Un des schémas les plus utilisés est celui dit de Horner, mais nous montrons dans le Chapitre 3 que d'autres schémas peuvent être plus efficaces.

Enfin, la reconstruction consiste à combiner les résultats des réductions d'argument mathématique et à base de table avec ceux des approximations polynomiales de façon à reformer le résultat final.

Chacune de ces étapes induit des erreurs d'approximation ou d'évaluation qu'il convient de quantifier pour garantir des bornes d'erreur. L'implémentation de ces étapes ainsi que le calcul de bornes d'erreur pour différentes précisions d'entrées/sorties sont des tâches plutôt laborieuses et sujettes à de nombreuses erreurs humaines [START_REF] Kupriianova | Towards a modern floating-point environment[END_REF]. C'est pourquoi on cherche à automatiser toujours plus les différentes phases de l'implémentation de fonctions, via des générateurs de code.

Génération de code pour les fonctions élémentaires

Sollya, un environnement pour le développement de codes numériques Sollya est un langage de script ainsi qu'une bibliothèque C pouvant effectuer certains calculs symboliques sur de nombreuses fonctions élémentaires, mais aussi calculer différentes approximations polynomiales comme celles de Taylor, de Remez ou encore FP-minimax [START_REF] Sylvain Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF][START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF]. Parmi ses fonctionnalités, Sollya fournit les routines suivantes :

• guessdegree essaye de deviner le degré minimal requis pour approximer une fonction sur un intervalle donné avec une certaine marge d'erreur ;

• taylor, remez et fpminimax peuvent générer des polynômes d'approximation pour les fonctions élémentaires [START_REF] Brisebarre | Efficient polynomial L ∞ -approximations[END_REF][START_REF] Sylvain Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF][START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF] ;

• implementpoly peut générer du code C qui implémente une évaluation polynomiale optimisée en double, double-double ou triple-double, suivant le schéma de Horner [START_REF] De | Optimizing polynomials for floating-point implementation[END_REF].

CGPE, générateur de schémas d'évaluation polynomiale

Bien que Sollya fournisse implementpoly pour évaluer un polynôme, seul le schéma de Horner est considéré. Or, on peut souhaiter atteindre des débits d'opérations plus élevés en utilisant d'autres schémas d'évaluation. C'est ce que propose CGPE, un logiciel en ligne de commande et une bibliothèque C++ pouvant générer du code efficace pour l'évaluation polynomiale sur différentes architectures [START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors -Polynomial evaluation-based algorithms and certified code generation[END_REF]. CGPE peut notamment générer des schémas d'évaluation exposant le plus de parallélisme d'instruction (ILP). CGPE est aussi capable de générer des schémas efficaces satisfaisant une certaine borne d'erreur. Dans ce cas, CGPE peut formellement vérifier cette propriété en générant une preuve Gappa, qui est présenté ci-dessous. L'analyse d'erreur est une tâche pénible et source d'erreurs, mais peut être automatisée grâce à des outils comme Gappa [START_REF] Daumas | Certification of Bounds on Expressions Involving Rounded Operators[END_REF][START_REF] Florent De Dinechin | Assisted Verification of Elementary Functions Using Gappa[END_REF][START_REF] Florent De Dinechin | Certifying the Floating-Point Implementation of an Elementary Function Using Gappa[END_REF]. Gappa peut automatiser et vérifier une analyse d'erreur formalisée. Il gère notamment les formats binary32, binary64 et binary128 [START_REF] Melquiond | De l'arithmétique d'intervalles à la certification de programmes[END_REF].

Le projet ANR MetaLibm

MetaLibm est le nom de multiples projets de génération automatique de code qui ont fusionné pour donner naissance au projet ANR Me-taLibm, au sein duquel cette thèse s'est déroulée. Un de ces projets est Metalibm-lutetia, conçu pour être utilisé par des développeurs avec ou sans expertise en calcul numérique. En effet, il permet de générer du code pour des expressions mathématiques générales en prenant en compte des spécifications fonctionnelles comme la précision ou l'intervalle d'évaluation. Kupriianova et Lauter ont développé un algorithme de découpage en sous-domaines qui permet d'optimiser l'approximation de ces expressions mathématiques par plusieurs polynômes [START_REF] Kupriianova | A domain splitting algorithm for the mathematical functions code generator[END_REF].

Le second projet est Metalibm-lugdunum, plutôt conçu pour les développeurs de bibliothèques mathématiques, mais pouvant aussi être utilisé comme un générateur d'implémentations efficaces pour des développeurs sans expertise numérique. Cet environnement repose sur un méta-langage permettant la description de plusieurs implémentations abstraites et factorisées, partageant des algorithmes communs. Ainsi, ces algorithmes communs peuvent être réutilisés à l'envi afin de gagner du temps de développement. Ce sont finalement les passes d'optimisation et de génération de code qui transcrivent le méta-code en implémentations efficaces pour différentes architectures et spécifications [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF][START_REF] Brunie | Code Generators for Mathematical Functions[END_REF]. Par exemple, le méta-code suivant décrit la multiplication x • x dans le méta-language, où x est la variable d'entrée de la fonction à implémenter. y = Multiplication(x, x, tag="y", precision = self.precision) Si la cible est du code portable C11 en précision binary32, Metalibm génère l'expression y = x * x;, tandis que si la cible est x86_avx2 en précision binary32 avec huit mots par vecteurs, ce sera le code carg = GET _ VEC _ FIELD _ ADDR(vec _ x); tmp = _ mm256 _ load _ ps(carg); y = _ mm256 _ mul _ ps(tmp, tmp); où vec _ x est l'équivalent vectoriel de x.

Le troisième projet est Metalibm-tricassium. Il concerne les filtres numériques pour notamment le traitement du signal. Ceux-ci ont en effet de nombreux points communs avec les polynômes ou les fonctions rationnelles [START_REF] Lopez | Optimal implementation of linear filters in fixed-point arithmetic[END_REF][START_REF] Volkova | Towards reliable implementation of digital filters[END_REF].

xii résumé substantiel en français

Programmation génératrice pour le calcul

Pour le calcul numérique, il existe de plus en plus de langages dédiés ou de compilateurs source à source ciblant différentes optimisations [20, 33, 35, 36, 48, 49, 54, 58-60, 80, 81, 83, 84, 112, 116, 117, 126, 131, 134, 145, 161, 162, 165, 166]. Cette liste est non-exhaustive et pourrait probablement être mise à jour tous les mois. Parmi ces projets, certains visent à améliorer voire à certifier la précision des calculs en arithmétique flottante ou virgule fixe, comme Sardana/Salsa [START_REF] Damouche | Improving the numerical accuracy of programs by automatic transformation[END_REF][START_REF] Ioualalen | Synthesis of arithmetic expressions for the fixed-point arithmetic: The Sardana approach[END_REF][START_REF] Ioualalen | Synthesizing accurate floating-point formulas[END_REF] ou Rosa/Daisy [START_REF] Darulova | Towards a Compiler for Reals[END_REF][START_REF] Darulova | Synthesis of fixed-point programs[END_REF][START_REF] Izycheva | On Sound Relative Error Bounds for Floating-Point Arithmetic[END_REF]. D'autres se focalisent sur de meilleurs compromis entre performances et précision en matériel [START_REF] Istoan | High-performance Coarse Operators for FPGAbased Computing[END_REF][START_REF] Pasca | High-performance floating-point computing on reconfigurable circuits[END_REF][START_REF] Thomas | A General-Purpose Method for Faithfully Rounded Floating-Point Function Approximation in FPGAs[END_REF], comme FloPoCo [START_REF] Florent De Dinechin | Automatic generation of polynomial-based hardware architectures for function evaluation[END_REF][START_REF] Florent De Dinechin | Generating high-performance custom floating-point pipelines[END_REF][START_REF] De | Designing Custom Arithmetic Data Paths with FloPoCo[END_REF], ou en logiciel (PeachPy [START_REF] Dukhan | PeachPy: A Python Framework for Developing High-Performance Assembly Kernels[END_REF][START_REF] Dukhan | PeachPy Meets Opcodes: Direct Machine Code Generation from Python[END_REF][START_REF] Dukhan | Methods for High-Throughput Computation of Elementary Functions[END_REF]), ou les deux [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF].

tables exactes pour les fonctions élémentaires

On a vu que les fonctions élémentaires sont omniprésentes dans de nombreuses applications de calcul haute performance. Cependant, leur évaluation est souvent une approximation. Pour atteindre de grandes précisions, on utilise souvent des propriétés mathématiques, des valeurs tabulées et des approximations polynomiales [START_REF] Defour | Elementary functions : algorithms and efficient implementation for correct rounding for the double precision[END_REF]Ch. 2]. Typiquement, chaque étape de l'évaluation combine des erreurs d'approximation (p. ex. en remplaçant une fonction élémentaire par un polynôme) et des erreurs d'évaluation en précision finie (p. ex. l'évaluation d'un tel polynôme). Les valeurs tabulées ne font pas exception, et contiennent généralement des erreurs d'arrondi dues à la transcendance de la plupart des fonctions élémentaires [START_REF] Markstein | IA-64 and elementary functions: speed and precision[END_REF].

Dans le Chapitre 2, nous présentons une méthode permettant de supprimer ces erreurs d'arrondi, qui peut être intéressante quand au moins deux termes sont tabulés dans chaque entrée de la table. Notre technique transfère toutes les potentielles erreurs d'arrondi dans un seul terme correctif, ce qui permet d'économiser asymptotiquement deux fois plus de bits que les méthodes de l'état de l'art. Pour les fonctions trigonométriques et hyperboliques sin, cos, sinh, cosh, nous montrons que les triplets pythagoriciens permettent la construction de telles tables en un temps raisonnable. Lorsque l'on vise l'arrondi correct en double précision pour ces fonctions, nous montrons aussi que cette méthode pourrait permettre de réduire, pendant la reconstruction, jusqu'à 29% les accès mémoire et jusqu'à 42% le nombre d'opérations flottantes.

Ces résultats ont fait l'objet de deux communications dans des conférences [START_REF] De | Range reduction based on Pythagorean triples for trigonometric function evaluation[END_REF][START_REF] De | Réduction d'argument basée sur les triplets pythagoriciens pour l'évaluation de fonctions trigonométriques[END_REF] et d'une publication dans un journal [START_REF] De | Exact Lookup Tables for the Evaluation of Trigonometric and Hyperbolic Functions[END_REF].

Évaluation à base de tables

Sans perte de généralité, on considère dans la suite uniquement l'évaluation de la fonction sin(x). Pour atteindre l'arrondi correct de cette fonction en un temps moyen raisonnable, on divise habituellement le schéma d'évaluation en deux phases. Une première phase rapide utilise des opérations d'une précision similaire à celle des opérandes pour garantir quelques bits supplémentaires par rapport à la précision visée. Si l'arrondi correct ne peut pas être déterminé après cette première phase, une phase lente basée sur une précision étendue est utilisée [START_REF] Daramy-Loirat | A library of correctly rounded elementary functions in double-precision[END_REF]. Chaque phase peut se décomposer en quatre étapes :

1. Une première réduction d'argument basée sur des identités mathématiques permet de réduire l'argument x en un argument

x * ∈ [0, π/4]. Précisons que x * est un nombre généralement irrationnel quand x est, lui, exact. La précision de x * conditionnant la précision du résultat final, il est alors indispensable de réaliser cette étape en utilisant un algorithme adapté [START_REF] Brisebarre | A New Range-Reduction Algorithm[END_REF][START_REF] Cody | Software manual for the elementary functions[END_REF][START_REF] Cody | Implementation and Testing of Function Software[END_REF][START_REF] Payne | Radian Reduction for Trigonometric Functions[END_REF]. L'argument réduit est alors éventuellement représenté par plusieurs nombres flottants dont la somme fournit une bonne approximation de x * . On trouve ici une première source d'erreur du processus d'évaluation.

2. L'intervalle de x * est encore trop grand pour approximer précisément f k (x * ) avec un polynôme de petit degré. Une deuxième réduction d'argument est alors effectuée, basée cette fois sur des données tabulées. Cette solution a été détaillée par Tang [START_REF] Tak | Table-Lookup Algorithms for Elementary Functions and Their Error Analysis[END_REF].

Elle consiste à diviser l'argument réduit x * en deux parties x * h et x * l , définies par :

x * h = x * × 2 p • 2 -p et x * l = x * -x * h . (1)
On se sert alors des p bits de poids fort de x h pour adresser une table de π/4 × 2 p entrées contenant des valeurs précalculées sin h ≈ sin(x * h ) et cos h ≈ cos(x * h ) et arrondies dans un format de destination, par exemple double-double. Le format doubledouble consiste à approximer un nombre z sur deux mots bi-nary64 qui ne se chevauchent pas. On utilise alors une des deux formules :

sin(x * ) = sin(x * h + x * l ) = sin h • cos(x * l ) + cos h • sin(x * l ) cos(x * ) = cos(x * h + x * l ) = cos h • cos(x * l ) -sin h • sin(x * l ). (2)
Cette étape introduit une deuxième source d'erreur, les approximations sin h et cos h .

3. Des approximations polynomiales de sin(x * l ) et cos(x * l ) sont alors évaluées. Ceci introduit une troisième source d'erreur.

4. Enfin, le résultat final est obtenu par une reconstruction qui utilise les termes obtenus dans les étapes 2 et 3 selon une des Équations [START_REF] Aad | Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC[END_REF]. Cette étape est encore une source d'erreur via les opérations en précision finie + et ×.

Il existe déjà des solutions satisfaisantes pour la première réduction d'argument, la génération et l'évaluation de polynômes d'approximation précis et efficaces, et la reconstruction [START_REF] Muller | Elementary Functions[END_REF].

xiv résumé substantiel en français La deuxième réduction d'argument a fait l'objet de propositions où l'objectif était d'équilibrer la première source d'erreur présente dans l'argument réduit et la deuxième source d'erreur présente dans les données tabulées. La méthode de Gal réduit les erreurs d'arrondis des valeurs tabulées, en autorisant un terme correctif corr dans les x * h [START_REF] Gal | Computing elementary functions: A new approach for achieving high accuracy and good performance[END_REF][START_REF] Gal | An accurate elementary mathematical library for the IEEE floating point standard[END_REF]. Pour chaque valeur de x * h , la variation corr est choisie de sorte que cos(x * h + corr) et sin(x * h + corr) soient très proches de nombres représentables en machine. Les données tabulées sont ainsi plus précises, ce qui simplifie les calculs intermédiaires. En revanche, il est nécessaire de stocker les termes correctifs. La recherche des valeurs corr a été améliorée par Stehlé et Zimmermann [START_REF] Stehlé | Gal's Accurate Tables Method Revisited[END_REF].

Dans le Chapitre 2, nous montrons comment aller plus loin et éliminer complètement la deuxième source d'erreur en tabulant des nombres entiers représentables sur un nombre flottant IEEE 754, qui donnent accès à des valeurs rationnelles proches de cos(x * h ) et sin(x * h ) ou de cosh(x * h ) et sinh(x * h ). Pour cela, nous nous appuyons sur les triplets pythagoriciens.

Tabulation de termes exacts

Notre objectif est de sélectionner des points pour lesquels il est possible de tabuler des valeurs exactes. Pour les fonctions trigonométriques et hyperboliques sinus et cosinus, nous utilisons les triplets pythagoriciens. La suite de la section présente le principe de notre approche, ainsi que les triplets pythagoriciens et la méthode proposée pour sélectionner les triplets pertinents pour la construction des tables.

Principe de la méthode

Notre approche consiste à construire des tables dont certaines valeurs sont stockées exactement. Plus particulièrement, rappelons qu'une table contient π/4 × 2 p entrées, et pour chaque entrée, les valeurs sin h et cos h doivent pouvoir être stockées sans erreur. Pour ce faire, nous imposons les conditions suivantes :

1. Les valeurs de sin h et cos h doivent être des nombres rationnels, c'est-à-dire :

sin h = s n /s d et cos h = c n /c d avec s n , s d , c n , c d ∈ N et s d , c d = 0.
2. Pour faciliter la reconstruction, les dénominateurs doivent être identiques : s d = c d . Cela permet de factoriser la division par les dénominateurs.

3.

Pour éviter l'évaluation de cette division durant la reconstruction, le dénominateur s d doit être égal à une même valeur k pour chaque entrée. Il est ainsi possible de l'intégrer dans les coefficients des polynômes. Cette valeur k peut alors être vue 0.2 tables exactes pour les fonctions élémentaires xv comme le plus petit commun multiple (PPCM) des dénominateurs de toutes les entrées.

4.

Pour réduire la taille des tables, s n et c n doivent être représentables en machine.

Après avoir trouvé des nombres qui vérifient ces quatre propriétés, la reconstruction devient :

sin(x * ) = s n • C(x * l -corr i ) + c n • S(x * l -corr i ),
avec :

• s n , c n et corr les valeurs tabulées définies par :

x * h = arcsin(s n /k) -corr i = arccos(c n /k) -corr i .

• C(x) et S(x) deux polynômes d'approximation définis pour x ∈

[-2 -p-1 + max i (corr i ), 2 -p-1 -min i (corr i )] par :

C(x) ≈ cos(x)/k et S(x) ≈ sin(x)/k.

Les valeurs sin h et cos h sont alors rationnelles. Mais seuls les deux numérateurs c n et s n seront stockés dans la table, les divisions étant précalculées dans les coefficients des polynômes.

Triplets pythagoriciens

Un triplet pythagoricien est un triplet d'entiers (a, b, c), où a, b et c vérifient :

a 2 + b 2 = c 2 ,
avec a, b, c ∈ N * . D'après le théorème de Pythagore, un triplet pythagoricien renvoie aux longueurs des côtés d'un triangle rectangle. Les valeurs de sin h et cos h peuvent alors être définies comme les quotients de ces longueurs. On en déduit qu'à tout triplet pythagoricien (a, b, c) peut être associé un angle θ ∈ [0, π/2] tel que sin(θ) = a/c et cos(θ) = b/c.

Un triplet pythagoricien pour lequel les fractions a/c et b/c sont irréductibles est appelé triplet pythagoricien primitif. Un triplet primitif et ses multiples renvoient à des triangles semblables, et représentent donc le même angle θ. Par conséquent, nous ne nous intéresserons ici qu'aux triplets pythagoriciens primitifs. Un exemple de triplet primitif est le triplet [START_REF] Abramowitz | Handbook of Mathematical Functions With Formulas[END_REF][START_REF] Anglin | Using Pythagorean Triangles to Approximate Angles[END_REF][START_REF] Tom | Introduction to analytic number theory[END_REF] qui renvoie à l'angle θ = arcsin(3/5) ≈ 0.6435rad, c'est-à-dire, environ 58°.

Construction et sélection de l'ensemble des triplets pythagoriciens

Il existe une infinité de triplets pythagoriciens primitifs, qui couvrent un large intervalle d'angles. Ceci est illustré sur la Figure 2.5 (p. [START_REF] Defour | A new scheme for table-based evaluation of functions[END_REF], qui montre tous les triplets pythagoriciens primitifs dont l'hypoténuse est strictement inférieure à 2 12 . L'ensemble des triplets primitifs possède une structure d'arbre ternaire [START_REF] Freek | Over Pythagorese en bijna-Pythagorese driehoeken en een generatieproces met behulp van unimodulaire matrices [On Pythagorean and quasi-Pythagorean triangles and a generation process with the help of unimodular matrices[END_REF][START_REF] Price | The Pythagorean Tree: A New Species[END_REF]. Ainsi, l'arbre de xvi résumé substantiel en français Barning-Hall [START_REF] Freek | Over Pythagorese en bijna-Pythagorese driehoeken en een generatieproces met behulp van unimodulaire matrices [On Pythagorean and quasi-Pythagorean triangles and a generation process with the help of unimodular matrices[END_REF] se construit à partir d'un triplet (a, b, c) dont on déduit trois noeuds fils en multipliant les matrices

  1 -2 2 2 -1 2 2 -2 3   ,   -1 2 2 2 -1 2 -2 2 3   et   1 2 2 2 1 2 2 2 3
  par ce triplet considéré comme un vecteur colonne. Tous les triplets primitifs peuvent donc être générés à partir du triplet [START_REF] Abramowitz | Handbook of Mathematical Functions With Formulas[END_REF][START_REF] Anglin | Using Pythagorean Triangles to Approximate Angles[END_REF][START_REF] Tom | Introduction to analytic number theory[END_REF], en prenant en compte la symétrie entre les triplets (a, b, c) et (b, a, c). Par exemple, après la première itération, on obtient les triplets [START_REF] Tom | Introduction to analytic number theory[END_REF][START_REF] Bernhart | Pythagoras' garden, revisited[END_REF][START_REF] Boehm | On de Casteljau's algorithm[END_REF], [START_REF] Boldo | A Simple Test Qualifying the Accuracy of Horner's Rule for Polynomials[END_REF][START_REF] Beadie | Refractive index measurements of poly(methyl methacrylate) (PMMA) from 0.4 -1.6 µm[END_REF][START_REF] Brisebarre | M, p, k)-Friendly Points: A Table-Based Method for Trigonometric Function Evaluation[END_REF], et [START_REF] Brunie | Modified Fused Multiply and Add for Exact Low Precision Product Accumulation[END_REF][START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF][START_REF] Cody | Implementation and Testing of Function Software[END_REF], et leur symétrique [START_REF] Bernhart | Pythagoras' garden, revisited[END_REF][START_REF] Tom | Introduction to analytic number theory[END_REF][START_REF] Boehm | On de Casteljau's algorithm[END_REF], [START_REF] Beadie | Refractive index measurements of poly(methyl methacrylate) (PMMA) from 0.4 -1.6 µm[END_REF][START_REF] Boldo | A Simple Test Qualifying the Accuracy of Horner's Rule for Polynomials[END_REF][START_REF] Brisebarre | M, p, k)-Friendly Points: A Table-Based Method for Trigonometric Function Evaluation[END_REF], et [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF][START_REF] Brunie | Modified Fused Multiply and Add for Exact Low Precision Product Accumulation[END_REF][START_REF] Cody | Implementation and Testing of Function Software[END_REF].

Parmi les triplets pythagoriciens primitifs, nous en sélectionnons un par entrée de la table, de telle sorte que chaque triplet (a, b, c) minimise la magnitude du terme corr défini par :

corr = θ -x * h , avec θ = arcsin(a/c) et x * h = i • 2 -p ,
où i représente l'indice de l'entrée dans la table (i 0). En effet, chaque écart x * h,k -θ k correspond au terme correctif corr k présent dans les évaluations polynomiales. En les réduisant, on réduit la taille de l'intervalle sur lequel on évalue les polynômes. Comme on l'a vu, au lieu de stocker a, b, c et corr, notre approche consiste à stocker A et B de la forme Nous avons implémenté cet algorithme en C++et nous avons inclus diverses optimisations pour accélérer les calculs et réduire la consommation mémoire. Par exemple, nous pouvons remarquer que le triplet dégénéré (0, 1, 1) permet de représenter l'angle θ = 0 • avec un terme correctif nul. Nous pouvons le choisir pour l'entrée d'indice 0 dans la table, ce qui nous permet de ne pas avoir à considérer cette entrée et ses triplets associés lors de la phase de génération.

Résultats numériques

La Table 2.1 (p.54) présente les résultats obtenus sur un serveur Linux CentOS (Intel Xeon CPU E5-2650 v2 @ 2.60 GHz, 125 Go de RAM). Ces résultats donnent les valeurs n et k min trouvées, les temps d'exécution, ainsi que le nombre de triplets et d'hypoténuses différents stockés. La table montre que l'on peut construire des tables indexées par p = 6 bits en 7 secondes. Pour p = 7 bits d'entrée, on passe alors à 31 secondes. Notons que dans cette dernière configuration, il est nécessaire de stocker en mémoire plus de 1.3 millions de triplets et 0.3 million d'hypoténuses. On observe que notre solution exhaustive n'est alors pas envisageable pour des valeurs de p 8.

Recherche heuristique

Afin de considérer, en un temps et une consommation mémoire raisonnables, des valeurs p plus grandes (p 8), nous avons cherché à caractériser les triplets et les PPCM déjà trouvés pour des p plus petits (p 7). La Table 2.5 (p.57) montre leur décomposition en facteurs premiers. En analysant cette table, On remarque que les PPCM obtenus sont toujours le produit de petits facteurs premiers, et plus particulièrement de petits nombres premiers pythagoriciens, c'est-à-dire, les nombres premiers de la forme 4n + 1 [START_REF] Fässler | Multiple Pythagorean number triples[END_REF]. Pour p 7, ces facteurs appartiennent à l'ensemble P des premiers pythagoriciens inférieurs ou égaux à 73 : 13, 17, 29, 37, 41, 53, 61, 73 . Nous avons alors développé une heuristique qui ne stocke que les triplets pythagoriciens primitifs dont l'hypoténuse est de la forme

P = 5,
c = i p r i i , avec p i ∈ P, et r i ∈ {0, 1} si p i = 5 r i ∈ N * sinon . (3)
Les résultats obtenus par cette heuristique sont consignés dans la Table 2.7 (p.58). La dernière colonne correspond au nombre d'hypoténuses comprises entre 2 n-1 et 2 n , et sélectionnées par l'heuristique. On observe que non seulement cette heuristique a une consommation mémoire nettement inférieure à notre algorithme exhaustif, avec optimisation, mais elle retrouve aussi les résultats obtenus par cet algorithme au moins jusqu'à p = 7, en des temps largement inférieurs. En effet, pour p = 7 par exemple, elle permet de ne stocker que 69 hypoténuses différentes et de générer la table en moins d'une seconde, contre plus de 0.3 million d'hypoténuses et 31 secondes pour notre algorithme exhaustif.

Le facteur limitant de cette heuristique algorithme n'est plus la recherche du PPCM, qui consiste à vérifier seulement quelques dizaines d'hypoténuses, mais la sélection des triplets. Nous avons donc développé une seconde heuristique utilisant une technique décrite par Fässler dans [START_REF] Fässler | Multiple Pythagorean number triples[END_REF]. Cette technique permet de construire des triplets pythagoriciens partageant une même hypoténuse. Les résultats de cette heuristique sont présentés dans la Table 2.9 (p.63). Ils montrent qu'on peut construire des tables dans la limite du stockage possible sur des nombres binary64 en un temps encore plus raisonnable.

Comparaisons avec les méthodes existantes

Nous avons présenté une réduction d'argument basée sur des points exacts ainsi qu'une méthode efficace pour générer ces points. Nous comparons cette solution aux solutions de Tang et de Gal. On considère un schéma d'évaluation de la fonction sinus en deux étapes qui cible l'arrondi correct en double précision. La phase rapide et la phase précise ont pour objectifs respectifs des précisions relatives d'au moins 2 -66 Concernant le nombre d'opérations flottantes, notre solution un gain jusqu'à 42% peut être noté. Enfin, on remarque que notre solution peut réduire le nombre d'accès mémoire. La phase rapide requiert 3 accès mémoire, tandis que l'approche de Tang en nécessite 4, ce qui représente une amélioration de 25%. Pour la phase précise, on passe de 6 accès pour Tang et 7 pour Gal à 5 accès pour notre méthode, c'est-à-dire, une amélioration jusqu'à 29%.

Conclusions et perspectives

Dans ce premier chapitre, nous présentons une nouvelle approche pour implanter la réduction d'argument pour l'évaluation de fonctions trigonométriques basée sur des valeurs tabulées. Notre méthode s'appuie sur les triplets pythagoriciens, ce qui permet de simplifier et d'accélérer l'évaluation de ces fonctions. Notre solution est la seule parvenant à éliminer l'erreur d'arrondi sur les approximations tabulées, pour la concentrer dans les approximations polynomiales. Nous réduisons ainsi jusqu'à 29% la taille des tables, le nombre d'accès mémoire, et le nombre d'opérations flottantes impliquées dans le processus de reconstruction.

Nous percevons deux axes de perspectives : Premièrement, il serait intéressant d'intégrer notre table de points exacts à une implémentation complète des fonctions sinus et cosinus pour observer son impact réel. Deuxièmement, une heuristique du même type que celle basée sur la technique de Fässler mais pour les fonctions hyperboliques serait souhaitable. Une piste serait d'utiliser les résultats décrits dans [START_REF] Rothbart | Pythagorean Triples: A New, Easy-to-Derive Formula With Some Geometric Applications[END_REF].

vectorisation de schémas polynomiaux

Les fonctions élémentaires sont souvent calculées à l'aide d'approximations polynomiales, dont l'efficacité dépend directement de celle du schéma d'évaluation sous-jacent. Le troisième chapitre de cette thèse montre que le schéma classiquement utilisé (Horner) est rarement le plus performant. En effet, d'autres schémas exploitent mieux les parallélismes des architectures modernes, en réduisant les dépendances de données. Ces résultats ont pour objectif d'être intégrés à un xx résumé substantiel en français générateur de code performant pour l'évaluation polynomiale dans le cadre de l'approximation de fonctions. Ils ont fait l'objet d'une publication à la conférence ComPAS 2016 [START_REF] De | Performances de schémas d'évaluation polynomiale sur architectures vectorielles[END_REF], tout en représentant une étape importante pour le Chapitre 4, qui a conduit à la soumission d'un article à la conférence ASAP 2018 [START_REF] De | Meta-implementation of vectorized logarithm function in binary floating-point arithmetic[END_REF] (en cours de relecture).

Contexte et Motivations

L'implémentation de fonctions élémentaires (sin, log, . . .) en logiciel passe souvent par une approximation polynomiale de degré raisonnable, sur un intervalle restreint [123, § 11.4]. L'évaluation d'un tel polynôme, c'est-à-dire l'ordre dans lequel les opérations sont effectuées, est fixé statiquement par un schéma d'évaluation dans le code, pour plusieurs architectures matérielles. Le schéma le plus utilisé est celui de Horner (illustré Figure 3.1 (p.75)), qui est essentiellement séquentiel. Sa latence est donc relativement élevée [START_REF] Munro | Optimal algorithms for parallel polynomial evaluation[END_REF]. Cependant, sélectionner manuellement un schéma plus rapide pour une architecture donnée peut s'avérer long et complexe.

Notre objectif est de générer automatiquement des schémas qui profitent au mieux des architectures modernes, notamment du parallélisme d'instructions (ILP) et du parallélisme de données. Ce chapitre présente une étude de la performance de différents schémas sur une architecture de type Intel Haswell 1 . Nous montrons deux résultats utiles : d'une part, des schémas exposant beaucoup d'ILP peuvent avoir à la fois une latence et un débit meilleurs que des schémas plus séquentiels, un résultat qui contredit les préconisations usuelles [123, § 11.5]. D'autre part, quand la production d'un résultat dénormalisé est très coûteuse, nous notons que l'intervalle d'évaluation prévu pour un polynôme donné est un facteur déterminant pour le choix d'un bon schéma.

Rappels sur l'évaluation polynomiale

Soit P un polynôme de degré n à coefficients flottants non nuls, comme définis dans la norme IEEE 754-2008 [START_REF][END_REF]. On veut évaluer P en un point x dans le même format. Dans cette thèse, on ne considère que les schémas d'évaluation sans adaptation de coefficients. (Plus de détails sur l'adaptation de coefficients peuvent être trouvés dans [START_REF] Eve | The evaluation of polynomials[END_REF][START_REF] Paterson | On the Number of Nonscalar Multiplications Necessary to Evaluate Polynomials[END_REF].)

Le schéma le plus courant est celui de Horner. Il permet d'évaluer P = n i=0 a i X i en n additions et n multiplications. Par exemple, pour n = 7, on aura :

P(x) = a 0 + x • (a 1 + x • (a 2 + x • (a 3 + x • (a 4 + x • (a 5 + x • (a 6 + x • a 7 )))))). ( 4
)
P(x) = (a 0 + x 2 • (a 2 + x 2 • (a 4 + x 2 • a 6 ))) + x • (a 1 + x 2 • (a 3 + x 2 • (a 5 + x 2 • a 7 ))). (5)
Ce schéma expose de l'ILP, comme illustré à la Figure

P(x) = ((a 0 + a 1 • x) + x 2 • (a 2 + a 3 • x)) + x 4 ((a 4 + a 5 • x) + x 2 (a 6 + a 7 • x)). (6)
Dans le cas d'un polynôme de degré n = 2 k -1, l'arbre d'évaluation d'Estrin est bien équilibré, mais il l'est un peu moins dans le cas général. Sa latence est en O(log n). Avec une architecture simplifiée infiniment parallèle ayant respectivement des latences de 3 et 5 cycles pour l'addition et la multiplication flottantes1 , pour n = 7 et en utilisant (4), ( 5) et [START_REF] Baker | Encyclopedia of Mathematics and its Applications[END_REF], on peut prévoir des latences de 56, 37 et 24 cycles, pour, respectivement, Horner, Horner-2 et Estrin. Avec un FMA en 5 cycles, on peut s'attendre, respectivement, à des latences de 35, 25 et 15 cycles. Il est donc important de choisir un bon schéma d'évaluation. Cependant, le nombre de schémas possibles au degré n croît rapidement avec n du fait de la combinatoire, ce qui rend impossible une recherche exhaustive dès que n 6 [144, ch.6]. D'autres schémas ainsi que la question de l'optimalité asymptotique sur parallélisme fini ou infini ont été étudiés dans [START_REF] Kiyoshi | On the Parallel Evaluation of Polynomials[END_REF][START_REF] Munro | Optimal algorithms for parallel polynomial evaluation[END_REF][START_REF] Muraoka | Parallelism Exposure and Exploitation in Programs[END_REF]. Dans le cadre de l'évaluation de fonctions élémentaires, nous avons comparé les performances de Horner à deux schémas classiques : Horner-2 et Estrin ; puis, à tous ceux de plus faible latence pour le degré 12 sur notre architecture simplifiée.

Comparaison des performances de schémas d'évaluation polynomiale

Dans cette section, nous analysons les résultats de mesures de performance obtenus pour les schémas de Horner, Horner-2 et Estrin, compilés avec GCC 5.2.0 (option -O2 activée), sur un processeur Intel Haswell bénéficiant des extensions AVX2 et FMA. Celles-ci opèrent sur 16 registres de 256 bits, en scalaire (sur un mot de poids faible) ou en vectoriel, et disposent d'un additionneur et de deux multi-xxii résumé substantiel en français plieurs/FMA [66, § 10.14]. GCC 5.2.0 est capable de générer des instructions vectorielles via une passe d'autovectorisation, que nous utilisons pour produire du code vectorisé.

Nous mesurons le débit inverse, c'est-à-dire le nombre moyen de cycles entre deux résultats. Ainsi, plus il est faible, plus le débit correspondant est grand. Pour un degré fixé, il dépend a priori du schéma et de l'architecture, que nous faisons varier en autorisant ou non l'utilisation du FMA. Sur ce type d'architecture, tant qu'aucun calcul ne dénormalise, nous observons que les schémas exposant le plus d'ILP ont un meilleur débit que les schémas séquentiels.

Afin que l'expérience soit représentative de l'évaluation de fonctions, nous faisons varier le degré du polynôme entre 3 et 32, et nous utilisons les coefficients de la série de Taylor de log(1 + x) en 0, calculés avec Sollya aucun des trois schémas ne souffre du nombre limité de registres, le séquenceur de GCC parvenant à conserver les résultats intermédiaires dans ces derniers. Cela est encourageant pour l'utilisation de schémas très parallèles dans le cadre de l'évaluation de fonctions.

En refaisant les mêmes expériences en double précision, on n'observe pas la perte de performance obtenue pour Estrin en simple précision. Cela conforte l'hypothèse que ceci est causé par un résultat intermédiaire dénormalisé. La diminution du débit avec le degré est alors quasi-linéaire, et le schéma d'Estrin a un débit toujours supérieur à Horner dès le degré 10 dans tous les cas. Nous concluons qu'une latence faible n'entraîne pas toujours une chute du débit, ce qui peut rendre des schémas comme Horner-2 ou Estrin plus intéressants si l'on veut privilégier la performance à la précision. Ce résultat est assez contre-intuitif car sur parallélisme fini, on pouvait s'attendre 1 Voir http://sollya.gforge.inria.fr/ et [START_REF] Sylvain Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF]. 0.3 vectorisation de schémas polynomiaux xxiii à une saturation des unités arithmétiques pour les schémas très parallèles.

Influence de l'intervalle d'évaluation

Dans la suite, on note RN(()x) l'arrondi dans une direction quelconque d'un réel x en simple précision. Les polynômes utilisés pour l'évaluation de fonctions élémentaires sont valides sur un petit intervalle. Les autres entrées sont soit traitées séparément, soit ramenées à l'intervalle d'évaluation. Cet intervalle est un voisinage du point d'approximation qui est souvent 0, c'est le cas pour ln(1 + x), exp(x), ou sin(x). Les entrées de petites magnitudes élevées à certaines puissances peuvent alors dénormaliser rapidement. Par exemple, pour

x ∈ [2 -12 ; 2 -5 ], x 16 dénormalise en simple précision dans environ (2 -126/16 -2 -12 )/2 -5 ≈ 13% des cas.

Les Figures 3.8 (p.81) et 3.9 (p.82) montrent les latences mesurées pour deux schémas d'élévation à la puissance 16 de nombres flottants simple précision dans le même intervalle. L'un correspond à une exponentiation naïve, l'autre à une exponentiation binaire. Leurs latences minimales sont de l'ordre de quelques cycles, mais quand

x devient suffisamment petit pour que x 16 dénormalise, on observe un palier pour les deux schémas. Les latences observées passent ainsi de quelques cycles à plus d'une centaine. D'après Agner Fog, cette pénalité provient de traitements micro-codés, pour toutes les opérations dont les opérandes sont normaux et dont le résultat est dénormalisé, ainsi que pour une multiplication dont un opérande est normal et l'autre dénormalisé [66, § 10.9]. Les paliers suivants observés pour la méthode naïve commencent quand les mônomes de degré inférieur dénormalisent. On remarque aussi que la version binaire ne retrouve une latence de quelques cycles que bien après que RN(()x 16 ) = 0, ce qui vient en partie contredire le manuel de Fog [66, § 10.9]. Cela suggère que le micro-code n'est exécuté que pour un intervalle précis de résultat. Par ailleurs, lorsque RN(()x 15 ) = 0, i.e. |x| 2 -10 , on peut voir un palier descendant pour la version naïve. Cela montre que la multiplication par zéro ne souffre pas de l'exécution de micro-code.

Ainsi, sur certaines architectures, le fait qu'un résultat intermédiaire dénormalise diminue grandement le débit des calculs, et notamment des évaluations polynomiales. Il est alors pertinent de s'assurer que, pour un schéma d'évaluation donné, l'intervalle d'entrée prévu n'engendre pas de dénormalisés. Ce résultat pourrait permettre d'améliorer la sélection automatique de schémas performants sur l'ensemble de son intervalle d'utilisation. Sur la Figure 3.10 (p.83), on remarque que tous les schémas de plus faible latence ont un meilleur débit que les schémas de Horner et Horner-2, et seuls une dizaine de schémas sont aussi bons qu'Estrin. Ce sont ceux qui impliquent le moins de multiplications. Assez logiquement, on observe aussi des paliers indiquant des débits plus faibles à mesure que le nombre de multiplications augmente. La Figure 3.11 (p.84) vient confirmer que le schéma d'Estrin peut être beaucoup plus performant que beaucoup de schémas de latence minimale, car il peut être en mesure de produire moins de résultats intermédiaires dénormalisés. En revanche, la Figure 3.12 (p.84) montre qu'Estrin est environ 12 fois moins bon que Horner et Horner-2 quand la dénormalisation intervient à partir de x 8 , alors qu'une cinquantaine de schémas générés restent plus performants que Horner et Horner-2, notamment parmi ceux qui impliquent peu de multiplications. Cela va dans le sens de nos conclusions précédentes, c'est-à-dire, qu'il faut s'assurer que le schéma que l'on choisit, pour une architecture qui lève des exceptions pour des résultats dénormalisés, ne produise pas lui-même de résultats dénormalisés sur son intervalle d'évaluation.

Une étude plus approfondie des meilleurs schémas générés par CGPE pourrait nous permettre de distinguer des motifs, réutilisables dans le cadre de l'approximation de fonctions par des polynômes univariés. Dans le cas contraire, une mesure automatisée des performances de ces schémas sur une architecture cible pourrait nous fournir des critères de choix. A posteriori, une borne sur la précision du résultat final pourrait être calculée afin de vérifier si un schéma sélectionné est conforme aux exigences demandées.

Conclusions

Dans le cadre de l'évaluation polynomiale pour l'approximation de fonctions, nous avons observé que certaines dénormalisations pouvaient fortement impacter les performances. C'est donc un facteur essentiel pour choisir un schéma performant sur tout son intervalle d'évaluation. Nous avons aussi vu que des schémas exposant beaucoup d'ILP comme le schéma d'Estrin ou ceux générés par CGPE ont souvent un meilleur débit que les schémas séquentiels comme Horner.

Ces résultats pourraient être la base de recherches plus approfondies. Il pourrait être intéressant d'intégrer des critères de choix dans CGPE pour éliminer les schémas qui, pour un intervalle d'évaluation et un ensemble de coefficients donnés, pourraient engendrer des 0.4 méta-implémentation de logarithmes xxv pertes de performances. Une piste serait d'utiliser Gappa à cette fin. Enfin, l'influence qu'a le choix d'un bon schéma sur les performances de l'évaluation de fonctions élémentaires mériterait d'être mesurée.

méta-implémentation de logarithmes

Afin de générer du code efficace, les environnements de développement devraient être en mesure d'exploiter toute source d'amélioration des performances. Par exemple, les extensions SIMD telles que SSE, AVX ou AVX-512 supportent des instructions vectorielles, permettant ainsi d'augmenter le débit d'applications exposant du parallélisme de données. En effet, ces instructions peuvent traiter actuellement jusqu'à 16 mots binary32 ou 8 mots binary64 avec un débit et une latence annoncés identiques aux instructions scalaires équivalentes. Le Chapitre 4 de cette thèse a pour objectif de montrer que la conception d'une implémentation de fonction élémentaire peut s'abstraire de plusieurs notions indépendantes des formats d'entrée/sortie ou de l'architecture visée en utilisant des générateurs de code tels que Metalibm-lugdunum. Dans le cadre de l'implémentation de fonctions élémentaires vectorielles, il n'est ainsi plus nécessaire d'écrire le code vectorisé à la main. En ajoutant du support pour la génération d'instructions vectorielles, c'est l'environnement Metalibm-lugdunum qui se charge entièrement de la vectorisation. Le développeur n'a alors plus qu'à concevoir des « méta-algorithmes » d'évaluation, vectorisables pour plus de performances.

Nous présentons dans ce chapitre une collaboration avec Nicolas Brunie1 , créateur et principal développeur de l'environnement Metalibm-lugdunum2 . Concrètement, nous présentons :

• une réduction d'argument classique mais réalisée sans branchements afin d'utiliser au mieux les ressources vectorielles de notre cible, AVX2 ;

• un « méta-algorithme » du logarithme naturel, avec arrondi fidèle, pour tout format de précision p 2 ;

• une méta-implémentation de cet algorithme dans l'environnement Metalibm-lugdunum, avec le support vectoriel nécessaire pour la cible AVX2 ;

• et des mesures de performances des implémentations générées pour les formats binary32 et binary64 sur une architecture AVX2.

Vectorisation de fonctions élémentaires

Depuis le milieu des années 1990, avec notamment le développement de l'imagerie en trois dimensions pour de nombreux domaines (médical, scientifique, architecture, divertissement. . . ), les jeux d'instruc-xxvi résumé substantiel en français tions SIMD se sont démocratisés [START_REF] Oberman | AMD 3DNow! technology: architecture and implementations[END_REF][START_REF] Peleg | MMX Technology Extension to the Intel Architecture[END_REF]. En effet, ces derniers proposent théoriquement de démultiplier le débit d'applications exposant un parallélisme de données massif en réalisant des instructions matérielles vectorielles, c'est-à-dire des instructions effectuant la même opération mais sur des données multiples.

Afin d'exploiter ce parallélisme pour l'évaluation des fonctions élémentaires, il convient d'utiliser des algorithmes qui se vectorisent efficacement. De tels algorithmes ne comportent habituellement pas de branchements ni de tables de correspondances qui dépendent des arguments. Dans ce chapitre, nous développons un algorithme vectorisable efficacement sur le jeu d'instruction SIMD AVX2, sans branchements mais avec une table de correspondances, pour l'évaluation du logarithme naturel en arithmétique flottante de précision arbitraire p 2. Cet algorithme est implémenté dans l'environnement Metalibm-lugdunum [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF][START_REF] Brunie | Code Generators for Mathematical Functions[END_REF] et permet de générer du code C scalaire ou vectoriel pour l'évaluation du logarithme naturel dans les formats binary32 et binary64, avec arrondi fidèle garanti. Cette précision est en effet suffisante pour certaines applications de physique théorique [START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF]. Il est possible, à faible coût, de générer des implémentations offrant un meilleur débit au prix d'une précision dégradée et non garantie, mesurable a posteriori, ce qui peut parfois être suffisant [START_REF] Hauth | Development and Evaluation of Vectorised and Multi-Core Event Reconstruction Algorithms within the CMS Software Framework[END_REF]. Comme notre algorithme ne dépend pas de la précision, il est envisageable de générer des implémentations pour d'autres formats comme binary16 ou binary128.

Conception d'un méta-logarithme vectorisable

On considère la fonction x → log(x) avec x un nombre flottant de précision p 2 tel que x > 0 et x ∈ {NaN, ±∞}. On pose

x = m • 2 e avec m ∈ [1, 2 -2 1-p ].
La réduction d'argument que nous utilisons est classique pour les implémentations scalaires. Soit r i une approximation de 1/m. Alors, en utilisant l'identité log(x) = e log(2) + log(m), le logarithme peut se calculer de la façon suivante :

log(x) = (e + τ) • log(2) -log(2 τ • r i ) + log(1 + u), avec u = r i • m -1, et τ = [m > √ 2]
. On introduit ici τ afin d'éviter une cancellation catastrophique qui peut arriver lorsque e = -1 et m ≈ 2, comme dans [START_REF] Revy | Automated Design of Floating-Point Logarithm Functions on Integer Processors[END_REF]. Cette réduction d'argument permet d'éviter des branchements conditionnels, ainsi que de profiter de l'approximation inverse rapide (rcp) souvent disponible sur les architectures vectorielles récentes. Quand x est un nombre dénormalisé, une technique couramment utilisée est une remise à l'échelle de x avant tout autre calcul [START_REF] Daramy-Loirat | A library of correctly rounded elementary functions in double-precision[END_REF]. On écrit donc x = x • 2 λ avec λ ∈ [0 . . p -1] le nombre de zéros en tête 0.4 méta-implémentation de logarithmes xxvii dans la représentation binaire de x. Alors on a x = m • 2 e 2 e min avec m ∈ [1, 2 -2 1-p ], et, finalement :

log(x) = (e + τ -λ) • log(2) -log(2 τ • r i ) + log(1 + u). (7)
La quantité r i devient donc une approximation de 1/m , tronquée sur i bits, et u = r i • m -1. La valeurlog(2 τ • r i ) est lue dans une table indexée par les i premiers bits fractionnaires de r i .

Nous montrons que u peut-être représenté exactement sous la forme d'un nombre flottant de précision p si le paramètre i est choisi de telle sorte que : Bien que les routines de la Libmvec se montrent plus performantes que nos implémentations, celles-ci ne sont précises qu'à 4 ulps quand les nôtres renvoient toujours un arrondi fidèle, à au plus 1 ulp du résultat exact. De plus, les routines de la Libmvec sont écrites une par une, directement en assembleur. L'avantage de notre méthode est de pouvoir générer en quelques secondes de nombreuses implémentations en C afin d'observer quels sont les meilleurs compromis, par exemple ici pour trouver quels paramètres donnent le débit le plus élevé.

i e max -2 and κ • (2 i + 2 -1 ) < 2 -1 , (8)

Conclusions

Dans ce chapitre, nous avons voulu montrer à travers l'exemple du logarithme les avantages offerts par la programmation générative. Nous avons effectivement montré comment générer de multiples implémentations vectorisées avec une garantie sur la précision du résultat final, et ce indépendamment des formats d'entrée/sortie. Des perspectives envisagées sont un support étendu du backend Intel vectoriel, notamment pour supporter les instructions AVX-512, le support de formats comme binary16 ou binary128, ou encore la métaimplémentation entièrement factorisée de logarithme en base β paramétrable.

conclusion générale

Cette thèse avait pour objectif le développement de techniques et d'algorithmes pour la génération de code performant appliquée à l'évaluation des fonctions élémentaires. Dans le Chapitre 2, nous avons présenté une méthode originale pour générer des tables de correspondances sans erreur utiles pour les fonctions trigonométriques et hyperboliques. Jusque là, aucune méthode, à notre connaissance, n'était parvenue à supprimer l'erreur sur les valeurs tabulées pour ces fonctions. conclusion générale xxix Dans le Chapitre 3, nous avons analysé l'impact sur le débit apporté par la vectorisation de schémas d'évaluation polynomiale. Premièrement, nous avons confirmé des résultats précédents montrant que certaines opérations en arithmétique flottante pouvaient sévèrement réduire le débit d'applications dès lors que des nombres dénormalisés sont impliqués. L'évaluation polynomiale fait partie de ces applications. Deuxièmement, nous avons observé que de nombreux schémas d'évaluation polynomiale exposant beaucoup de parallélisme d'instruction pouvaient avoir un plus grand débit que des schémas classiques (tels Horner) sur des architectures Intel Haswell.

Finalement, dans le Chapitre 4, nous avons défendu le concept de « méta-libms », c'est-à-dire des environnements de développement permettant de décrire une fonction mathématique (à divers niveaux d'abstraction) plutôt que de fournir des implémentations très spécifiques, et qui sont ensuite capable de emphgénérer du code source optimisé pour certains critères, comme les formats d'entréesortie, la précision, l'architecture cible, ou bien la vectorisation. Nous avons montré, à travers notre méta-implémentation du logarithme, nos contributions à l'un de ces environnments : Metalibm-lugdunum. Ces contributions incluaient notamment un plus grand support des générateurs de code pour les architectures Intel x86 SIMD, ce qui nous a permis de générer des implémentations du logarithme entièrement vectorisées en C portable ou pour les architectures AVX2. L'avantage de cette méthode est qu'elle permet de factoriser du code commun à de nombreuses implémentations, ce qui autorise des gains en temps de développement potentiellement très importants.

Perspectives

Dans le Chapitre 2, nous avons présenté une heuristique basée sur des algorithmes de Fässler, bien plus performante que notre méthode exhaustive, mais seulement pour les fonctions trigonométriques. Une heuristique similaire pour les fonctions hyperboliques serait intéressante et pourrait probablement se baser sur les travaux de Rothbart et Paulsell [START_REF] Rothbart | Pythagorean Triples: A New, Easy-to-Derive Formula With Some Geometric Applications[END_REF].

Par ailleurs, intégrer nos tables et mesurer leurs performances dans les implémentations en précision intermédiaires (4096 bits) de Johansson [START_REF] Johansson | Efficient Implementation of Elementary Functions in the Medium-Precision Range[END_REF] pourrait être intéressant pour estimer leur impact dans des programmes réels.

Dans le Chapitre 3, nos mesures de performances de schémas d'évaluation polynomiale nous font envisager plusieurs perspectives à court et à long terme. D'abord, il serait intéressant d'étudier l'impact sur les performances du choix d'un schéma pour l'évaluation d'une fonction élémentaire. Des tests préliminaires semblent montrer un faible impact avec des polynômes de petit degré. Cependant, au sein d'un environnement comme Metalibm-lutetia [START_REF] Brunie | Code Generators for Mathematical Functions[END_REF], se basant sur l'évaluation de nombreux polynômes, les conclusions pourraient être différentes.

xxx résumé substantiel en français Ensuite, il pourrait être intéressant de développer CGPE pour que celui-ci puisse éliminer des schémas si, sur un intervalle cible, ceux-ci peuvent produire des résultats dénormalisés. Gappa pourrait être utilisé pour implémenter une telle heuristique. CGPE pourrait alors garantir qu'un schéma donné ne subira pas de pertes de performances dues à des résultats intermédiaires dénormalisés, par exemple.

Une autre perspective serait de fusionner les générateurs de schémas d'évaluation polynomiale de CGPE dans les différents projets Metalibm. L'interface PythonCGPE présentée au Chapitre 4 permet déjà accéder à des générateurs de CGPE dans Metalibm-lugdunum, mais mériterait peut-être d'être développée.

Finalement, dans le Chapitre 1, section 1.3, mais aussi dans le Chapitre 4, nous avons présenté et défendu les concepts de « méta-libm » et plus généralement de « programmation génératrice » qui semblent prometteurs pour produire des bibliothèques mathématiques sur mesure pour les applications haute performance. De plus, cela pourrait réduire les temps de développement et de maintenance en favorisant la factorisation de codes numériques.

De potentielles perspectives pourraient donc être : la fusion de méta-implémentations similaires comme les logarithmes dans différentes bases, en suivant [START_REF] Revy | Automated Design of Floating-Point Logarithm Functions on Integer Processors[END_REF], dans Metalibm-lugdunum ; le support des formats binary16 et binary128 et des fonctions en précision mixte. Metalibm-lugdunum propose déjà des opérateurs FMA en précision mixte utiles pour, par exemple, l'apprentissage automatique [START_REF] Brunie | Modified Fused Multiply and Add for Exact Low Precision Product Accumulation[END_REF]. En ce sens, il pourrait être intéressant de dresser une liste des opérateurs et formats utiles mais manquants dans les différents projets Metalibm afin de les implémenter ; l'inventaire des méta-blocs qu'il serait intéressant de factoriser. Par exemple, dans Metalibm-lugdunum, il existe une méta-réduction d'argument de Payne et Hanek [START_REF] Payne | Radian Reduction for Trigonometric Functions[END_REF] qui pourrait probablement être abstraite un peu plus afin de pouvoir la réutiliser dans d'autres méta-fonctions ; et, finalement, un sujet de recherche à moyen terme pourrait être la détection automatique et en boîte noire de comportements asymptotiques dans Metalibm-lutetia. En effet, cet environnement fournit déjà un support pour détecter automatiquement et en boîte noire certaines propriétés de périodicité et de symétrie [START_REF] Brunie | Code Generators for Mathematical Functions[END_REF]. Cependant, en premier lieu, une approche en boîte blanche semble plus accessible, car pouvant raisonner sur des propriétés algébriques connues. Une telle analyse pourrait permettre d'améliorer l'algorithme de découpage en sous-domaines utilisé par Metalibm-lutetia. In 2012, a new elementary particle of mass 125 GeV was discovered independently by two CERN experiments: CMS (Compact Muon Solenoid), and ATLAS (A Toroidal LHC ApparatuS) [START_REF] Aad | Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC[END_REF][START_REF] Chatrchyan | Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC[END_REF]. 1 A bursting fever quickly spread in the general media that CERN was about to announce that this new particle could well be the famous Higgs boson. The Higgs boson is an elementary particle that was expected by the Standard Model of physics since the 1960's but that had never been experimentally observed with enough confidence. Up to now, in 2018, many experiments around the world have gone along with the hypothesis that indeed, this new particle is the one that physicists had been looking for for the past 50 years [START_REF] Aad | Measurements of the Higgs boson production and decay rates and coupling strengths using pp collision data at √ s = 7 and 8 TeV in the ATLAS experiment[END_REF][START_REF] Chatrchyan | Observation of a new boson with mass near 125 GeV in pp collisions at √ s = 7 and 8 TeV[END_REF][START_REF] Chatrchyan | Study of the Mass and Spin-Parity of the Higgs Boson Candidate via Its Decays to Z Boson Pairs[END_REF][START_REF]Evidence for the spin-0 nature of the Higgs boson using AT-LAS data[END_REF][START_REF] Khachatryan | Constraints on the spin-parity and anomalous HVV couplings of the Higgs boson in proton collisions at 7 and 8 TeV[END_REF][START_REF] Khachatryan | Precise determination of the mass of the Higgs boson and tests of compatibility of its couplings with the standard model predictions using proton collisions at 7 and 8 TeV[END_REF][START_REF] Sirunyan | Measurements of properties of the Higgs boson decaying into the four-lepton final state in pp collisions √ s = 13 TeV[END_REF]. This discovery was both comforting the Standard Model and disturbing our understanding of mass. It was such a breakthrough that the following year, Peter Higgs and François Englert, two of the six theoretical physicists who predicted the existence of the particle, received the 2013 Nobel Prize in Physics for the theoretical discovery they had made in 1964.

But what few media related was that the experimental discovery would probably have taken much more time without the computing power provided by CERN's Worldwide LHC Computing Grid. The CERN experiments produced 13 petabytes (13 × 10 15 B) of raw data in 2010 [START_REF] Brumfiel | High-energy physics: Down the petabyte highway[END_REF]. To give a perhaps more striking comparison, this represents a pile of Blu-Ray discs about twice as high as the Eiffel Tower. In addition, these data require complex processing to precisely reconstruct the high energy particle collisions, using state-of-the-art mathematical models. Such complex data processing is automated by deferring all the computations to specialized software on the Grid. However, as pointed out by CMS scientists, the time spent by the software evaluating elementary functions such as the exponential or the logarithm can represent a good share of the overall execution time [START_REF] Hauth | Development and Evaluation of Vectorised and Multi-Core Event Reconstruction Algorithms within the CMS Software Framework[END_REF][START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF].

In physics, another, perhaps simpler example where the sine function is used is the refraction of light at the interface of two media m 1 and m 2 in the Snell-Descartes law:

n 1 sin θ 1 = n 2 sin θ 2

INTRODUCTION

where n i and θ i are, respectively, the index of refraction and the angle of incidence of the medium m i . This well-known physical phenomenon is illustrated in Figure 1 at the interface between air and Poly(methyl methacrylate) (PMMA). Even with a rather high measuring error, we can see that θ air ≈ 60°and that θ PMMA ≈ 35°. Using the free-software "basic calculator" GNU bc, anyone 1 can devise an approximation of the refractive index of PMMA in the visible spectrum, assuming it does not depend on the wavelength and that the refractive index of air is 1:

echo 'pi=4 * a(1); n _ pmma = s(60 * pi/180)/s(35 * pi/180); \ n _ pmma' | bc --mathlib
This shell command line outputs a refractive index of 1.5099 . . . with 16 more digits. This value lays inside the bounds 1.485 and 1.510 found in the literature for the visible spectrum [START_REF] Beadie | Refractive index measurements of poly(methyl methacrylate) (PMMA) from 0.4 -1.6 µm[END_REF]. Simply said, the above command line asks bc to approximate π and output the ratio between sin(60°) and sin(35°). But under these apparently simple operations, bc had to approximate two sines (plus an arctangent to approximate π) and a few multiplications and divisions. As the result Figure 1 -Refraction of light at the interface between air and PMMA. Credit: Zátonyi Sándor, https://commons.wikimedia.org/ is displayed almost instantaneously, the naive user I was before the beginning of my Ph.D. would have thought that this computation was easy, and with only one final rounding. I had never been so wrong!

challenges of this thesis

The underlying challenges of this thesis are the design of efficient software implementations of elementary functions in floating-point arithmetic, with a particular focus on their optimization for current 1 Well, perhaps only GNU users, to be honest.

vector micro-architectures providing SIMD instructions. Indeed, after the basic arithmetic operations that are addition, subtraction, multiplication, division, and square root, elementary functions such as the sine and cosine functions, the exponential or the logarithm functions, are the second most basic resource in scientific computing. Their computation or approximation has a long history. It dates back to at least Babylonian mathematics, for which many tables of reciprocals and square roots have been found. The famous Plimpton 322 clay tablet (circa 1800 BC), which lists 15 Pythagorean triples, might even have served to obtain sine and cosine approximations, but this assertion is still debated [START_REF] Knuth | Ancient Babylonian Algorithms[END_REF]. From the end of XVIII th century to the beginning of XX th century, approximation theory was thoroughly developed, progressively leading to optimal polynomial and rational approximations to continuous functions [START_REF] Steffens | The history of approximation theory: from Euler to Bernstein[END_REF]. As polynomial and rational approximations are only made of elementary arithmetic operations, they represent a very efficient way to "compute" transcendental functions. However, even using these mathematical objects, efficiently and accurately approximating a function is already a pretty hard task to do, as we will see in Chapter 2 for trigonometric and hyperbolic functions, or in Chapter 4 for the logarithm function.

In the current context, most modern general-purpose processors are not designed to accurately approximate even basic elementary mathematical functions. On some general-purpose processors though, there exist hardware instructions (that is, code that the silicon can understand) that can approximate the sine or the cosine of a value. But we see at least three issues with these hardware instructions: 3. They are usually split into many micro-operations or µ-ops, i.e. basic hardware instructions such as memory loads or arithmetic instructions [START_REF] Fog | Instruction tables: Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD and VIA CPUs[END_REF][START_REF] Fog | The microarchitecture of Intel, AMD and VIA CPUs: An optimization guide for assembly programmers and compiler makers[END_REF]. This makes them very slow compared to more basic instructions.

Instead of using these problematic hardware instructions, an accurate computing software tool like bc or sollya1 will call software routines that may have more latency but may also be more accurate, if not correctly rounded. In Chapter 1, we will define such terms as correctly rounded, but for the time being, we only try to give the reader some intuitions about the evaluation of elementary functions. Although rather expensive, these software routines terminate in a very short time that is not perceptible to the human eye. We are talking in milliseconds here, but this is already much more than what a processor can do for basic tasks. Adding two numbers, for example, only takes a few nanoseconds, hence being about 1 000 000 faster.

Therefore, finding the right trade-off between performance and accuracy of elementary mathematical functions can be of paramount INTRODUCTION importance when millions of such computations are involved, e.g. in the CERN CMS experiment software. Hence, a current challenge is to develop techniques and tools to automate the process of writing mathematical libraries (called libm's) optimized for a given microarchitecture, instead of providing general purpose libraries. This way the user can generate his/her own libm for his/her specific problem, that is, for example, optimized for custom input ranges or hardware.

In this direction, a particular effort started in 2014 as the ANR1 MetaLibm project,2 whose goal is to develop such generation tools for mathematical functions and digital filters [START_REF] Brunie | Code Generators for Mathematical Functions[END_REF]. This project has involved PEQUAN (LIP6, Sorbonne Université / CNRS), Socrate (Citi, INSA Lyon / Inria), CERN, and DALI (LIRMM, Université de Perpignan Via Domitia / CNRS) teams. Members of MetaLibm tackled the above goal from different angles, described in more detail in Chapter 1. A first angle consists in open-ended code generators. Particularly, we consider a "function" (in every sense) as a black box to be implemented in floating-point arithmetic, and we produce source code valid for a specific input range [START_REF] Kupriianova | Towards a modern floating-point environment[END_REF][START_REF] Kupriianova | Metalibm: A Mathematical Functions Code Generator[END_REF]. A second angle addresses the design of function implementations that are (more or less) compliant with the standard libm, and optimized for different software and hardware [START_REF] De | Hardware Implementations of Fixed-Point Atan2[END_REF][START_REF] De | Meta-implementation of vectorized logarithm function in binary floating-point arithmetic[END_REF][START_REF] Quirin | A new open-source SIMD vector libm fully implemented with high-level scalar C[END_REF][START_REF] Le | Computing floating-point logarithms with fixed-point operations[END_REF][START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF]. A third angle is to focus on the particular case of digital filter implementations in fixed-point arithmetic [START_REF] Volkova | Reliable evaluation of the Worst-Case Peak Gain matrix in multiple precision[END_REF][START_REF] Volkova | Reliable verification of digital implemented filters against frequency specifications[END_REF][START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF]. This thesis takes place in the context of this project. Our work targets the generation of high performance codes under architectural constraints, for the implementation of mathematical functions. The contributions are detailed below.

contributions and outline of this thesis

This thesis is organized in six chapters, including this introduction, a state-of-the-art, and a conclusion. The end of this section is devoted to the brief presentation of each scientific chapter.

Chapter 1 -Computation of Elementary Functions

This chapter gives a brief overview of computer arithmetic. Particularly, it defines what we call floating-point arithmetic, and related objects and functions used in this thesis. Then, it presents state-ofthe-art designs for the evaluation of elementary functions, and code generation frameworks that can help automate performance tuning, with a particular focus on the techniques and tools used in this thesis.

Chapter 2 -Exact Lookup Tables for Elementary Functions

Various approaches have been developed to build trigonometric and hyperbolic functions. Among these, let us cite the table-lookup based approaches. These tables generally hold tabulated values of the functions to be implemented. However these values are approximations embracing rounding errors. Hence, their use makes it more complicated to guarantee a certain accuracy on the output.

This chapter presents a new approach to design "exact lookup tables", that is, where error has been removed, together with efficient heuristics to build them. These are based on mathematical objects called Pythagorean triples. First, some mathematical background on Pythagorean triples allows us to formally prove a lower bound on the main computed values. Second, we use relations existing between these functions and Pythagorean triples to give exact representations of tabulated values instead of intrinsically rounded values that have been used up to now. In comparison to Tang's and Gal's methods, two well-known ones, and when targeting correct rounding in double precision, we estimate theoretical gains of up to 29% and 42% in memory accesses and floating-point operations, respectively, during the reconstruction step.

This work has led to two publications in conferences, at Compas 2015 [START_REF] De | Réduction d'argument basée sur les triplets pythagoriciens pour l'évaluation de fonctions trigonométriques[END_REF] and ASAP 2015 [START_REF] De | Range reduction based on Pythagorean triples for trigonometric function evaluation[END_REF], and one publication in an international journal, IEEE Transactions on Computers [START_REF] De | Exact Lookup Tables for the Evaluation of Trigonometric and Hyperbolic Functions[END_REF].

Chapter 3 -Vectorizing Polynomial Schemes

Polynomial approximations take an important part in the evaluation of elementary functions [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. They can often be combined with tablelookup methods. The efficiency of their evaluation on a given microarchitecture can greatly impact the performance of the whole function implementation.

This chapter presents an analysis of the performances of classic polynomial evaluation schemes, namely Horner's and Estrin's, as well as of automatically-generated polynomial evaluation schemes, in the context of auto-vectorization. Polynomial evaluation schemes correspond to implementations of polynomial approximations. Autovectorization is the ability for a compiler to automatically transform scalar instructions into vectorized instructions that can treat multiple data simultaneously. We observe that on Intel Haswell architectures, the production of subnormal numbers as a result of arithmetic operations can increase the latency of such operations by a factor up to 50.

This work has led to one publication in a conference, at Compas 2016 [START_REF] De | Performances de schémas d'évaluation polynomiale sur architectures vectorielles[END_REF].

Chapter 4 -The Power of Meta-Implementations

Vectorization can increase the throughput of compute-intensive applications on supporting architectures [START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF]. The diversity of microarchitectures (Skylake, Haswell, Sandy Bridge, or Xeon Phi, . . . for x86 micro-architectures, for example) complicates the development

INTRODUCTION

of efficient implementations that exploit at best the specificities of the underlying hardware. To tackle this issue, we favor the use of code generation tools to automate the implementation process optimized for specific hardware.

This chapter concerns meta-implementations of elementary functions, vectorization, and code generation. We present the design and implementation process of a vectorizable natural logarithm within the Metalibm-lugdunum code generator framework. Particularly, it presents our contributions to the Metalibm-lugdunum framework, which consist in the extension of the Intel x86 backends, the introduction of reusable meta-blocks, the integration of PythonCGPE that enables the automatic generation of efficient polynomial schemes, and the vectorized meta-logarithm itself. Careful optimizations to avoid penalties due to subnormal numbers or to help the code generator produce vectorized code are presented for our meta-implementation. Comparisons with the free GNU Libmvec available in glibc 2.22 are reported and show that meta-implementations can be competitive while bringing decisive advantages such as easier code maintenance and reduced development times.

This work was done in collaboration with Kalray, and it has led to a submission to ASAP 2018 [START_REF] De | Meta-implementation of vectorized logarithm function in binary floating-point arithmetic[END_REF], which is currently under review. In this chapter, our goal is to give the reader the keys to understanding the rest of this thesis. Hence, let us recall a few basic bricks about how computers do what they are meant to do: computing. Since the late 1950's, transistors have formed the core components of electronic circuits. Transistors are active electronic components, which means they need an external source of energy to work. They can be used as electronic switches, in one of two states commonly called "on" and "off". These two available states make Boolean logic relatively easy to implement through logic gates. In turn, Boolean logic allows to build universal Turing machines, such as our modern computers [START_REF] Nisan | The elements of computing systems: building a modern computer from first principles[END_REF]. A universal Turing machine is a special kind of theoretical computer that can compute anything that is computable, provided an infinite amount of memory. Thus, if we forget the "infinite memory" requirement of a truly universal Turing machine, modern computers, because they implement Boolean logic, are fundamentally able to compute anything that is computable. Computability theory is a fascinating but complex field of research that is not the subject of this thesis. Hence, we will simply admit that the elementary functions we deal with in this thesis are all computable for all their (computable) inputs. The interested reader can find more details in [START_REF] Ker | Complexity Theory of Real Functions[END_REF].

Back to our two-state transistors, Boolean logic makes binary arithmetic easy to implement, which is why nowadays, the vast majority of computers use binary digits (bits, whose value can be either 0 or 1) to convey information and perform basic computations such as addition or multiplication [START_REF] Cowlishaw | Decimal Floating-Point: Algorism for Computers[END_REF]. However, this has not always been the case. Indeed, many civilizations have embraced the radix-10 number system probably because they have been used to counting on their two pentadactyl ("five-finger") hands. 1 Therefore, most calculating tools and mechanical computers, such as abaci or the Pascaline, used radix 10 or radix 5. In the early development of two-state electronic computation of elementary functions calculators, many early computers would still internally use decimal digits coded on four bits (a representation format generally called Binary Coded Decimal or BCD), mostly for commercial applications [75, p. 20, 31]. Most handheld calculators, e.g. the Texas Instruments TI-89 and TI-92 calculators, Casio calculators, and Hewlett-Packard calculators, are still using decimal arithmetic internally [START_REF] Cowlishaw | Decimal Floating-Point: Algorism for Computers[END_REF][START_REF] Motorola M | FAMILY Programmer's Reference Manual[END_REF]. Financial software often do the same, usually through programming language facilities. Such an arithmetic is still useful when exact representation of radix 10 decimal numbers is required, because binary arithmetic generally cannot represent these numbers exactly, as we will see in Section 1.1.

The memory of a computer is generally addressed by words made of a finite number of bits. Thus, the set of representable real numbers on one word is essentially finite. Therefore, only finite subsets of real numbers can fit in a machine word, and allow approximations of real-number arithmetic. Among the many representation systems used for these approximations on modern computers, the one required by the IEEE 754 standard for floating-point arithmetic 1 has been the most pervasive since its adoption in 1985 [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. Hence we begin this chapter with a short introduction to floating-point arithmetic, and especially its IEEE 754-2008 standardization, because it is the arithmetic we use all along this thesis. Then we present general ideas for the design and implementations of efficient and accurate elementary functions, as well as pitfalls of approximating transcendental functions. Finally, we focus on state-of-the-art efforts that deal with assisted or automatic generation of high-performance code for the evaluation of elementary functions.

floating-point arithmetic

We are accustomed to using fixed-point arithmetic on a daily basis: when calculating a bill, or measuring volumes or weights, the radix point usually keeps a fixed position. Consider the following fixedpoint computation in radix 10: 42 + 0.17 = 42.17.

Here, the decimal point is fixed before and after the addition: it always indicates where the integer part of a number ends and where its fractional part begins. However, when dealing with very large or very small numbers, it is much more common to use the scientific notation, where the radix point can float. Consider Coulomb's constant. It is the proportionality factor in Coulomb's law, which describes the force of interaction between two electrically charged particles. The constant is noted k e and it can be made exact [START_REF] Dugdale | Essentials of Electromagnetism. Macmillan physical science[END_REF]:

k e = 8 987 551 787.368 176 4 N m -2 C -2 .
1 IEEE stands for the Institute of Electrical and Electronics Engineers. More info can be found at https://ieee.org/about/ and in [START_REF][END_REF]77] for the 754 standards.

Notice how the scientific notation 8.987 551 787 368 176 4 × 10 9 N m -2 C -2 presents the same information, but is easier to read and understand. Indeed, since the position of the decimal point is embedded in the exponent 9, we immediately know the order of magnitude of the constant compared to the fixed point notation, which requires counting the number of decimal digits. Moreover, in many fields, it is often more useful to know the order of magnitude of a numerical result than all its digits. Finally, multiplication using the scientific notation is easier: only a multiplication with small significands is needed, while the exponents are simply added together. For instance, two elementary charges1 that are one meter apart are subject to a force of interaction, the magnitude F of which is given by Coulomb's law:

F = k e × e 2 ≈ (8.99 × 10 9 ) × (2.57 × 10 -38 ) = (8.99 × 2.57) × 10 9-38 = 23.1043 × 10 -29 N.
Next, a renormalization step and a rounding step can be performed so that F be rewritten F ≈ 2.31 × 10 -28 . In this example, the renormalization step shifted the decimal point and adjusted the exponent value (-29 was incremented to -28) so that the significand 2.31043 stay between 1 and 10. Then the rounding step kept only three significant figures, rounding the significand to the nearest 3-digit number 2.31.

Floating-point arithmetic is in many ways very similar to scientific notation. It consists in representing a real number x as the product of a nonnegative significand (or mantissa) m 0 and an integral power (the exponent) e of a radix β. To be able to represent negative numbers, a sign bit s ∈ {0, 1} is added to the representation so that finally:

x = (-1) s • m • β e . ( 1.1)
In this thesis we mainly consider β = 2, also known as binary floating-point arithmetic, in floating-point environments conforming to the IEEE 754-2008 standard [START_REF][END_REF]. However, whenever we believe it can be useful, we try to give intuitive examples in decimal arithmetic (β = 10). In the following, we present only the IEEE 754-2008 features that are useful in this thesis:

• floating-point formats and representation,

• required operations and rounding modes,

• and recommended correctly rounded functions. • its precision p 2, which is the maximum number of significant digits of m;

• the minimum exponent e min ;

• and the maximum exponent e max .

Since the 2008 revision requires that e min = 1 -e max for all formats, the above list can be reduced to p and one of the exponents extrema. Furthermore, floating-point numbers that conform to such a format are required to be normalized so that their representation be unique. In this case, for a precision-p floating-point number x = (-1)

s • m • 2 e : • either m = 1.m 1 • • • m p-1
and e min e e max , and x is called a normal number;

• or m = 0.m 1 • • • m p-1 and e = e min , and x is called a subnormal number.

For instance, in a radix-10 floating-point format with precision 3 and e min = -2, the number -123 = -1.23 • 10 2 is normal, while the number 0.0042 = 4.2 • 10 -3 = 0.42 • 10 -2 is subnormal. Note that subnormal numbers have at most p -1 significant digits. Subnormal numbers allow for a gradual flush towards zero, preserving some properties but sometimes at much more expensive computation costs as we will see in Chapter 3.

The two formats that we use in this thesis are the binary32 (sometimes called single precision) and the binary64 (sometimes called double precision) floating-point formats. The parameters specifying these formats as well as the binary128 format are summarized in Table 1.1.

As we write this thesis, the latest version of the C programming language (ISO C11) provides only binary32 and binary64 formats as

S E T (sign) (biased exponent) (trailing significand field)
1 bit w bits p -1 bits Most significant bits Least significant bits basic types with the float and double keywords, respectively [START_REF]International Organization for Standardization[END_REF].

Although a double-extended format can be available on some systems and compilers via the long double type specifier, ISO C11 only requires it to be at least as precise as the binary64 format. However, some compilers such as GCC (at least since version 4.8) provide the binary128 format on some target architectures as a language extension, via the __ float128 or _ Float128 types. On x86_64 architectures, which have no hardware support for binary128, arithmetic operations on binary128 words are usually carried out in software. 1Finally, to fully characterize a floating-point format, the IEEE 754-2008 standard also specifies binary interchange formats, i.e. the representation of binary floating-point numbers in machine words:

• a sign bit S ∈ {0, 1};

• a w-bit biased exponent E ∈ [0 . . 2 w -1]
, where [a . . b] denotes the set of integers no less than a and no greater than b;

• and a (p -1)-bit trailing significand field

T = m 1 m 2 • • • m p-1 .
This is summarized in Figure 1

• the normal number (-1) S • (1 + T • 2 -p ) • 2 E-e max if E 1 and T = 0;
• or the subnormal number (-1) S • T • 2 e min +1-p if E = 0. Thus the signed zero ±0 is represented when T = 0.

In the case E = 2 w -1, the interchange format represent special values:

• the signed infinity ±∞ if T = 0;

• or "not a number" (NaN) if T = 0. NaNs2 are used to signal invalid operations without interrupting the program flow, and therefore propagate through arithmetic computations 3 .

In this thesis, we indifferently denote by F the different sets of values that can be represented by binary interchange formats.

Before the wide adoption of the IEEE 754 standard, there were various floating-point arithmetics that were often incompatible with one another [123, § 1.1]. For instance, the UNIVAC 1107 used radix-2 with a precision p = 27 and an exponent of width w = 8 biased by 128 [168], whereas the IBM System/360 provided an environment with β = 16, p = 24 or 56 and w = 7 with a bias of 64 for its single and double precision formats, respectively [START_REF] Ibm | A Programmer's Introduction to IBM System/360 Assembler Language[END_REF]. Furthermore, these architectural differences were the cause of major portability and results reproducibility issues, as well as difficulties in writing numerical algorithms [START_REF] Kahan | Why do we need a floating-point arithmetic standard?[END_REF]. The need for standardization was met with an involvement of the main manufacturers and academics in the redaction of the first IEEE 754 standard for binary floating-point arithmetic proposal, which was accepted in 1985 [77]. This version of the standard was quickly adopted by hardware manufacturers so that, in 2018, most if not all floating-point units are IEEE-754 compliant. The 2008 version presented in this section is a major revision of the 1985 one. It standardized the Fused Multiply and Add (FMA) operation, as well as quadruple precision (binary128), and included decimal floatingpoint arithmetic, among other changes. A minor revision is being prepared for release in 2018. 1

IEEE 754-2008 required operations and rounding attributes

Basic arithmetic operations such as addition, subtraction, multiplication, division, square root, and FMA are required to be correctly rounded by the IEEE 754-2008 standard. An operation is said to be correctly rounded if it always return the rounding of its infinitely precise result according to a rounding-direction attribute. There are four required rounding-direction attributes for the correct rounding property in IEEE 754-2008 binary floating-point arithmetic, illustrated in Figure 1.2: roundtiestoeven (rn) Return the floating-point number closest to the infinitely precise result. If the mathematical result is exactly halfway between two floating-point numbers, return the floating-point number that has an even least significant digit.

roundtowardpositive (ru) Return the floating-point number closest to and no less than the infinitely precise result.

roundtowardnegative (rd) Return the floating-point number closest to and no greater than the infinitely precise result.

roundtowardzero (rz) Return the floating-point number closest to and no greater in magnitude than the infinitely precise result.

Hence, whenever an atomic operation such as a floating-point addition is performed, a rounding occurs. For a computed result ŷ corre- sponding to the mathematical, infinitely precise result y, the rounding error can be quantified by the absolute error

1 See http://754r.ucbtest.org/. 0 ∞ Floating-point numbers (a) RoundTiesToEven 0 ∞ (b) RoundTowardPositive 0 ∞ (c) RoundTowardNegative 0 ∞ (d) RoundTowardZero
δ = |y -ŷ| .
Often, we may use the ulp function to quantify the absolute error of the operation. ulp stands for "unit in the last place". Several definitions of the ulp function coexist, but we choose Definition 1, from [123, p. 33]:

Definition 1 (Goldberg's definition extended to real numbers). If |x| ∈ [2 e , 2 e+1
), then ulp(x) = 2 max(e,e min )-p+1 .

For example, in binary32,

42 ∈ [2 5 , 2 6 ) =⇒ ulp(42) = 2 -18 .
This means that the smallest binary32 number greater than 42 is 42 + ulp(42) = 42 + 2 -18 , and that their midpoint (which is not a binary32 number) is 42 + 1 2 ulp [START_REF] Defour | Cache-Optimised Methods for the Evaluation of Elementary Functions[END_REF]. Note that ulp(0) = 2 e min -p+1 = 2 -149 for the binary32 format. Using Figure 1.2, one can get the intuition that the error δ of a correctly-rounded arithmetic operation returning x is such that • |δ| 1 2 ulp(x) for the RN attribute; • or |δ| ulp(x) for the RU, RD or RZ attributes.

Sometimes, e.g. for a sequence of operations, we may measure the absolute error, or the relative error of the computed result [START_REF] Goldberg | What Every Computer Scientist Should Know About Floating-point Arithmetic[END_REF]. For the latter, we use Definition 2: Definition 2 (Relative error). If x ∈ R * is approximated by a finiteprecision floating-point number x, then the relative error of x with respect to x is

x -x

x .

This metric informs us about the approximate number of significant digits in the computed result of an operation: if the relative error is close to β -n , then there are roughly n radix-β significant digits in the result. For example, x = 0.1 is not representable in binary32. Indeed, RN(0.1) = 0.100 000 001 490 116 119 384 765 625.

Hence, the relative error due to rounding 0.1 in binary32 is

0.1 -RN(0.1) 0.1 ≈ 1.4999 • • • × 10 -8 .
This means that RN(0.1) has roughly eight (slightly less, actually) significant decimal digits. This metric is often used in the context of elementary functions, in conjunction with the absolute error. In the 2008 standard, elementary functions are part of recommended operations, which are the topic of the next subsection. The major issue arising with transcendental functions is their correct rounding, as it is hard to tell in advance how much accuracy is required to ensure this property [START_REF] Vincent Lefèvre | Toward Correctly Rounded Transcendentals[END_REF]. Hence, as we will see in Section 1.2, many libraries have been trading accuracy for performance by providing faithfully-rounded elementary functions. Faithful rounding is the property for an arithmetic operation whose infinitely precise result is x to return one of RU(x) or RD(x). Note that it is not a rounding attribute. Therefore the absolute error δ of a faithfully rounded result x is such that |δ| ulp(x).

The difficulty to correctly round elementary functions is related to the Table Maker's Dilemma, which we briefly explain in the next section. Then we present classic designs of elementary functions evaluation, as well as the recent efforts to design efficient vectorized functions by exploiting the available low-level parallelism of modern processors architectures.

evaluation of elementary functions

As we have seen in the Introduction, elementary functions are pervasive in many scientific or financial fields, and their evaluation of-

logp1 log e (1 + x) [-1, +∞] pow(x, y) x y [-∞, +∞] × [-∞, +∞] sin sin(x) (-∞, +∞) cos cos(x) tan tan(x) sinPi sin(π • x) cosPi cos(π • x) asin arcsin(x) [-1, +1] acos arccos(x) atan arctan(x) [-∞, +∞] atanPi arctan(x)/π atan2(y, x) arctan(y/x) [-∞, +∞] × [-∞, +∞] atan2Pi(y, x) arctan(y/x)/π sinh sinh(x) [-∞, +∞] cosh cosh(x) tanh tanh(x) asinh arsinh(x) [-∞, +∞] acosh arcosh(x) [+1, +∞] atanh artanh(x) [-1, +1]
ten require tradeoffs between accuracy and performance. However, since most are transcendental functions, they are subject to the Therefore, the accuracy required to correctly round an approximation of a transcendental number can be much higher than the target precision.

For instance, let us take the cosine function cos(x), with, say, the RU rounding direction attribute. This example was given in [START_REF] Vincent Lefèvre | Toward Correctly Rounded Transcendentals[END_REF]. The first decimal digits of cos(1.149) are 0.409 400 000 428 5

• • • . (1.2)
Therefore, in precision p = 4 decimal digits, the correct rounding towards positive (RU) of cos(1.149) is 0.4095. But because the exact value is a transcendental number, it cannot be computed in a finite number of steps. Hence, let us assume that an intermediate result ŷ was first computed on twice as many digits as p with an absolute error δ such that |δ| 10 -6 . In this thesis, we abstract absolute error notation by writing δ i to denote an absolute error less than or equal to β i in magnitude, e.g. |δ -6 | 10 -6 in this example. Therefore let δ = δ -6 . Then we have

cos(1.149) = ŷ + δ -6 .
If, e.g. ŷ = 0.409 399 91, then δ -6 = -8.042 858 708 1 • • • × 10 -8 . We see that RU( ŷ) = 0.4094, which is not the correctly rounded result. Moreover, without knowing the exact value of δ -6 , we know that cos(1.149) ∈ ŷ -10 -6 , ŷ + 10 -6 = [0.409 398 91, 0.409 400 91] .

From this equation, we deduce that the correct rounding (RU) of cos(1.149) is either 0.4094 or 0.4095, but as is, we cannot decide. The absolute error of the intermediate result ŷ is too high. From the digits in (1.2), we see that computing an approximation within an error bound δ -10 such that |δ -10 | 10 -10 would be enough to correctly round cos (1.149). But for other inputs, the maximum error bound may be less than δ -10 ! Hence a natural question comes to mind: "What is the maximum absolute error bound δ -m on ŷ such that for all possible x, RU( ŷ) = RU(cos(x))?" This is what the TMD is about: solving the TMD roughly means finding the minimal integer m answering such a question for a given format and rounding attribute, for all possible inputs of a function.

For the binary32 format, exhaustive search can solve the TMD in reasonable time. But as the TMD is exponentially hard to solve with the size of the inputs, exhaustive search is much more compute intensive for the binary64 format. Fortunately, clever algorithms have been designed to solve the TMD for many elementary functions in binary64 [START_REF] Gouicem | Conception et implantation d'algorithmes efficaces pour la résolution du dilemme du fabricant de tables sur architectures parallèles[END_REF][START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF][START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF][START_REF] Vincent Lefèvre | Toward Correctly Rounded Transcendentals[END_REF]. It is now known that approximating most functions on roughly three times the format precision will allow for correct rounding. Yet, this extra-precision is seldom needed and has a substantial overhead. That is why correct rounding is often achieved using a strategy in several steps known as Ziv's onion-peeling, which can be explained by the stochastic or random model.

In the random model, it is assumed that an elementary function behaves like a uniform pseudo-random bit generator, i.e. that every bit of a result has a probability 1/2 to be either a zero or a one. This is not true, but experiments show that this simple model is in fact quite accurate [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. A direct consequence of this model is that any bit string of length n has a probability 2 -n to appear at any position in the result. In particular, if {d n } represents a bit string made of n bits set to d, then the bit strings {0} n , {1} n , 1 {0} n-1 and 0 {1} n-1 have probability 2 -n to happen after the bits that are to be stored. According to the random model, this means that correct rounding may be achieved for any rounding direction with probability 1 -2 -n , provided a temporary result is computed with n digits more than the targeted precision. For example, in binary64, if an intermediate result is computed on 66 bits, which means 13 additional bits compared to the target precision p = 53, the random model predicts that about 1 result in 2 13 = 8192 will be rounded with no guarantee. This corresponds to a guaranteed rounding rate near 99.999%.

Ziv followed this model to implement his correctly rounded mathematical library [START_REF] Ziv | Fast Evaluation of Elementary Mathematical Functions with Correctly Rounded Last Bit[END_REF]. The principle is intuitive. First, compute an approximation with n extra bits of precision, which allows correct rounding most of the time. If correct rounding is not possible, compute a new approximation on more bits. As the algorithm was not proven to always terminate in a reasonable amount of operations, Ziv limited the extra precision of his last step, so that the probability of rounding incorrectly was astronomically small [START_REF] Florent De Dinechin | Towards the Post-Ultimate Libm[END_REF].

The CR-Libm effort combined Ziv's strategy with the solutions to the TMD to provide a binary64 correctly-rounded libm [START_REF] Daramy-Loirat | A library of correctly rounded elementary functions in double-precision[END_REF]. They pushed for the concept of only two steps in Ziv's algorithm: a first "quick" step that allows correct rounding most of the time, and, following the TMD solutions, an "accurate" step that guarantees correct rounding for hard-to-round cases. Such design is now very common for correctly-rounded transcendental functions. The decision of running the accurate phase relies on a simple test involving the compu-tation of a tight error bound [START_REF] De | On Ziv's rounding test[END_REF][START_REF] Torres | Tools for the Design of Reliable and Efficient Functions Evaluation Libraries[END_REF]. In the next section, we will see different methods of approximations that allow such "quick" and "accurate" steps in elementary function evaluations.

Algorithms for the evaluation of elementary functions

As mentioned earlier, elementary functions recommended by the IEEE 754 standard are mostly transcendental. This entails that they generally need an unbounded number of arithmetic operations to be evaluated. Therefore, different techniques have been developed to approximate such functions with a finite number of arithmetic operations. In this section, we present different algorithms that serve the purpose of approximating elementary functions.

Approximation methods

One of the simplest and earliest methods to approximate a function is to store precomputed values in mathematical tables, and to look up the value corresponding to an input when needed. This technique was used centuries ago by the Babylonians, in Greece, and in Ancient India, to get fast approximations of the basic trigonometric functions. The precomputation step was obviously the most difficult one. It may rely on mathematical identities, e.g. cos(2θ) = 1 -2 sin(θ) 2 , applied to already known values such as sin(π/4) = √ 2/2 [START_REF] Abramowitz | Handbook of Mathematical Functions With Formulas[END_REF], or on rational approximations, e.g. using Pythagorean triples, as we will see in Chapter 2. During the Renaissance, many handbooks of mathematical functions provided mathematical tables that were used to approximate a function by hand, possibly using interpolations. They generally provided values for equally spaced inputs, as a linear approximation would provide approximately the same accuracy whichever the input. For instance, Bernegger and Meyer's mathematical tables provided sine, tangent, and secant values for all arc minutes in

[0°, 90°] [11].
For digital computers, table-lookup algorithms have been used in many software and hardware implementations [START_REF] Defour | A new scheme for table-based evaluation of functions[END_REF][START_REF] De | Multipartite table methods[END_REF][START_REF] Michael | Approximating elementary functions with symmetric bipartite tables[END_REF][START_REF] Tak | Table-Driven Implementation of the Logarithm Function in IEEE Floating-Point Arithmetic[END_REF][START_REF] Tak | Table-Lookup Algorithms for Elementary Functions and Their Error Analysis[END_REF][START_REF] Wang | k)-Friendly Points: A Table-Based Method to Evaluate Trigonometric Functions[END_REF]. Used in conjunction with mathematical range reductions (see Section 1.2.2), they make a trade-off between computation costs and memory. The larger the table, the smaller the approximation error. But for performance reasons and because L1 data (L1d) caches are small, tables typically do not hold more than a few thousand rows [START_REF] Defour | Cache-Optimised Methods for the Evaluation of Elementary Functions[END_REF]. Although the mathematical range reduction can narrow the approximation interval so that a medium degree polynomial may be used, it is sometimes desirable to further reduce the interval of approximation so as to have very few arithmetic operations to perform. Table-based range reduction schemes allow for such additional reduction. They simply rely on a time/memory trade-off by storing precomputed values of the function(s) needed to reconstruct the final result.

In 1991, Tang suggested to use regularly spaced tables in hardware implementations of elementary functions that are already reduced to medium size intervals [START_REF] Tak | Table-Lookup Algorithms for Elementary Functions and Their Error Analysis[END_REF]. In the three examples given in his article, he uses equally spaced values for 2 x on [-1, 1] and log(x) on [START_REF] Aad | Measurements of the Higgs boson production and decay rates and coupling strengths using pp collision data at √ s = 7 and 8 TeV in the ATLAS experiment[END_REF][START_REF] Aad | Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC[END_REF], while a different yet regular tabulation is used for sin(x) on [0, π/4]. Indeed, for the trigonometric function, Tang's tabulated values are indexed by two variables j and k such that the input "breakpoints" are of the form 2 -j (1 + k/8), for j ∈ [1 . . 4] and k ∈ [0 . . 7]. By choosing the breakpoint that is closest to the input, the interval of approximation is reduced to [-1/32, 1/32]. This small interval allows for polynomial approximations of small degree given the targeted precision. At the same time, the unequal breakpoint spacing may account for finer-grained values when close to zero, where the sine reaches its maximum rate of change.

On one hand, such a tabulation allows for optimizations in the computation of the reduced arguments and in addressing the table, since the breakpoints are easily computable. On the other hand, it is not really able to optimize the accuracy of all tabulated values. Indeed, since most elementary functions suffer from the TMD, images of rational values may fall very close to a midpoint of their storage format. Therefore the error on such tabulated values, provided they are correctly rounded to the nearest in their storage format, may be as much as 0.5 ulp [START_REF] Goldberg | What Every Computer Scientist Should Know About Floating-point Arithmetic[END_REF]. Other table-based range reductions have been suggested in order to increase the accuracy of tabulated values, such as Gal's accurate tables.

In 1985, Gal suggested a method of tabulation for elementary functions that allows a virtually increased accuracy for tabulated values [START_REF] Gal | Computing elementary functions: A new approach for achieving high accuracy and good performance[END_REF]. By using the random model described in Section 1.2.1, Gal shows that for any function f conforming to this model and for any p > 0, there is approximately one x in any set of 2 n values such that f(x) has a string of n identical bits after its first p-th bits. Gal's probabilistic hypothesis was comforted by picking pseudorandom samples around regularly-spaced values so that a table of almost regularly-spaced values could be built with virtually extended accuracy. Although the precomputation cost is exponential in the number of virtual bits gained for each tabulated value, Gal and Bachelis later showed that it was reasonable to build tables that reached 9-extra bits of accuracy [START_REF] Gal | An accurate elementary mathematical library for the IEEE floating point standard[END_REF]. This means that each tabulated value can be stored with p-bit precision while having a virtual precision of p + 9 bits, which makes the maximum relative error of each value shrink down to 2 -10 ulp.

In 2005, an improvement to the computation of Gal's accurate tables was suggested by Stehlé and Zimmermann [START_REF] Stehlé | Gal's Accurate Tables Method Revisited[END_REF]. By cleverly correlating Gal's table entries and the search for hard-to-round cases of an elementary function, they reduce the relative error upper bound to 2 -20 ulp for each tabulated value. For the trigonometric functions, which require two simultaneously-bad cases for each entry, their algorithm achieves a heuristic complexity of about 2 n/2 instead of the 2 n complexity of Gal's search based on the random model. This allowed them to build accurate tables for the functions 2 x and sin(x) on the interval [START_REF] Aad | Measurements of the Higgs boson production and decay rates and coupling strengths using pp collision data at √ s = 7 and 8 TeV in the ATLAS experiment[END_REF][START_REF] Aad | Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC[END_REF].

Other methods of tabulations have been developed, especially for hardware implementations. For example, we can cite the "(M, p, k)friendly points" method for the trigonometric functions [START_REF] Brisebarre | M, p, k)-Friendly Points: A Table-Based Method for Trigonometric Function Evaluation[END_REF][START_REF] Wang | k)-Friendly Points: A Table-Based Method to Evaluate Trigonometric Functions[END_REF], or bipartite tables [START_REF] Michael | Approximating elementary functions with symmetric bipartite tables[END_REF], and the generalized multipartite tables [START_REF] De | Multipartite table methods[END_REF] for any elementary function. More detail about the "(M, p, k)-friendly points" are given in Chapter 2. Table-based methods can provide a first rough approximation, but they quickly require exponential memory if one wants more accurate values for any input. Hence, combining small tables with efficient numerical approximations has been a design choice for many elementary functions for a long time.

Power series are among the first general numerical methods to approximate most elementary functions [123, ch.11]. The most common approximation functions are polynomials, but sometimes rational approximations (i.e. a quotient of two polynomials) are preferred. Perhaps the most famous rational approximations are Padé approximants, as they often give better approximations than Taylor series [START_REF] Baker | Encyclopedia of Mathematics and its Applications[END_REF].

However, polynomial approximations are often preferred in libms as they do not involve any division, which is often an expensive operation. There are various kinds of polynomial approximations, among which: taylor series They interpolate the function to approximate at a single point of multiplicity n + 1, and have the advantage of conserving parity properties. However, their numerical accuracy deteriorates the closer to the boundaries of the approximation interval [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. These approximations are commonly used for elementary functions implementations. Taylor series are generally chosen for their simplicity and parity-saving property, which then allows to halve the number of polynomial coefficients. Minimax polynomials are often chosen for their optimality in terms of accuracy, but Chebyshev are sometimes preferred because Remez' algorithm is slower than a simple Chebyshev interpolation for large degrees. As said above, there exist many other kinds of polynomial approximations. We may also cite Bernstein polynomials [START_REF] Steffens | The history of approximation theory: from Euler to Bernstein[END_REF], pervasively used for Bézier curves in computer graphics, although they were not used in this thesis.

chebyshev
When it comes to evaluating a polynomial, Horner's rule is "one of the first things a novice programmer is usually taught" [90, p. 486]. Actually, it is a rather sequential way of evaluating a polynomial p(x). Suppose that p is of degree n and that its coefficients are all nonzero. Then Horner's rule gives the following polynomial evaluation scheme:

p(x) = ((p n x + p n-1 ) x + • • • ) x + p 0 .
Because of its elegant simplicity, but also thanks to its numerical accuracy and many other properties [START_REF] Boldo | A Simple Test Qualifying the Accuracy of Horner's Rule for Polynomials[END_REF] Horner's rule is pervasively used for elementary function evaluations.

On hardware architectures that benefit from Instruction-Level Parallelism (ILP) through the use of pipelining, Horner's scheme can suffer from performance issues because each operation requires the intermediate result that is produced immediately before, hence causing pipeline stalls [START_REF] Hennessy | Computer Architecture. A Quantitative Approach[END_REF]. Therefore, various polynomial evaluation schemes have been proposed for univariate [START_REF] Estrin | Organization of Computer Systems: The Fixed Plus Variable Structure Computer[END_REF][START_REF] Eve | The evaluation of polynomials[END_REF][START_REF] Knuth | The Art of Computer Programming[END_REF][START_REF] Kiyoshi | On the Parallel Evaluation of Polynomials[END_REF][START_REF] Munro | Optimal algorithms for parallel polynomial evaluation[END_REF][START_REF] Ya | Methods of computing values of polynomials[END_REF][START_REF] Paterson | On the Number of Nonscalar Multiplications Necessary to Evaluate Polynomials[END_REF], bivariate [START_REF] Mouilleron | Automatic Generation of Fast and Certified Code for Polynomial Evaluation[END_REF][START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors -Polynomial evaluation-based algorithms and certified code generation[END_REF], or multivariate polynomials [START_REF] Boehm | On de Casteljau's algorithm[END_REF][START_REF] Schumaker | Efficient evaluation of multivariate polynomials[END_REF], to optimize their performance.

As we will see in more detail in Chapter 3, it can sometimes be preferable to automatically generate efficient polynomial evaluation schemes, although the exponential search-space makes it still a complex and open research problem [START_REF] Knuth | The Art of Computer Programming[END_REF][START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors -Polynomial evaluation-based algorithms and certified code generation[END_REF].

Finally, let us also cite the rather famous CORDIC algorithm [START_REF] Pramod | 50 years of CORDIC: Algorithms, architectures, and applications[END_REF][START_REF] Palmer | Numeric data processor[END_REF][START_REF] Volder | The CORDIC trigonometric computing technique[END_REF]. COordinate Rotation DIgital Computer was designed and extended to compute trigonometric and hyperbolic functions using simple arithmetic operations (namely, addition, subtraction and bit shift) and table lookups. It was developed in 1959 for hardware that did not have any multiplier, and is still often used on pocket calculators and reconfigurable hardware (FPGAs). On highly-constrained hardware, with little memory dedicated to numerical code, CORDIC gives a good ratio between provided functions and code size, thanks to its versatility [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF][START_REF] Nave | Implementation of Transcendental Functions on a Numerics Processor[END_REF][START_REF] Tak | Table-Lookup Algorithms for Elementary Functions and Their Error Analysis[END_REF].

All the above methods alone have quite limited achievable accuracy for the large input ranges offered by IEEE 754 floating-point numbers. This is why they are often combined with each other or with range reductions, which we briefly present in the next section.

Range Reductions

Polynomial or rational approximations are usually designed to be accurate on small intervals. This way, the degree of the involved polynomials can generally stay low enough for the evaluation to be fast. Additionally, table lookups should not hold more than a few thousands values for performance reasons. Therefore, whenever a function is to be efficiently approximated by a table-lookup or a numerical approximation (or a combination of both), range reductions are often performed. Range reductions allow to shrink the interval on which these approximations take place, so that they can stay efficient. They generally involve mathematical properties of the function

x ∈ F Range reduction Mathematical properties Table-lookups Numerical approximations Reconstruction f(x) ∈ F Figure 1.
3 -Overview of a typical evaluation scheme for an elementary function x → f(x). Note that it results in f(x), and not in f(x), as it is an approximation in finite-precision arithmetic.

that is approximated, so that one or more (simpler) functions have to be approximated on small intervals. The final result is generally computed by a reconstruction step, which roughly "reverts" the range reduction, by combining all intermediate results from the range reduction and from the numerical approximations. Then, a typical evaluation scheme combining these four steps can be sketched as in Figure 1.3. The mathematical properties used for range reduction generally include:

parity If a function defined on R is even or odd, then it is sufficient to evaluate it on R 0 and to change the sign of the result when needed.

periodicity If a function f defined on R is such that

∃C ∈ R, ∀k ∈ Z, ∀x ∈ R, f(x + k • C) = f(x),
then it is sufficient to evaluate it on [-C/2, C/2] after performing an additive reduction by subtracting k • C from x. This reduction step must be carefully done especially if C cannot be stored exactly in memory [START_REF] Brisebarre | A New Range-Reduction Algorithm[END_REF][START_REF] Cody | Software manual for the elementary functions[END_REF][START_REF] Cody | Implementation and Testing of Function Software[END_REF][START_REF] Payne | Radian Reduction for Trigonometric Functions[END_REF].

symmetry If two functions f and g defined on R are such that

∃C ∈ R, ∀x ∈ R, f(x + C) = g(x),
then g(x -C) can be evaluated instead of f(x) if it is advantageous. This can allow for further range reductions, such as for the trigonometric functions.

sinh(x) = e x -e -x 2 ; sinh(n log 2) = 2 n-1 -2 -n-1 ; cosh 2 x -sinh 2 x = 1.
The logarithm function, for example, can be mathematically transformed into log(x) = log(m • 2 e ) = log(m) + e log(2), so as to only have to approximate the log function on [1, 2) for normal numbers [START_REF] Markstein | IA-64 and elementary functions: speed and precision[END_REF].

In this section, we have seen that in the context of high performance libms, the whole evaluation scheme of an elementary function is often split into several steps: range reductions, table-lookups, numerical approximations and the reconstruction of the final result. As compute-intensive applications using elementary functions (such as the CERN CMS experiment software) often present Data-Level Parallelism (DLP), they try to exploit as much parallelism offered by computers as possible. At the hardware level, exploiting DLP can be achieved with Vector, SIMD and GPU architectures [START_REF] Fog | The microarchitecture of Intel, AMD and VIA CPUs: An optimization guide for assembly programmers and compiler makers[END_REF][START_REF] Hennessy | Computer Architecture. A Quantitative Approach[END_REF]. Hence, a current challenge is to develop both efficient and accurate vectorized libms, which is the topic of the next section.

SIMD vectorization of elementary functions

Exhibiting arithmetic parallelism is easy in many applications such as image processing, finite element analysis or 3D graphics. As Moore's law accurately predicted that the number of transistors per area unit would double each 18 months [START_REF] Moore | Lithography and the future of Moore's law[END_REF], hardware designers really began to take advantage of these extra transistors by supporting arithmetic parallelism in the 1990s. Such support was partly enabled by the democratization of Single Instruction Multiple Data (SIMD) architectures [START_REF] Oberman | AMD 3DNow! technology: architecture and implementations[END_REF][START_REF] Peleg | MMX Technology Extension to the Intel Architecture[END_REF].

SIMD architectures provide arithmetic and logic operators that are atomically applied to multiple data. Common SIMD architectures such as Intel's or AMD's are implemented as vector extensions to the Instruction Set Architecture (ISA), and are also known as SIMD Within A Register (SWAR) units. These units operate on vector registers whose length currently ranges from 64 to 512 bits, as can be seen in Table 1.3. For instance, Intel's Advanced Vector Instructions (AVX) and AVX2 operate on 256-bit registers that can be split into e.g. These extensions can increase the throughput of arithmetic-intensive programs, but they usually do not implement vectorized elementary functions, i.e. elementary functions that would return a SIMD vector of outputs from a vector of inputs. However, manufacturers or independent developers have distributed mathematical libraries that provides such functions. They naturally rely on SIMD instructions provided by the hardware, or on the auto-vectorizing capabilities of modern compilers. Most of them are implemented in the C or C++ languages. The commercial Intel SVML library 1 provides intrinsic support for vectorized mathematical functions, but it is not freely available. Intel's Math Kernel Libraries provides vectorized mathematical functions optimized for large arrays of inputs. In this library, each function is implemented in three different flavors: high accuracy (HA), low accuracy (LA), and enhanced performance (EP). The maximum measured absolute error of the various routines ranges from 0.5 ulp for basic functions in their HA flavor to several thousand ulps -i.e. without any significant digit -for some EP routines. 2 Other, non-commercial implementations have been proposed. In 2015, Intel contributed the GNU Libmvec to glibc 2.22 [START_REF] Senkevich | Glibc Wiki[END_REF]. GNU Libmvec claims reasonable testing to pass 4-ulp maximum relative error, on each function domain and only in rounding-to-nearest mode. The latest version (from glibc 2.26) is only optimized for performance, although HA and LA flavors have been mentioned for future developments. It is written in x86 assembly though, which reduces its portability to other architectures and makes it more difficult to maintain and extend. For now, the functions may or may not handle special cases, and scalar callbacks to the original GNU libm functions may be used, which can lower the average performance. AMD also contributed the proprietary but free ACML_MV (no longer available) and the AMD LibM. 3 AMD ACML_MV was announced returning re- sults with a maximum error of at most 1 ulp, but with unpredictable behavior for subnormal inputs. The CERN's VDT is another open source effort aiming at performance [START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF]. It is written in C++ and relies on modern compilers autovectorizing abilities to generate vectorized routines for large arrays. In terms of accuracy, the logarithm function is announced with at most two bits differing from GNU libm results in single and double precision, eight bits being the overall maximum difference for the arccos function in double precision.

Influenced by the Cephes library,1 like the VDT library, Pommier's SSE2 and Garberoglio's AVX3 versions of "Mathfun" headers offer SIMD routines for the sin, cos, exp, and log functions in single precision. Like EP implementations, Mathfun does not claim high accuracy. However the authors write that the trigonometric functions are measured to be faithful on the range [-8192, 8192]. Their performance benchmarks show an average reciprocal throughput close to 30 cycles per element (CPE) for the SSE version, while the AVX version is said to be close to a 1.7x speedup compared to SSE routines.

Two recent projects are Yeppp! [START_REF] Dukhan | Methods for High-Throughput Computation of Elementary Functions[END_REF], and SLEEF [START_REF] Shibata | Efficient evaluation methods of elementary functions suitable for SIMD computation[END_REF]. Yeppp! targets EP implementations by relying on the PeachPy assembly kernel generator [START_REF] Dukhan | PeachPy: A Python Framework for Developing High-Performance Assembly Kernels[END_REF][START_REF] Dukhan | PeachPy Meets Opcodes: Direct Machine Code Generation from Python[END_REF], while the SIMD Library for Evaluating Elementary Functions (SLEEF) provides HA and LA routines [START_REF] Shibata | Efficient evaluation methods of elementary functions suitable for SIMD computation[END_REF].

In 2016, Lauter proposed an open-source auto-vectorizable mathematical library written in high level scalar C [START_REF] Quirin | A new open-source SIMD vector libm fully implemented with high-level scalar C[END_REF]. Its main feature is probably to be system and hardware independent, since it does not rely on any other libraries and only makes use of standard C constructs. This makes this library available on any architecture for which an auto-vectorizing backend is implemented in the compiling toolchain. The library was developed with both binary32 and binary64 formats in mind, but until now efforts have been concentrated on the binary64 routines only. It currently provides nine elementary functions such as trigonometric functions, exponential and logarithm, and claims ongoing work for functions harder to vectorize such as the power function (x, y) → x y or the bivariate arctangent (x, y) → arctan(y/x). Each implementation targets enhanced performance at the expense of accuracy: a maximum error of 8 ulps was measured for the whole set of functions. The implementation, which is self-sufficient as it requires no other library, relies on compiler autovectorizers like the CERN's VDT. For this reason, conditional branches and table-lookups are entirely avoided.

Note that none of these mathematical libraries provide a guaranteed accuracy, since, to our knowledge, documented accuracies are obtained by non-exhaustive testing. Furthermore, writing high performance libms is a tedious, error-prone task that often leads to monolithic, general-purpose libms. As for any monolithic library, one can think that general-purpose libms are to tailored software what COTS1 circuits are to full-custom ASICs2 in hardware design: they usually meet the requirements, but something else could fulfill specific needs much better. As the knowledge about floating-point arithmetics and elementary functions grows, the current challenge is now to automatize the generation of code for libms, or provide tools that help developers in this time-consuming task. This is the topic of the next section.

code generation for elementary functions

The design and implementation of a whole libm usually involve high development costs [START_REF] Kupriianova | Towards a modern floating-point environment[END_REF]. Furthermore, real-world applications using mathematical functions are so diverse that they often need specific levels of accuracy over specific input intervals with respect to their targeted performance (e.g. throughput or latency) [START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF]. For instance, different research fields such as earthquake engineering or acoustics may use the same trigonometric functions but with very different input intervals. The first may use angles in the interval [-π, π] as in [140, p. 290], while the second might use much wider intervals such as 0, 10 4 as in [START_REF] Dai | Compensation of Loudspeaker Nonlinearity in Acoustic Echo Cancellation Using Raised-Cosine Function[END_REF]. Such diverse needs push for tailored implementations of each function [START_REF] De | Matériel et logiciel pour l'évaluation de fonctions numériques: précision, performance et validation[END_REF]. As this possibly means an exponential number of different implementations, recent efforts have been put on automatic code generation [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF][START_REF] Brunie | Code Generators for Mathematical Functions[END_REF][START_REF] De | Code generation for argument filtering and argument reduction in elementary functions[END_REF][START_REF] Florent De Dinechin | Automatic generation of polynomial-based hardware architectures for function evaluation[END_REF][START_REF] Kupriianova | Towards a modern floating-point environment[END_REF][START_REF] Quirin | Arrondi correct de fonctions mathématiques -Fonctions univariées et bivariées, certification et automatisation[END_REF][START_REF] Lauter | Semi-Automatic Floating-Point Implementation of Special Functions[END_REF][START_REF] Revy | Automated Design of Floating-Point Logarithm Functions on Integer Processors[END_REF][START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors -Polynomial evaluation-based algorithms and certified code generation[END_REF] to let the machines handle the exponential development costs.

The results of these efforts have been to reduce development costs down to a few man-months for a whole custom libm. Many software tools that can generate or help at generating code for elementary function evaluation have been developed over the recent years. Now, these efforts have been aiming at bringing development costs further down to a few man-hours. We present here such tools that were used for this thesis.

Tools geared towards code generation

Sollya: an environment for the development of numerical codes Sollya 3 is an extensible scripting language and a C library that can perform symbolic computations on many elementary functions, as well as compute different polynomial approximations. Among its features, Sollya provides the following algorithms:

• guessdegree tries to guess the minimal polynomial degree that is required to approximate a function on some interval within an error bound;

• taylor, remez and fpminimax try to generate Taylor, Remez, or FP-minimax polynomial approximations, for elementary or special functions [START_REF] Brisebarre | Efficient polynomial L ∞ -approximations[END_REF][START_REF] Sylvain Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF][START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF];

• implementpoly can generate C code that implements a polynomial object in optimized double, double-binary64 or triple-binary64 precision following Horner's scheme [START_REF] De | Optimizing polynomials for floating-point implementation[END_REF].

A code generator for polynomial evaluation: CGPE

Although Sollya provides implementpoly to generate optimized-precision C codes for polynomial evaluation, only the Horner's scheme is considered. As we will see in Chapter 3, one might want to achieve higher computing throughput or lower latency by using different evaluation schemes. This is the purpose of CGPE. CGPE 1 is an interactive software and a C++ library that can generate efficient code for polynomial evaluations on multiple targets [START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors -Polynomial evaluation-based algorithms and certified code generation[END_REF]. It is able to generate polynomial schemes for bivariate polynomials and return the ones that exhibit the most Instruction Level Parallelism (ILP), i.e. schemes in which many arithmetic operations are independent, hence that could be run in parallel. Given bivariate polynomial coefficients, input intervals and an error bound, CGPE can also try to generate an efficient polynomial evaluation scheme that satisfies the error bound for the polynomial on its input intervals. In this case, CGPE can formally verify such error bound using Gappa, which is briefly presented in the next section.

Formalization and code certification

Finite precision implementations usually document their approximation errors or provide guarantees on the relative error on the outputs compared to the mathematical results. Such guarantees can be obtained by a tight analysis of evaluation and approximation errors. This analysis is tedious and error-prone, but can be automated with tools such as Gappa [START_REF] Daumas | Certification of Bounds on Expressions Involving Rounded Operators[END_REF][START_REF] Florent De Dinechin | Assisted Verification of Elementary Functions Using Gappa[END_REF][START_REF] Florent De Dinechin | Certifying the Floating-Point Implementation of an Elementary Function Using Gappa[END_REF].

Gappa 2 combines interval arithmetic and logical mathematical rewritings to automate and verify intelligent forward error analyses. It can manipulate many kinds of real numbers, including binary32, binary64 or binary128 floating-point numbers [START_REF] Melquiond | De l'arithmétique d'intervalles à la certification de programmes[END_REF]. For example, the following Gappa script is able to determine a tight interval in which lies the binary32 square of a binary32 approximation of a real number in [-0.1, 0.1]. which means that x * x is in the interval [0, 10737419b-30] in dyadic notation, that is, in [0, 0x1.47ae16p-7] in hexadecimal notation.

The MetaLibm ANR project

MetaLibm is the name of multiple automatic code generation projects1 that eventually merged to form the MetaLibm ANR project,2 in which this thesis takes place. All of these projects have been aiming at generating source code for the evaluation of elementary functions or for digital filters, but they differ in their respective approaches.

Metalibm-lutetia, a black-box generator

One of them is the Metalibm-lutetia project, which adopts a top-down approach. It was designed for wide adoption by developers with or without numerical expertise. Indeed, it is a framework that tries to generate code for general mathematical expressions, with I/O precision, range and target accuracy specifications.

First, it tries to automatically detect periodicities or symmetries that can allow for range reductions. Then, for the remaining interval on which the function must be approximated, it applies an algorithm of domain splitting: the interval is split into as many subintervals as needed so that the function can be approximated by different numerical approximations on each of these intervals. For polynomial approximations, a halting condition for splitting can be a maximum value for the polynomial degree. Kupriianova and Lauter developed a split-then-merge algorithm that uses such a condition to perform domain splitting and then merge as many subintervals as possible, in order to optimize the generated code [START_REF] Kupriianova | A domain splitting algorithm for the mathematical functions code generator[END_REF].

Metalibm-lugdunum, a libm-oriented generator

The second project is Metalibm-lugdunum. 3 This framework enables to describe the implementation of a function using a meta-language in a Python syntax. More particularly the descriptions consist in abstracting evaluation schemes as a set of high-level building blocks. Then, the MetaLibm backend can optimize and translate the whole description for different target architectures (C11, x86-intrinsics C11, OpenCL, VHDL, . . . ) and numerical formats (floating-point, fixedpoint, . . . ) [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF][START_REF] Brunie | Code Generators for Mathematical Functions[END_REF].

As an example, the following listing describes the multiplication x • x in this meta-language, where x is the input variable of our program. y = Multiplication(x, x, tag="y", precision = self.precision)

Then, if the generic C11 architecture is targeted, Metalibm will simply generate the statement y = x * x, while if x86 architectures with AVX2 are targeted, it will produce the following piece of C code, where vec _ x is the vector equivalent of x. Metalibm-lugdunum is more oriented towards libm developers as it follows a bottom-up approach. However, it can optimize generated implementations from abstract evaluation schemes for various goals, such as (among others) target accuracy, input interval, vector size, or compliance to the IEEE 754 standard. Therefore, it is also designed to be used as a toolbox, e.g. by theoretical physicists needing a tailor-made, efficient OpenCL implementation of the logarithm or the exponential functions.

As abstract evaluation schemes can be target-agnostic, Metalibm pushes for factoring shared algorithms across architectures, different I/O precisions, or target accuracies. Such design aims at shrinking the development time for mathematical libraries. Facilities such as automatic testing for performance and accuracy, and formal proof generation are included, which represents an additional time-saving opportunity for meta-libm developers and final users.

Metalibm-tricassium, a filter-oriented generator

Digital signal processing (DSP) numerical algorithms have much in common with polynomial and rational functions. Therefore, the Meta-Libm project also aims at generating efficient and accurate source codes for the implementation of digital filters on various architectures [START_REF] Lopez | Optimal implementation of linear filters in fixed-point arithmetic[END_REF][START_REF] Volkova | Towards reliable implementation of digital filters[END_REF].

Generative programming for numerical algorithms

As said by Rompf et al., there has been an ever-growing appeal to "go meta" [START_REF] Rompf | Go Meta! A Case for Generative Programming and DSLs in Performance Critical Systems[END_REF]. For numerical programs, the number of domain specific languages (DSL) and source-to-source compilers is skyrocketing. They now target various optimizations goals for improved, automatically generated or transformed, numerical programs [20, 33, 35, 36, 48, 49, 54, 58-60, 80, 81, 83, 84, 112, 116, 117, 126, 131, 134, 145, 161, 162, 165, 166]. This list is far from being exhaustive and could probably be updated every month. We will just present some of the goals that such "meta-optimizers" target. Without going into much detail, some of them aim at improving and certifying the accuracy of floating-point or fixed-point programs such as Sardana/Salsa [START_REF] Damouche | Improving the numerical accuracy of programs by automatic transformation[END_REF][START_REF] Ioualalen | Synthesis of arithmetic expressions for the fixed-point arithmetic: The Sardana approach[END_REF][START_REF] Ioualalen | Synthesizing accurate floating-point formulas[END_REF] or Rosa/Daisy [START_REF] Darulova | Towards a Compiler for Reals[END_REF][START_REF] Darulova | Synthesis of fixed-point programs[END_REF][START_REF] Izycheva | On Sound Relative Error Bounds for Floating-Point Arithmetic[END_REF]. Other explore better tradeoffs between efficiency and accuracy for hardware [START_REF] Istoan | High-performance Coarse Operators for FPGAbased Computing[END_REF][START_REF] Pasca | High-performance floating-point computing on reconfigurable circuits[END_REF][START_REF] Thomas | A General-Purpose Method for Faithfully Rounded Floating-Point Function Approximation in FPGAs[END_REF], such as FloPoCo [START_REF] Florent De Dinechin | Automatic generation of polynomial-based hardware architectures for function evaluation[END_REF][START_REF] Florent De Dinechin | Generating high-performance custom floating-point pipelines[END_REF][START_REF] De | Designing Custom Arithmetic Data Paths with FloPoCo[END_REF], software (PeachPy [START_REF] Dukhan | PeachPy: A Python Framework for Developing High-Performance Assembly Kernels[END_REF][START_REF] Dukhan | PeachPy Meets Opcodes: Direct Machine Code Generation from Python[END_REF][START_REF] Dukhan | Methods for High-Throughput Computation of Elementary Functions[END_REF]), or both [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF].

E X A C T L O O K U P TA B L E S F O R E L E M E N TA R Y F U N C T I O N S
"It is a bit uncertain who first realized that [Primitive Pythagorean Triples] had a beautiful family tree, but a 1934 paper in Swedish could be the earliest (Berggren 1934 [10])."

Bernhart and Price [START_REF] Bernhart | Pythagoras' garden, revisited[END_REF].

As we have seen in Chapter 1, elementary mathematical functions are pervasively used in many applications such as electronic calculators, computer simulations, or critical embedded systems. However, their evaluation is usually an approximation. To reach high accuracies, they usually make use of mathematical properties, precomputed tabulated values and accurate polynomial approximations [START_REF] Defour | Elementary functions : algorithms and efficient implementation for correct rounding for the double precision[END_REF]Ch.2]. Each step of the evaluation process generally combines error of approximation (e.g. replacing an elementary function by a polynomial) and error of evaluation on finite-precision floating-point arithmetic (e.g. the evaluation of such a polynomial). Whenever they are used, tabulated values usually embed rounding error inherent to the transcendence of most elementary functions [START_REF] Markstein | IA-64 and elementary functions: speed and precision[END_REF]. Several interesting techniques have been developed to reduce this rounding error as much as possible, yet, to our knowledge, none has ever managed to remove it entirely.

In this chapter, we present a general method that enables to build error-free tabulated values that is worthwhile whenever at least two terms are tabulated in each table row. Our technique transfers all potential rounding errors on tabulated values to one small corrective term, which allows to asymptotically save twice as many bits as stateof-the-art techniques. For the trigonometric and hyperbolic functions sin, cos, sinh, and cosh, we show that Pythagorean triples can lead to such tables in reasonable time and memory usage. When targeting correct rounding in double precision for the same functions, we also demonstrate that this method might save memory accesses and floating-point operations by up to 29% and 42%, respectively, during the reconstruction step.

exact lookup tables for elementary functions

introduction

The algorithms developed for the evaluation of elementary functions, such as the logarithm, the exponential, or the trigonometric and hyperbolic functions, can be classified into at least two categories. The first category concerns algorithms based on small and custom operators combined with tabulated values that target small accuracies, typically less than 30 bits [START_REF] Defour | A new scheme for table-based evaluation of functions[END_REF][START_REF] De | Multipartite table methods[END_REF][START_REF] Debjit | Faithful Bipartite ROM Reciprocal Tables[END_REF]. The second category concerns algorithms that usually target single or double precision (24 or 53 bits of precision) with an implementation on general processors that relies on the available hardware units [START_REF] Daramy-Loirat | A library of correctly rounded elementary functions in double-precision[END_REF][START_REF] Ziv | Fast Evaluation of Elementary Mathematical Functions with Correctly Rounded Last Bit[END_REF].

Implementations of those functions that target correct rounding are usually divided into two or more phases [START_REF] Muller | Elementary Functions[END_REF]. A quick phase is first performed: it is based on a fast approximation that provides a few more bits than the targeted format, which makes correct rounding possible most of the time at a reasonable cost. When correct rounding is not possible, a slower accurate phase is executed. The quick phase uses operations with a precision slightly greater than the targeted precision, while the accurate phase is based on a larger precision. For instance, in the correctly-rounded library CR-Libm, the quick phase for the sine and cosine functions in double precision targets 66 bits while the accurate phase corresponds to 200 bits [34, § 9]. Then, in order to guarantee that an implementation actually computes correctly-rounded results, a proof of correctness must be built. This proof is based on the mandatory number of bits required to ensure correct rounding, which is linked with the search for the worst cases for the TMD [START_REF] Lefèvre | Hardest-to-Round Cases[END_REF].

The design of such correctly-rounded implementations is a real challenge, as it requires to control and limit every source of numerical error [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF]. Indeed these implementations involve various steps, such as range reduction, polynomial evaluation, or reconstruction, during all of which errors may occur. Since those errors accumulate and propagate up to the final result, any solution or algorithm that reduces them can have an impact on the simplicity of the proof of correctness and the performance of the implementation.

Overview of this chapter

This chapter presents works on the trigonometric functions sin and cos that were published and presented in [START_REF] De | Range reduction based on Pythagorean triples for trigonometric function evaluation[END_REF], their extension to the hyperbolic functions sinh and cosh that were published in [START_REF] De | Exact Lookup Tables for the Evaluation of Trigonometric and Hyperbolic Functions[END_REF], as well as a faster heuristic to build trigonometric exact lookup tables that has, to our knowledge, never been published before. Our main contributions are twofold:

1. a general method that eliminates rounding errors from tabulated values and transfers them into a single corrective term to be added to the reduced argument, which, in most cases, already contains rounding errors;

q ≡ 0 mod 4 q ≡ 1 mod 4 q ≡ 2 mod 4 q ≡ 3 mod 4 sign(x) = 1 sign(x) = -1 Figure 2.
1 -Graphical representation of the first range reduction for the trigonometric functions.

2. and a formally justified algorithm that shows how to apply this method to trigonometric and hyperbolic functions.

In Section 2.2, we recall classic evaluation schemes for the trigonometric and hyperbolic functions. Section 2.3 details the properties of the proposed tables and demonstrates how it is possible to remove two sources of error involved in this step. Section 2.4 presents our theoretical approach to build exact lookup tables for trigonometric and hyperbolic functions using Pythagorean triples. Then, Section 2.5 presents some experimental results that show that we can precompute exact lookup tables up to 12 indexing bits reasonably fast with the help of suitable heuristics. A quantitative comparison with classic approaches is made in Section 2.6. They show that our method can lower the table sizes and the number of floating-point operations performed during the reconstruction step. Finally, an example with a toy table is also presented in Section 2.7. We conclude and discuss perspectives in Section 2.8.

table-based evaluation schemes 2.2.1 Usual framework

In order to better understand the challenges that we face when designing evaluation schemes for elementary functions, let us detail each step of a classic function evaluation process with a focus on methods based on table lookups through the example of the trigonometric and hyperbolic sine and cosine. For this purpose, let y be a machine-representable floating-point number, taken as the input to the considered functions. The design of these evaluation schemes is usually very similar to the four-step process that we have seen in Chapter 1 (p.22):

1. A first range reduction, based on mathematical properties, narrows the domain of the function to a smaller one. For the trigonometric functions, properties of periodicity and symmetry lead to evaluating f q ∈ {± sin, ± cos} at input

|x| = |y -q • π/2| ∈ [0, π/4],
where q ∈ Z and the choice of f q depends on the function to evaluate, q mod 4, and sign(x) [START_REF] Muller | Elementary Functions[END_REF]. Figure 2.1 sums up how q mod 4 and sign(x) may change the function to evaluate.

For the hyperbolic functions, let

|x| = |y -q • ln(2)| ∈ [0, ln(2)/2] ,
with q ∈ Z. Addition formulas can then be used together with the analytic expressions for the hyperbolic functions involving the exponential [START_REF] Markstein | IA-64 and elementary functions: speed and precision[END_REF], which gives, for the sine:

sinh(y) = 2 q-1 -2 -q-1 • cosh |x| + (-1) sign x • 2 q-1 + 2 -q-1 • sinh |x| .

2.

A second range reduction, based on tabulated values, further reduces the range on which polynomial evaluations are to be performed. The argument x is split into two parts, x h and x , such that:

x = x h + x with |x | 2 -p-1 . (2.1)
The term x h is the p-bit value of the form

x h = i • 2 -p with i = x • 2 p , ( 2.2)
where a denotes the rounding of a floating-point number a to the nearest integral value. 1 The integer i is used to address a table of n = κ • 2 p-1 + 1 rows, with κ ∈ {π/2, ln(2)}. This table holds precomputed values of either trigonometric or hyperbolic sines and cosines of x h , which we indifferently name S h and C h .

3.

Meanwhile, polynomial approximations of the trigonometric or hyperbolic sine and cosine on the interval [-2 -p-1 , 2 -p-1 ], denoted by P S and P C , are evaluated at input x .

4. Finally, a reconstruction step allows to compute the final result using the precomputed values retrieved at step 2 and the computed values from step 3. For the trigonometric sine, if we assume that q ≡ 0 mod 4 and sign(x) = 1 so that f q = sin, one has to perform the following reconstruction:

sin(y) = S h • P C (x ) + C h • P S (x ),
1 Tie cases are generally rounded to the nearest even integral value, e.g. 2.5 = 2.

while the reconstruction for the hyperbolic sine is:

sinh(y) = 2 q-1 -2 -q-1 • (C h • P C (x ) + S h • P S (x )) + (-1) sign x • 2 q-1 + 2 -q-1 • (S h • P C (x ) + C h • P S (x )) .
Satisfactory solutions already exist to address the first range reduction [START_REF] Brisebarre | A New Range-Reduction Algorithm[END_REF][START_REF] Daumas | Modular Range Reduction: a New Algorithm for Fast and Accurate Computation of the Elementary Functions[END_REF][START_REF] Payne | Radian Reduction for Trigonometric Functions[END_REF], the generation and evaluation of accurate and efficient polynomial evaluation schemes [START_REF] Brisebarre | Efficient polynomial L ∞ -approximations[END_REF][START_REF] De | Optimizing polynomials for floating-point implementation[END_REF][START_REF] Mouilleron | Automatic Generation of Fast and Certified Code for Polynomial Evaluation[END_REF], and the reconstruction step. The interested reader can find more details in [START_REF] Muller | Elementary Functions[END_REF].

All these steps generally add errors of approximation and/or evaluation. Here, we address the second range reduction based on tabulated values for the trigonometric and hyperbolic sine and cosine. At this step, each tabulated value embeds a rounding error. Our objective is to remove the error on these values, and consequently transfer the error due to range reduction into the reduced argument used in the polynomial evaluations. Precisely the proposed method relies on finding precomputed values with remarkable properties that simplify and accelerate the evaluation of these functions. These properties are threefold:

1. Each pair of values holds the exact images of a reference argument under the functions, that is, without any rounding error.

For this purpose, we require these values to be rational numbers;

2. The numerator and denominator of each rational value are exactly representable in a representation format available in hardware, e.g., binary64;

3. And finally, these rational values all share the same denominator. This enables to incorporate the division by this denominator directly into the polynomial coefficients.

These three properties lead to tabulated values that are exact and representable on single machine words. For the trigonometric and hyperbolic functions, Pythagorean triples give rational numbers that satisfy such properties.

Focus on table-lookup range reductions

The second range reduction uses fast but inaccurate approximations. In practice it is often implemented using table lookups. 1 This section presents three solutions that address this step with the sine function to illustrate them. However, these methods are general and also apply to the trigonometric cosine as well as to the hyperbolic functions.

C h ≈ cos(x h ) S h ≈ sin(x h ) x x h x

Tang's tables

Tang proposed a general method to implement elementary functions that relies on hardware-tabulated values [START_REF] Tak | Table-Lookup Algorithms for Elementary Functions and Their Error Analysis[END_REF]. Given the reduced argument x as in Equation (2.1), Tang's method uses the upper part

x h to retrieve two tabulated values S h and C h that are good approximations of sin(x h ) and cos(x h ), respectively, rounded to the destination format. Figure 2.2 gives a visual representation of the second range reduction with an example of Tang's tables for the trigonometric functions. If an implementation targets correct rounding, then those approximations are generally stored as floating-point expansions.

A floating-point expansion of size n usually consists in representing a given number as the unevaluated sum of n non-overlapping floatingpoint numbers so that the rounding error be reduced compared to a single floating-point number [START_REF] Richard | Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates[END_REF]. If we denote by RN i (x) the rounded value of a number x to the nearest floating-point number of precision i, let us note ε -i its relative rounding error such that |ε -i | 2 -i , then, for the trigonometric functions, we have:

S h = RN 53j (sin(x h )) = sin(x h ) • (1 + ε -53j ) and C h = RN 53j (cos(x h )) = cos(x h ) • (1 + ε -53j ),
where j is the size of the binary64 expansions used to represent S h and C h .

In parallel to the extraction of the values S h and C h , the evaluation of two polynomial approximations P S (x) and P C (x) can be performed. They respectively approximate the sine and cosine functions over the interval covered by x , namely [-2 -p-1 , 2 -p-1 ]. Finally, the result of sin(x) is reconstructed as follows:

sin(x) ≈ S h • P C (x ) + C h • P S (x ).
Tang's method is well suited for hardware implementations on pipelined and out-of-order architectures. It takes advantage of the capability on these architectures to access tabulated values in memory and to perform floating-point computations concurrently. Once the argument x is split into the two parts x h and x , memory units can provide the two tabulated values S h and C h , while floating-point units evaluate the polynomials P S and P C . As the degree of the polynomials decreases when the table size increases, the objective is to find parameters so that the polynomial evaluations take as long as memory accesses, on average [START_REF] Defour | Cache-Optimised Methods for the Evaluation of Elementary Functions[END_REF].

Tang's tables store images of regularly spaced inputs, rounded to the destination format. This rounding error is problematic when seeking a correctly rounded implementation in double precision, since worst cases for trigonometric and hyperbolic functions require more than 118 bits of accuracy [START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF][START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF]. A solution consists in storing values on expansions of size 3 and using costly extended-precision operators [START_REF] Daramy-Loirat | A library of correctly rounded elementary functions in double-precision[END_REF].

Gal's accurate tables

In Tang's method, S h and C h are approximations of sin(x h ) and cos(x h ), respectively. They are rounded according to the format used in the table and the targeted accuracy for the final result. To increase the accuracy of these tabulated values, Gal suggested a method to transfer some of the errors due to rounding over the reduced argument [START_REF] Gal | Computing elementary functions: A new approach for achieving high accuracy and good performance[END_REF]. This consists in introducing small, exact corrective terms on the values x h , hereafter denoted by corr. Figure 2.3 illustrates how such a table can be used for the trigonometric functions. For each input entry x h of the table, one corr term is carefully chosen to ensure that both sin(x h + corr) and cos(x h + corr) are "very close" to floating-point machine numbers. In [START_REF] Gal | An accurate elementary mathematical library for the IEEE floating point standard[END_REF], Gal and Bachelis were able to find corr values such that

S h = RN 53j (sin(x h + corr)) = sin (x h + corr) • 1 + ε -10-53j and C h = RN 53j (cos(x h + corr)) = cos (x h + corr) • 1 + ε -10-53j C h ≈ cos(x h + corr) S h ≈ sin(x h + corr) corr x x h
x for each row of an 8-bit-indexed table. This corresponds to 10 extra bits of accuracy for both tabulated values of sine and cosine compared to Tang's tabulated values, thanks to a small perturbation corr on the values x h . For the trigonometric and hyperbolic functions, such corr values are produced by a pseudo-random sampling with an expected sample size of 2 18 [START_REF] Stehlé | Gal's Accurate Tables Method Revisited[END_REF]. Indeed, using a statistically pseudo-random model, Gal expects these corr values to have about 18 bits of significand each. This method entails storing the corrective terms corr along with the values S h and C h , which will usually make Gal's accurate tables larger than Tang's tables when targeting correct rounding. The value sin(x) is reconstructed as follows:

sin(x) = S h • P C (x -corr) + C h • P S (x -corr),
where P C (x) and P S (x) are polynomials approximating cos(x) and sin(x), respectively, on the interval covered by x -corr.

Gal's solution requires a rather expensive search in order to find one corrective term for each entry x h , within a search space that grows exponentially with the number of extra bits for S h and C h . Stehlé and Zimmermann proposed an improvement based on lattice reduction and worst cases to accelerate the precomputation [START_REF] Stehlé | Gal's Accurate Tables Method Revisited[END_REF]. They were able to increase the accuracy of Gal's tables by 11 extra bits, which translates to 21 extra bits compared to Tang's tables.

Gal's technique virtually increases the accuracy of tabulated val-

ues by a few extra bits, slightly reducing the error bounds in the reconstruction step. However, there are still errors in the reduced argument as well as in tabulated approximations.

Brisebarre et al.'s (M, p, k)-friendly points

In 2012, Brisebarre, Ercegovac, and Muller proposed a new method for trigonometric sine and cosine evaluation with a few table lookups and additions in hardware [START_REF] Brisebarre | M, p, k)-Friendly Points: A Table-Based Method for Trigonometric Function Evaluation[END_REF]. It was improved and implemented by Wang and the same authors in 2014 [START_REF] Muller | Elementary Functions[END_REF][START_REF] Wang | k)-Friendly Points: A Table-Based Method to Evaluate Trigonometric Functions[END_REF]. Their approach consists in tabulating four values a, b, z, and x, defined as:

z = 1/ a 2 + b 2 and x ≈ arctan(b/a),
where a and b are small particular integers. The reconstruction then corresponds to:

sin(x) = b • cos(x -x) + a • sin(x -x) • z.
The values (a, b, z, x) are chosen amongst specific points with integer coordinates (a, b) called (M, p, k)-friendly points. These points are recoded using the canonical recoding, which minimizes the number of nonzero digits in table entries [START_REF] Reitwiesner | Binary Arithmetic[END_REF]. More precisely, a and b must be positive integers less than M, so that the number z = 1/ √ a 2 + b 2 has less than k nonzero bits in the first p bits of its canonical recoding. Therefore, with this method, a and b can be stored exactly, but not z nor 1/z in general. However these requirements make hardware multiplications by a, b, and z much more acceptable, since they are equivalent to just a few additions. Compared to other hardware evaluation schemes, the authors claim that their solution can reduce the area on FPGAs for 24-bit accuracy by 50%. Overall, this is a suitable method for competitive low and medium-precision hardware approximations. But although it could be an interesting question to investigate, this technique does not seem specifically suited for software implementations nor for high precisions.

Furthermore, as far as our knowledge goes, there have not been any similar application of this approach to the evaluation of the hyperbolic sine and cosine. For these functions, we may imagine a similar method using small integers a > b, so that z = 1/ √ a 2 -b 2 makes (a, b) an (M, p, k)-friendly point. By defining x ≈ atanh(b/a), this method might also be suitable for the evaluation of hyperbolic functions. However, as the original authors have not made any statement about such an application, further investigations remain to be done to prove its feasibility.

exact lookup tables for elementary functions

design of an exact lookup table

After the first range reduction, the reduced number x is assumed to be an irrational number. Therefore, it has to be rounded to some finite precision, which means that after the second range reduction, only

x contains a rounding error. Gal's method adds an exact corrective term corr to x that allows to increase the accuracy of transcendental tabulated values. Instead, we suggest not to worry about inexact corrective terms, as long as they make tabulated values exactly representable. This way, the error is solely concentrated in the reduced argument x -corr used during the polynomial evaluations. Now, let us characterize these exactly representable values and then explain how the evaluation scheme can benefit of this property. For this purpose, let f and g be the functions to approximate, κ ∈ {π/2, ln(2)} the additive constant for the first range reduction, and y an input floating-point number. As seen in the introduction, the reduced argument x obtained after the first range reduction is such that x = |y -q • κ| with q ∈ Z and 0 x κ/2.

As y is rational, q is an integer, and κ is irrational, then unless q = 0, x must be irrational, and it has to be rounded to a floating-point number x, with a precision j greater than the targeted format such that: x = x • (1 + ε -j ). We should mention that x is generally represented as a floating-point expansion of size 2 or more to reach an accuracy of at least j bits. As seen in Equation (2.1), the second range reduction splits x into two parts, x h and x , such that x = x h + x and |x | 2 -p-1 .

As shown in Equation (2.2), the value x h is then used to compute an address i in the table T made of n rows. Let us assume that this table T holds exact precomputed values of k • f(x h + corr i ) and

k • g(x h + corr i ),
where k is an integer scale factor that makes both tabulated values integers, and corr i is an irrational corrective term, precomputed for each table entry i, such that |corr i | 2 -p-1 . To build such values, we assume the following properties:

1. The functions f and g are right invertible on the domain [0, (n -

1/2) • 2 -p ].
This allows corrective terms to be determined.

Corrective terms are such that |corr

i | 2 -p-1 .

3.

For each table entry i, the values f(x h + corr i ) and g(x h + corr i ) are rational numbers with the same denominator, that is:

f(x h + corr i ) = η i k i and g(x h + corr i ) = γ i k i with η i , γ i ∈ Z and k i ∈ Z * .
4. Let k = lcm(k 0 , . . . , k n-1 ), i.e. the least common multiple (LCM) of the denominators k i . Let

F i = k • f(x h + corr i ) and G i = k • g(x h + corr i ).

It is clear that the values

F i = k • η i /k i and G i = k • γ i /k i are
integers, but we will also assume that they are representable as one floating-point machine number each, e.g. they fit on 53 bits if using the binary64 format.

With numbers satisfying those properties, the reconstruction step corresponds to

f(x) = G i • P f/k (x -corr i ) + F i • P g/k (x -corr i )
for some i in [0 . . n -1], where

• P f/k (x) and P g/k (x) are polynomial approximations of the functions f(x)/k and g(x)/k, respectively, on the interval

-2 -p-1 -max i (corr i ), 2 -p-1 -min i (corr i ) ,
• and F i , G i , and RN j (corr i ) are gathered from T . The integers F i and G i can be stored without error as one floating-point number each. The third tabulated value RN j (corr i ) is a rounding of the corrective term corr i such that:

x h + corr i ∈ f -1 F i k ∪ g -1 G i k .
It can be seen that the new reduced argument x -corr i covers a wider range than just x . In the worst case, this interval can have its length doubled, making the polynomial approximations potentially more expensive. Actually, those polynomials approximate functions that usually do not variate much on such reduced intervals, so that their degree need rarely be increased to achieve the same precision as for an interval that may be up to twice as narrow. This was tested with some of our tables with Sollya's guessdegree function [START_REF] Sylvain Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF]. In unfortunate cases for which the degree should be increased, rarely more than one additional monomial will be required, which barely adds two floating-point operations (one addition and one multiplication). Furthermore, such additional operations can often benefit from instruction parallelism in the evaluation scheme.

Also note that instead of considering polynomial approximations of f(x) and g(x) directly, we propose to incorporate the division by k into the polynomial coefficients. Therefore, we consider approximations P f/k (x) and P g/k (x) of f(x)/k and g(x)/k, respectively, which avoid the prohibitive cost of the division and the associated rounding error.

application to sine and cosine functions

The proposed lookup table for the second range reduction brings several benefits over existing solutions. However, building such tables of error-free values is not trivial, since the integers F i , G i , and k are not always straightforward, especially for transcendental functions. For the trigonometric and hyperbolic sine and cosine, we rely on some useful results on Pythagorean triples [START_REF] Sierpi | Pythagorean triangles[END_REF]. These objects are briefly described in the first part of this section. Then, we present a method to efficiently generate Pythagorean triples using the Barning-Hall tree. Finally, we detail the selection criteria for the generated Pythagorean triples to build exact lookup tables for the trigonometric and hyperbolic functions.

Pythagorean triples and their properties

Pythagorean triples are a set of mathematical objects from Euclidean geometry which have been known and studied since ancient Babylonia [30, ch. 6]. There exist several definitions of Pythagorean triples that usually differ on authorized values. We choose the following one:

Definition 3. A triple of positive integers (a, b, c) is a Pythagorean triple if and only if a 2 + b 2 = c 2 .
A Pythagorean triple (a, b, c) for which a, b, and c are coprime is called a primitive Pythagorean triple (PPT). In the following, we will refer to the set of PPTs as PPT. Recall that we are looking for several rational numbers which hold exact values for the trigonometric and hyperbolic functions. Yet a PPT and its multiples share the same rational values a/b, b/c, . . . For example, the well known PPT (3, 4, 5) and all its multiples can be associated to the ratio a/b = 3/4. Therefore we can restrict our search to primitive Pythagorean triples only, and apply a scale factor afterwards if needed.

According to the fundamental trigonometric and hyperbolic identities, we have

∀x ∈ R, cos(x) 2 + sin(x) 2 = 1 cosh(x) 2 -sinh(x) 2 = 1. (2.
3) It follows from Definition 3 and Equation (2.3) that all Pythagorean triples can be mapped to rational values of trigonometric or hyperbolic sine and cosine. Indeed, let (a, b, c) be a Pythagorean triple. Without loss of generality, we can assume that b = 0. Then we have: The consequences of these properties are twofold:

a 2 + b 2 = c 2 ⇐⇒ a c 2 + b c 2 = 1 ⇐⇒ c b 2 - a b 2 = 1.
1. A PPT and its multiples share the same angles θ and ϕ.

2. Any Pythagorean triple can be mapped to the sides of a right triangle (possibly degenerated into a segment), whose hypotenuse is the third item c of the triple. This is illustrated in Figure 2 Hence, in the following, the word "hypotenuse" is used to refer to the third item of a Pythagorean triple, while the word "legs" is used to refer to its first and second items. Also, we often treat Pythagorean triples as row or column vectors, which will prove to be a useful way to manipulate these objects.

Construction of primitive Pythagorean triples

The set of primitive Pythagorean triples PPT is countably infinite [START_REF] Freek | Over Pythagorese en bijna-Pythagorese driehoeken en een generatieproces met behulp van unimodulaire matrices [On Pythagorean and quasi-Pythagorean triangles and a generation process with the help of unimodular matrices[END_REF][START_REF] Hall | 232. Genealogy of Pythagorean Triads[END_REF][START_REF] Price | The Pythagorean Tree: A New Species[END_REF]. Figure 2.5 represents the subset of all PPTs with a hypotenuse c < 2 12 . It shows that PPTs rapidly cover a wide range of angles over [0, π/2] as c increases. But, as our exact lookup tables need one triple per row, one may ask if there will always be at least one PPT that matches each row, no matter how narrow the entry interval of a row can be.

Actually, a result due to Shiu states that the set of trigonometric angles covered by Pythagorean triples is dense in [0, π/4]. This is formalized by Theorem 1.

Theorem 1 (Shiu 1983 [154]). ∀θ ∈ [0, π/4], ∀δ > 0, there exists a primitive Pythagorean triple with a corresponding trigonometric angle θ , such that θ -θ < δ. Proof. (From [START_REF] Shiu | The Shapes and Sizes of Pythagorean Triangles[END_REF]. We quote his proof because an analogous one will be used for Corollary 1.) Let (a, b, c) ∈ PPT. It is well-known that, assuming that a is even, (a, b, c) can be rewritten (2mn, m 2 -

n 2 , m 2 + n 2 )
where m, n are coprime integers satisfying m > n > 0.

It follows that

tan(θ ) = m 2 -n 2 2mn = 1 2 m n - n m .
Let us write t = tan(θ ) and r = m/n so that we have r 2 -2tr -1 = 0 and hence

r = t + t 2 + 1 = tan(θ ) + sec(θ ).
Note that we do not have r = t -√ t 2 + 1 because r > 0. It is now easy to prove the theorem. Let 0 θ π/4 and u = tan(θ) + sec(θ). One can choose a sequence of rational numbers r 0 , r 1 , r 2 , . . . converging to u, where

r k = m k n k , k = 0, 1, 2, . . .
such that m k and n k are positive and coprime [START_REF] Anglin | Using Pythagorean Triangles to Approximate Angles[END_REF].

Let x = 2m k n k , y = m 2 k -n 2 k , z = m 2 k + n 2
k with corresponding angle θ k . The angles θ 0 , θ 1 , θ 2 , . . . tend to θ, therefore the θ k 's can approximate θ arbitrarily closely, as required.

By symmetry of the Pythagorean triples (a, b, c) and (b, a, c), the density of the aforementioned set can be extended to the interval [0, π/2]. A corollary is the density of the set of hyperbolic angles covered by Pythagorean triples in R + . This is stated by Corollary 1 below.

Corollary 1. ∀ϕ ∈ R + , ∀δ > 0, there exists a primitive Pythagorean triple with an associated hyperbolic angle ϕ , such that ϕ -ϕ < δ.

Proof. By analogy, we replace tan(ϕ ) by sinh(ϕ ), and sec(ϕ ) by cosh(ϕ ) in Shiu's proof of Theorem 1. We have

sinh(ϕ ) = m 2 -n 2 2mn = 1 2 m n - n m .
Let us write s = sinh(ϕ ) and r = m/n so that we have r 2 -2sr -

1 = 0 and hence r = s + s 2 + 1 = sinh(ϕ ) + cosh(ϕ ).
Let ϕ 0 and u = sinh(ϕ) + cosh(ϕ). One can choose a sequence of rational numbers r 0 , r 1 , r 2 , . . . converging to u, where

r k = m k n k , k = 0, 1, 2, . . .
such that m k and n k are positive and coprime. Let

x = 2m k n k , y = m 2 k -n 2 k , z = m 2 k + n 2
k with corresponding hyperbolic angle ϕ k . It follows that the hyperbolic angles ϕ 0 , ϕ 1 , ϕ 2 , . . . tend to ϕ and that the ϕ k 's can approximate ϕ arbitrarily closely, as required.

Although we are now certain that there will always be an infinite number of PPTs for each row of our tables, these theorems do not give any bounds on the size of the PPTs. In practice, we will see that the PPT sizes allow us to build tables in double precision indexed by usual values of indexing bits (up to 12). Also, we still have to find an easy means to generate PPTs efficiently. We choose to use the Barning-Hall tree [START_REF] Freek | Over Pythagorese en bijna-Pythagorese driehoeken en een generatieproces met behulp van unimodulaire matrices [On Pythagorean and quasi-Pythagorean triangles and a generation process with the help of unimodular matrices[END_REF][START_REF] Hall | 232. Genealogy of Pythagorean Triads[END_REF], which exhibits a ternary structure that links any PPT to three different PPTs. It allows to recursively generate PPTs quite efficiently: From any PPT represented as a column vector, the Barning-Hall tree allows to compute three new PPTs by multiplying the former with the three matrices

  1 -2 2 2 -1 2 2 -2 3     -1 2 2 -2 1 2 -2 2 3     1 2 2 2 1 2 2 2 3   .
(2.4)

It has been proven that all PPTs can be generated from the root [START_REF] Abramowitz | Handbook of Mathematical Functions With Formulas[END_REF][START_REF] Anglin | Using Pythagorean Triangles to Approximate Angles[END_REF][START_REF] Tom | Introduction to analytic number theory[END_REF] with increasing hypotenuse lengths [START_REF] Price | The Pythagorean Tree: A New Species[END_REF].

For every generated PPT (a, b, c), we also consider its symmetric PPT (b, a, c), because it may be interesting for three reasons: (b) Whenever arcsin(b/c) ≈ π/4, one also has arcsin(a/c) ≈ π/4, and both triples may either fall into two different rows of the exact lookup table, which can help reduce the value k (since they share a common hypotenuse), or fall into the same subinterval, which can help find a better corrective term.

2.

For the hyperbolic functions, since the reduced argument x lies in [0, ln(2)/2], interesting PPTs must satisfy a/b ∈ sinh([0, (n + Hence the first step of PPT generation using the Barning-Hall tree is the following: multiplying the matrices in Equation (2.4) by the root (3, 4, 5) taken as a column vector, one gets the three new PPTs [START_REF] Tom | Introduction to analytic number theory[END_REF][START_REF] Bernhart | Pythagoras' garden, revisited[END_REF][START_REF] Boehm | On de Casteljau's algorithm[END_REF], [START_REF] Boldo | A Simple Test Qualifying the Accuracy of Horner's Rule for Polynomials[END_REF][START_REF] Beadie | Refractive index measurements of poly(methyl methacrylate) (PMMA) from 0.4 -1.6 µm[END_REF][START_REF] Brisebarre | M, p, k)-Friendly Points: A Table-Based Method for Trigonometric Function Evaluation[END_REF], and [START_REF] Brunie | Modified Fused Multiply and Add for Exact Low Precision Product Accumulation[END_REF][START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF][START_REF] Cody | Implementation and Testing of Function Software[END_REF], and their symmetric counterparts [START_REF] Bernhart | Pythagoras' garden, revisited[END_REF][START_REF] Tom | Introduction to analytic number theory[END_REF][START_REF] Boehm | On de Casteljau's algorithm[END_REF], [START_REF] Beadie | Refractive index measurements of poly(methyl methacrylate) (PMMA) from 0.4 -1.6 µm[END_REF][START_REF] Boldo | A Simple Test Qualifying the Accuracy of Horner's Rule for Polynomials[END_REF][START_REF] Brisebarre | M, p, k)-Friendly Points: A Table-Based Method for Trigonometric Function Evaluation[END_REF], and [START_REF] Brunie | Contributions to computer arithmetic and applications to embedded systems[END_REF][START_REF] Brunie | Modified Fused Multiply and Add for Exact Low Precision Product Accumulation[END_REF][START_REF] Cody | Implementation and Testing of Function Software[END_REF]. In the following, note that we always consider the "degenerated" PPT (0, 1, 1), because it gives us an exact corrective term for the first table entry, without making the LCM k grow.

1/2) • 2 -p ]) ⊃ [0, √ 2/4]. Therefore, if a 4b/ √ 2,

Selection of primitive Pythagorean triples

For each row i of the table indexed by x h • 2 p , we want to select exactly one PPT (a, b, c) with a corresponding value θ = arcsin(a/c) or ϕ = arsinh(a/b) such that:

|x h -θ| < 2 -p-1 for trigonometric functions, or |x h -ϕ| < 2 -p-1 for hyperbolic functions.
The existence of such values θ and ϕ is a consequence of Theorem 1 and its Corollary 1.

In the following, we define our corrective terms, which we denote by corr, as: corr = θ -x h for trigonometric functions, or corr = ϕ -x h for hyperbolic functions.

Once one PPT has been selected for each row i, a naive solution would consist in storing exactly each a i , b i , c i in the lookup table T , plus an approximation of corr on as many bits as necessary. Instead, as presented in Section 2.3, we suggest to store two integers C h and S h of the form

C (i) h = b i c i • k and S (i) h = a i c i • k (2.5)
for the trigonometric functions, or

C (i) h = c i b i • k and S (i) h = a i b i • k (2.6)
for the hyperbolic functions, where k ∈ N * is the same for all table entries. In order to reduce the memory footprint of the table T , we look for small table entries C h and S h . This entails looking for a value k that is rather small, which can be simplified as a search for a small value in a set of least common multiples.

We showed in [START_REF] De | Range reduction based on Pythagorean triples for trigonometric function evaluation[END_REF] that a straightforward approach, consisting in finding the minimum of the set of all possible LCMs, was quickly out of reach with current technology, as soon as tables had more than a few rows. For example, for the trigonometric functions, Figure 2.7 shows the number of PPTs per entry for p = 7, i.e. when there are 102 entries, such that the generated hypotenuses were less than 2 18 . 18 corresponds to the minimum number of bits for the hypotenuses so that there be at least one PPT per entry. In this figure, the mean number of PPTs per entry is about 413, and the standard deviation is close to 33. Because we are looking for a small least common multiple k of one PPT hypotenuse per 413 symbols: this corresponds to 413 101 ≈ 10 173 combinations! This is obviously way too many -there are only about 10 80 atoms in the observable universe [39, pp. 37-43] -to permit an exhaustive search. Indeed, even if we could test a billion combinations per second, it would take approximately 10 156 years to find the right one, which is about 10 146 as much as the age of our universe.

It can be observed in Figure 2.7 that the number of PPTs for the table entry 0, that is, near the angle 0, is equal to 106 whereas the mean is around 413. Actually, there are advantages of not even considering this entry, by selecting the degenerated PPT (0, 1, 1). This PPT corresponds to the angle θ = 0 • or ϕ = 0 exactly, with a corrective term corr = 0. We see three advantages of selecting this PPT:

1. The search space is reduced: we need not worry about the first entry;

2. Its hypotenuse is equal to 1, which will not impact the search of LCM, as lcm(k, 1) = 1.

3.

Its corresponding corrective term is exact and null: corr = 0.

In Figure 2.8, we have reported the number of PPTs per entry for p = 7, such that the generated hypotenuses were less than 2 14 . Here, 14 corresponds to the minimum number of bits for the hypotenuses so that there be at least one PPT per entry, with the first entry excluded from the search space. In this figure, the mean number of PPTs per entry is about 26 and the standard deviation is close to 5. Therefore, if the first entry is excluded, the search space is reduced to about number of combinations to compute the smallest LCM in a reasonable time. To further reduce the space search, the LCM will not be computed but rather be looked for amongst generated hypotenuses.

To reduce the computational time, the solution we propose in [START_REF] De | Range reduction based on Pythagorean triples for trigonometric function evaluation[END_REF] consists in looking for an LCM directly amongst generated values: hypotenuses (for the trigonometric functions), or bigger legs (for the hyperbolic functions) [START_REF] De | Exact Lookup Tables for the Evaluation of Trigonometric and Hyperbolic Functions[END_REF]. We call this solution "exhaustive search" in the sequel, as it always gives the smallest LCM. This claim may be counter-intuitive since one may ask why the smallest LCM would necessarily appear as the denominator side of some PPT. We prove this claim in [START_REF] De | Exact Lookup Tables for the Evaluation of Trigonometric and Hyperbolic Functions[END_REF]. It is done in three steps, which give Theorems 3 and 4. First, we start by proving that for any table, there is a set of PPTs that fills it, for which the denominator sides (either hypotenuses or bigger legs) have the least LCM possible:

Lemma 1. If k ∈ N * is
the least element in the set of possible LCMs for a table of n rows, then there is a set of primitive Pythagorean triples {T i = (a i , b i , c i )} i∈{0,...,n-1} , with T i belonging to row i, such that: lcm(c 0 , . . . , c n-1 ) = k for the trigonometric functions, and lcm(b 0 , . . . , b n-1 ) = k for the hyperbolic functions.

Proof. By construction. For readability, and without loss of generality, we will only consider the LCM for the trigonometric functions. The set of possible LCMs for a table T of n rows is non-empty as a consequence of Theorem 1. Thus, by the Well-Ordering Principle, we know that k exists [5, p. 13].

Since k exists, let T 0 , . . . , T n-1 be a set of Pythagorean triples that fill the table T , such that lcm(c 0 , . . . , c n-1 ) = k. Note that the T i s need not necessarilly be primitive. If there exists i ∈ {0, . . . , n -1} such that T i is non primitive, then there exists an integer α > 1 and

T i ∈ PPT such that T i = α • T i . Thus, we have k = lcm(c 0 , . . . , α • c i , . . . , c n-1 ).
Hence α divides k. Since k is the least element of the set of possible LCMs, k/α < k cannot be a valid LCM. Therefore the prime factors of α must be shared with other c j . Hence, it is possible to write k = lcm(c 0 , . . . , c i , . . . , c n-1 ).

As mentioned in Section 2.4.1, T i and T i share the same angles θ and ϕ, and consequently they belong to the same subdivision of the table T . By repeating this process c i by c i as long as necessary, one can construct a set of primitive Pythagorean triples {T i } i∈{0,...,n-1} such that k = lcm(c 0 , . . . , c n-1 ), which concludes the proof.

Second, let us recall two important results concerning primitive Pythagorean triples.

Theorem 2. Let n be a positive integer. First n is the hypotenuse of a primitive Pythagorean triple if and only if all of its prime factors are of the form 4k + 1. Second, when n > 2, n is a leg of a primitive Pythagorean triple if and only if n ≡ 2 mod 4.

Proof. For the first result, see [START_REF] Sierpi | Elementary Theory of Numbers: Second English Edition[END_REF]ch. XI.3,[START_REF] Vella | When is n a member of a Pythagorean triple?[END_REF]. For the second result, see [START_REF] Albert | Recreations in the theory of numbers: The queen of mathematics entertains[END_REF]p. 116,[START_REF] Sierpi | Elementary Theory of Numbers: Second English Edition[END_REF]ch. II.3]. Now, for the trigonometric functions, we formulate Theorem 3, which states that the smallest LCM k is the hypothenuse of a PPT. Theorem 3. If k ∈ N * is the least element of the set of possible LCMs for a table of n rows for trigonometric functions, then k is the hypotenuse of a primitive Pythagorean triple.

Proof. Lemma 1 taught us that each T i could be made primitive. Thus, by Theorem 2, we have:

∀i ∈ {0, . . . , n -1} , ∃ β p,i p∈P π ∈ N P π , c i = p∈P π p β p,i
where P π is the set of Pythagorean primes, i.e. the set of prime numbers congruent to 1 modulo 4.1 By construction, we have

k = lcm(c 0 , . . . , c n-1 ) = p∈P π p max i β p,i
which means that k is a product of Pythagorean primes. By Theorem 2, k is the hypotenuse of at least one primitive Pythagorean triple, which concludes the proof.

Finally, for the case of hyperbolic functions, Theorem 4 shows that the smallest LCM k is a leg of a PPT. Theorem 4. If k ∈ N * is the least element of the set of possible LCMs for a table of n rows for hyperbolic functions, then k is a leg of a primitive Pythagorean triple.

Proof. By Theorem 2, we only need to prove that k ≡ 2 mod 4.

Using Lemma 1, we know that each T i = (a i , b i , c i ) can be made primitive. Thus, using the aforementioned equivalence, we have ∀i ∈ {0, . . . , n -1} , b i ≡ 2 mod 4.

(2.7) Hence, either every b i is odd, or there exists i ∈ {0, . . . , n -1} such that b i is even. In the first case, the LCM k of all b i 's is obviously odd too. In the latter case, b i is even implies that it is a multiple of 4 as a consequence of Equation (2.7), which entails that k is also a multiple of 4. Therefore, by construction,

k = lcm(b 0 , . . . , b n-1 ) ≡ 2 mod 4,
which concludes the proof.

Theorems 3 and 4 justify why we call our search amongst generated hypotenuses or legs "exhaustive", because by searching exhaustively amongst increasing PPT hypotenuses or legs, the smallest LCM will be found eventually. However, we will see in the next section that such an exhaustive search is exponentially expensive with p. Therefore we will also detail the design and implementation of efficient heuristics to build exact lookup tables.

implementations and numeric results

In this section, we present how we implemented the proposed method to generate exact lookup tables as described in Section 2.4. We have designed three solutions to look for a small common multiple k. The first solution, presented in Section 2.5.1, is based on an exhaustive search and allows us to build tables indexed by up to 7 bits in a reasonable amount of time. The second solution uses a heuristic approach, which reduces memory consumption during the execution and the search time. Presented in Section 2.5.2, it allows us to build tables indexed by 10 bits much faster than the exhaustive search. Finally, the third solution, originally proposed as a perspective in [START_REF] De | Exact Lookup Tables for the Evaluation of Trigonometric and Hyperbolic Functions[END_REF], uses a much more efficient heuristic for the trigonometric functions. It allows us to build tables up to the limit of binary64-representable entries, that is, indexed by 14 bits, in a few seconds. This third solution might also be adaptable to the hyperbolic functions, as explained in Section 2.5.3.

Exhaustive search

This section first describes the design of our exhaustive algorithm, then discusses the different choices that we must make when several PPTs are available inside one table row for the found LCM. Finally, numeric results are presented for various precisions p along with two toy exact lookup tables for the trigonometric and hyperbolic functions.

Algorithm

As seen in Section 2.4.3, we restrain our search to the set of generated hypotenuses c or to the set of generated legs b, for the trigonometric and hyperbolic functions, respectively. In other words, we require that the LCM k be the hypotenuse or a leg of one of the generated PPTs. Moreover, for the hyperbolic functions, one has sinh(ln(2)/2) < 1, so that the search can safely be limited to the set of bigger legs only.

To perform such an exhaustive search, we designed a C++ program that takes as input the number p of bits used to index the table T and tries to generate an exact lookup table. To achieve this, the algorithm looks for the smallest integer k, amongst the generated values, that is a multiple of at least one PPT hypotenuse or leg per entry. By Theorems 3 and 4, this search is guaranteed to find the smallest LCM. The algorithm is the following: Start with a maximum denominator size (in bits) of n = 4 and then follow these four steps:

1. Generate all PPTs (a, b, c) such that c 2 n for the trigonometric functions, or c 2 n / cos ((N -1/2) • 2 -p ), where N is the number of table entries, for the hyperbolic functions. The latter inequality guarantees that every interesting leg for the search step will be generated. It is during this step that the Barning-Hall matrices from Equation (2.4) are used.

2.

Store only the PPTs that belong to a table entry, i.e. the potentially interesting ones.

3.

Search for the LCM k among the PPT hypotenuses c or legs b that lie in [2 n-1 , 2 n ], for the trigonometric or the hyperbolic functions, respectively.

4.

If such a k is found, build values (S h , C h , RN j (corr i )) for every row using Equation (2.5) or (2.6) and return an exact lookup table. Otherwise, set n ← n + 1 and go back to step 1.

Possible options for selecting only one triple per row

Sometimes, whenever k is found, several PPTs are possible for a same table row. To build the final table, however, only one PPT per row is to be used. Then the selection of only one of the available PPTs depends on which goal we want to achieve. Said another way, the choice of one PPT over another may simply rely on an optimization function to minimize or maximize. Here, we make the choice of selecting the PPT for which arcsin(a/c) or arsinh(a/b) is the closest to x h , as this minimizes the ratio between corr and x . This virtually increases the precision of the corrected reduced argument x -corr.

But we see at least three alternatives to this solution:

1. We could try to minimize the value max i (corr i ) -min j (corr j ). This could be interesting to store smaller corrective terms δ i = corr i -min j (corr j ), and incorporate min j (corr j ) into the polynomial approximations. But minimizing such a value might be much more expensive as its time complexity is n-1 i=0 m i , where m i is the number of possible triples for the entry i.

2.

In a way similar to Gal's strategy explained in Section 2.2.2, we could choose the PPT for which the corrective term has as many identical bits as possible after its least significant stored bits. If there are enough, the precision of the corrective term is also virtually extended, which could allow for smaller tables.

3. Finally, we could relax the constraints on the k i 's to allow greater values (i.e. non optimal) of the LCM k, so that the corrective terms may have even lower magnitudes or be even closer to machine numbers.

The last option can be done in different fashions, but we will only describe one to explain the idea to the reader. As we have seen in Section 2.4.3, the greater p, the lower the magnitudes of the corrective terms. Therefore it becomes possible to build tables for lower values of p with exponentially lower magnitudes of corrective terms. For instance, imagine that we want to build a hyperbolic table indexed by 7 bits. Instead of building the one for which k is optimal, let us first build a table indexed by 10 bits. Now let us call i 7 and i 10 the row indices of both tables, respectively. Then, we only need to keep the rows i 10 from the table indexed by 10 bits such that:

i 10 • 2 -10 = i 7 • 2 -7 ,
which means the set

i 10 = i 7 • 2 3 | i 7 ∈ 0 . . 2 6 ln 2 = {i 10 = 8i 7 | i 7 ∈ [0 . . 44]} .
Since the corrective terms in the table indexed by 10 bits have a magnitude bounded by 2 -11 , then so will the corrective terms of our improved table indexed by 7 bits, whereas they were only bounded by 2 -8 for a classic table. Therefore, just by doing this, it is possible to reduce the bound on the corrective terms by a factor 8 while keeping exactly stored values of sinh and cosh. The reader may object that for some configurations of tables, say indexed by p bits and p bits such that p < p , the last row cannot be built because max i p < 2 p -p max i p . This case happens e.g. for p = 10 and p = 13 for the hyperbolic functions: max i 13 = 2 12 ln 2 = 2839 and max i 10 = 2 9 ln 2 = 355, which gives no solution for max i 13 = 2 3 • 355 = 2840. In such a case, selecting the last row of the table indexed by p bits would suffice to build the last row of the table indexed by p bits. This, however, may theoretically increase the upper bound on the magnitudes of the improved corrective terms up to 2 -p-1 because of the last row. Instead, it is possible to build a custom table indexed by p bits with max i p = 2 p -p max i p . In our example, this would correspond to a table with 2841 rows. This technique could increase the precision of each corrective term at the cost of greater PPTs. As long as the exact stored values are still machine-representable, then this method remains interesting for software implementations. However, for hardware implementations, this might be less interesting since it involves more storage bits than a regular table.

Numeric results

Tables 2.1 and 2.2 show the results obtained for the trigonometric and hyperbolic functions, respectively. Timings were measured on a server with an Intel Xeon E5-2650 v2 @ 2.6 GHz processor (16 physical cores) with 125 GB of RAM running on GNU/Linux. For the number p of bits that is targeted, the tables describe the value k that was found, followed by the number n of bits used to represent k (that is, n = log 2 (k) ), the time (in seconds) taken by our program to run, and the numbers of PPTs and denominators considered during the LCM search. As can be seen, it was possible to find k and to build tables indexed by up to p = 7 bits in a reasonable amount of time. However, it is clear that the number of dynamic memory allocations, which are mainly used to store the triples and the denominators, grows exponentially with p. Consequently, it was not possible to find k for p 8 with our hardware configuration. Table 2.3 describes an exact lookup table for the trigonometric functions when p = 4, where k = 5525 and the absolute value of the corrective terms is at most 0x1.d6513b89c7237p-6, that is, ≈ 0.0287 for input index i = 5. Table 2.4 presents an exact lookup table for the hyperbolic functions when p = 5, where k = 10 080 and the absolute value of the corrective terms is at most 0x1.7c74108520aebp-7, that is, ≈ 0.0116 for input index i = 3. Figures 2.9 stored triples by dots and the limits of each entry interval by dashed lines.

Characterizing exhaustive results: a first heuristic

To build tables indexed by a larger number of bits, a more efficient solution must be found. In order to drastically reduce the search space, we have developed two heuristics, one for the trigonometric functions and one for the hyperbolic functions. These heuristics try to characterize the denominator sides (hypotenuses or bigger legs) to either keep or reject PPTs during the generation step. k Prime factorization

144 2 4 • 3 2 840 2 3 • 3 • 5 • 7 10 080 2 5 • 3 2 • 5 • 7 180 180 2 2 • 3 2 • 5 • 7 • 11 • 13 1 081 080 2 3 • 3 3 • 5 • 7 • 11 • 13
For this purpose, let us observe the decomposition in prime factors of each k found using the exhaustive search. Such decompositions are given in Tables 2.5 and 2.6 for both the trigonometric and the hyperbolic tables. These factorizations show that every k in the table is a composite number divisible by relatively small primes. Furthermore, for the trigonometric functions, all those small primes are of the form 4n + 1, as per Theorem 2.

Therefore, each heuristic that we propose follows a simple rule: for the trigonometric functions, only store PPTs with a hypotenuse of the form: For the hyperbolic functions, the heuristic only stores PPTs with a bigger leg of the form:

c = p∈P π p r p with        r p ∈ N * if p = 5 r p ∈ {0, 1} if 13 p 73 r p = 0 if p > 73 (2.8)
b = p∈P p r p with        r p ∈ N * if p < 5 r p ∈ {0, 1} if 5 p 23 r p = 0 if p > 23 (2.9)
where P is the set of all primes:

P = {2, 3 , 5, 7, 11, 13, 17, 19, 23, . . . } .
Results, timings, and numbers of considered triples and potential LCMs (hypotenuses or legs) for this heuristic are given in Tables 2.7 and 2.8 for the trigonometric and the hyperbolic functions, respectively. As can be seen, this algorithm considers a number of potential LCMs several orders of magnitude less than the exhaustive search solution. This reduces drastically the memory usage and execution times. For instance, for the trigonometric functions, when p = 7, only 3453 triples are stored, compared to the 1 347 954 triples for the exhaustive algorithm. In this first case, the execution time was reduced from 31 seconds to 0.4 seconds. And for the hyperbolic functions, when p = 7, only 1743 triples are stored, while there were 147 748 when using the exhaustive algorithm. In this second case, the execution time was reduced from 328 seconds to 0.3 seconds.

With this heuristic, the bottleneck is no longer the memory but the selection of PPTs during their generation. Indeed, this selection is carried out on the elements of an exponentially-growing set. And checking if a given integer satisfies either Equation (2.8) or Equation (2.9) requires checking if it is multiple of certain prime numbers, which is a rather expensive test involving integer division.

A much faster heuristic using Fässler's algorithms

In [START_REF] De | Exact Lookup Tables for the Evaluation of Trigonometric and Hyperbolic Functions[END_REF], we concluded that it could be interesting to further characterize potential small LCMs in order to speed up the table precomputation process even more. On one hand, since precomputing large tables (p 10) was still exponentially expensive, we were unable to generate them for p > 10 in a reasonable time. On the other hand, such large tables are rarely useful because of the exponential memory footprint they involve compared to the small computing benefits they bring in return. However, in this section, we will detail another heuristic that allows us to go a little bit further, for two reasons. First, we find this theoretical problem mathematically interesting in itself. Second, and more importantly, as we have seen in Section 2.5.1, solving it could help us improve the tables for lower values of p. Furthermore, the theoretical background of this heuristic gives us lower bounds for the optimal k.

How to generate Pythagorean triples sharing a common hypotenuse?

Before we present our second heuristic, let us observe the variations of the previously found LCMs bit-lengths n with respect to the table index sizes p. They are represented in Figure 2.11 and Figure 2.12, for the trigonometric and hyperbolic functions, respectively. Using the gnuplot 4.6 stats command1 , we performed ordinary least squares (OLS) linear regressions on both datasets, showing very good linear correlation coefficients between n and p. We can see that if the extrapolation is meaningful, which seems a plausible hypothesis, then it would be impossible to store tabulated values as exact binary64 floating-point numbers starting from p ≈ 13 or 14 for the trigonometric functions and from p ≈ 15 or 16 for the hyperbolic functions.

Therefore, we choose to set our ultimate objective to generating trigonometric and hyperbolic tables in a reasonable amount of time, and indexed by at most those values (p = 14 or p = 16).

In the following, we focus only on the trigonometric functions. We propose to directly compute Pythagorean triples (PTs) sharing a common hypotenuse instead of generating huge sets of PPTs and then selecting the relevant ones. To do this, we will use some of Fässler's inspiring mathematical results taken from [64,

a = 2xy, b = x 2 -y 2 , c = x 2 + y 2 with x > y > 0, gcd(x, y) = 1, x + y ≡ 1 mod 2.
(2.10)

Second, as already stated in Theorem 2, the hypotenuse c is necessarily of the form

c = p β 1 1 • p β 2 2 • • • p β n n (2.11)
with pairwise different primes p i ≡ 1 mod 4. Fässler assumes that each β i 1, and we follow his assumption. Then, he reminds us about a famous theorem: Theorem 5 (Fermat's theorem on sums of two squares). A prime number p is of the form 4m + 1 if and only if it can be written as a unique sum of two squares p = x 2 + y 2 , up to a permutation of the squares:

∃(x, y) ∈ N 2 | p = x 2 + y 2 ⇐⇒ p ≡ 1 mod 4.
With these results in mind, Fässler uses the following identity, known as either the Brahmagupta-Fibonacci identity or the Diophantus identity [START_REF] Stevin | L'arithmétique de Simon Stevin de Bruges, revue[END_REF]§III.19]

(r 2 + s 2 ) • (u 2 + v 2 ) = (ru ± sv) 2 + (rv ∓ su) 2
(2.12)

to show that since each prime factor p i of c can be represented as a sum of two squares, then it is recursively possible to write c as 2 n-1 sums of two squares where n is the number of prime factors of c. Furthermore, the representations given by Equation (2.12) comply with all the requirements of Equation (2.10), so that all the 2 n-1 sums of two squares can be used to generate 2 n-1 PPTs sharing the same hypotenuse c using Equation (2.10). For instance, there are only two PPTs that share the common hypotenuse 65 = 5 • 13: First, we write 5 = 1 2 + 2 2 and 13 = 2 2 + 3 2 . Using Equation (2.12), we have:

65 = 1 2 + 2 2 • 2 2 + 3 2 = (2 ± 6) 2 + (3 ∓ 4) 2 = 8 2 + 1 2 = 4 2 + 7 2 .
Then using Equation (2.10), we can compute the only two PPTs (a, b, c) that share the hypotenuse c = 65, namely [START_REF] Brisebarre | Efficient polynomial L ∞ -approximations[END_REF][START_REF]Evidence for the spin-0 nature of the Higgs boson using AT-LAS data[END_REF][START_REF] Fog | Instruction tables: Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD and VIA CPUs[END_REF] and [START_REF] Debjit | Faithful Bipartite ROM Reciprocal Tables[END_REF][START_REF] Damouche | Improving the numerical accuracy of programs by automatic transformation[END_REF][START_REF] Fog | Instruction tables: Lists of instruction latencies, throughputs and micro-operation breakdowns for Intel, AMD and VIA CPUs[END_REF]. Actually, Fässler acknowledges that these results were progressively developed by Fermat, Euler, Lagrange and Legendre. A corollary of these results is what Fässler names the "Formula of Legendre":

Corollary 2 (Formula of Legendre). If c has the form in Equation (2.11), then there are exactly M Pythagorean triples, where

M = (2β 1 + 1) • (2β 2 + 1) • • • (2β n + 1) -1 2 .
Since we know the number of Pythagorean triples that are needed for a given table, we are now able to find a more precise lower bound of k thanks to Corollary 2.

Improved lower bound for the optimal least common multiple Corollary 2 gives us the exact number of Pythagorean triples (PTs) sharing the same hypotenuse c. Hence the lower bound on k can be improved by requiring that the number M of PTs sharing the same hypotenuse k be at least 2 p-1 π/2 . For instance, for p = 5, our trigonometric table has 26 entries, so we need to find at least 25 PTs sharing the same hypotenuse c. There is one less triple needed because we exclude the first entry from the search, as seen in Section 2.4.3.

Let us assume for a moment that k = p 1 • p 2 • • • p n , i.e. that each

β i = 1 for 1 i n.
From Corollary 2, we have (3 n -1) /2 triples sharing the same hypotenuse c. In order to have at least 25 PTs with the least k, we need the first n Pythagorean primes such that:

3 n -1 2 25 ⇐⇒ n log (51) log(3) ⇐⇒ n 4.
So if each β i = 1 and n = 4, there are exactly (3 4 -1)/2 = 40 PTs. Thus the minimum value of such a hypotenuse is the product of the first four Pythagorean primes: c 5 • 13 • 17 • 29 = 32 045. Actually, our first assumption that each β i = 1 may be too simple. Indeed, for n 2, if we remove the greatest Pythagorean prime factor of c, and add x to β 1 = 1 so that now β 1 = 1 + x, we have:

M = (2x + 3) • 3 n-2 -1 2 .
Hence, solving M 3 n -1 2 for x, one finds x 3. This means that for any n ∈ N * , if β 1 4, then there are at least as many PTs sharing the hypotenuse c 1 = 5 β 1 2 i n-1 p i as there are PTs sharing the hypotenuse c 2 = 1 i n p i . Solving c 1 < c 2 for p n , we find p n > 5 β 1 with β 1 4, i.e. p n 625. The least p n satisfying this inequality is p 55 = 641. Therefore,

c 1 = 5 4 • 2 i 54 p i < c 2 = 1 i 55 p i
is the least example showing than c 2 may not be a lower bound for k. Fortunately, we have no practical need for such a consideration because M = 3 54 -1 2 is greater than 10 25 , which is way greater than our practical use: M generally does not exceed 10 5 for our tables. Hence we will compute a lower bound on k with each β i set to 1.

β i = 0 ∀i log(2N + 1) log (3)
and β i = 1 otherwise.

Then our heuristic follows these six steps:

1. Generate an arbitrarily small set K of possible LCMs k, using 1 β 1 3 for p 1 = 5 and β i ∈ {0, 1} otherwise, according to the lower bound found for k.

2. Set c ← min {c ∈ K}.

3. Generate all PTs sharing the same hypotenuse c = i p β i i , by composing Equation (2.12); 4. If there is at least one PT in each table entry, generate the final table and terminate;

5. Else, remove c from K.

6.

If K = ∅, return an error. Else go back to Step 2.

The results using this heuristic are summarized in Table 2.9. As opposed to the linear model prevision that we drew from Figure 2.11, you can see that our heuristic was not able to find an LCM k that could be stored exactly in a single binary64 floating-point number for p = 13. However, we can see that this method allows to build exact trigonometric tables up to p = 12 within reasonable amounts of time and memory usage, and that it finds the same LCMs k as both the exhaustive method and the first heuristic up to p = 10, which comforts the results found for p ∈ {11, 12}.

In this section, we did not address the generation of exact lookup tables for the hyperbolic functions. For this problem, an efficient algorithm that generates all the Pythagorean triples sharing a common bigger leg would be needed. A first lead might be to follow the results presented in [START_REF] Rothbart | Pythagorean Triples: A New, Easy-to-Derive Formula With Some Geometric Applications[END_REF] by Rothbart and Paulsell. An open question remains to be answered: does this heuristic actually give optimal results for binary64 exact lookup tables whenever we use unlimited values for the β i 's?

comparisons with other methods

We have presented a range reduction for the trigonometric and hyperbolic functions based on exact lookup tables and a method to efficiently build such tables. In order to compare this solution with the other solutions presented in Section 2.2.2, we consider a twophase evaluation scheme for the trigonometric and hyperbolic functions that targets correct rounding for the rounding to nearest in double precision. The quick and the accurate phases target a relative error less than 2 -66 and 2 -150 , respectively [START_REF] Lefèvre | Hardest-to-Round Cases[END_REF][START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. We choose p = 10 which corresponds to 805 rows in the trigonometric table and 356 rows in the hyperbolic table.

In order to ease comparisons, we consider only the number of memory accesses required by the second range reduction and the number of floating-point operations (FLOPs) involved in the reconstruction step, and table sizes. We will consider that expansion algorithms are used whenever high accuracy is required as it is the case in the correctly rounded library CR-Libm [START_REF] Daramy-Loirat | A library of correctly rounded elementary functions in double-precision[END_REF]. Let us recall that an expansion of size n consists in an unevaluated sum of n nonoverlapping floating-point numbers that represents a given number with a larger precision than the one available in hardware [START_REF] Richard | Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates[END_REF]. We will take Table 2.10 extracted from [34, § 2.3, 101, § 3.2, 102] as the reference costs for those algorithms. The notation E n stands for an expansion of size n in double precision, so that, with this notation, E 1 represents a regular binary64 word.

E 2 = E 1 + E 2 11 E 2 = E 2 + E 2 12 E 3 = E 1 + E 3 16 E 3 = E 3 + E 3 27
E 2 = E 1 × E 2 20 6 E 2 = E 2 × E 2 26 11 E 3 = E 1 × E 3 47 19 E 3 = E 3 × E 3 107 43
2.6.1 Costs for the trigonometric functions Tang's solution

In order to reach an accuracy of 66 bits, Tang's solution requires access to tabulated values S h and C h that are stored as expansions of size 2. These values need to be multiplied by the two results of the polynomial evaluations, which can be considered stored as expansion of size 2 as well. The total cost of the quick phase becomes: 4 doubleprecision memory accesses, 2 multiplications E 2 × E 2 , and 1 addition

E 2 + E 2 , that is, 34 FLOPs (64 without FMA).
In case the quick phase failed to return correct rounding, the accurate phase is launched. This phase requires access to 2 extra tabulated values to represent S h and C h as expansions of size 3. Those values need to be multiplied by the two results of the polynomial evaluations, which can be considered stored as expansion of size 3 as well. The total cost of the accurate phase becomes: 2 extra memory accesses, 2 multiplications E 3 × E 3 , and 1 addition E 3 + E 3 . This corresponds to 6 memory accesses, 113 FLOPs (241 without FMA), and a 38 640 byte table.

Gal's solution

Using Gal's method, the corrective terms allow around 63 bits of accuracy, and Stehlé and Zimmermann's improvement allows to reach 74 bits. By considering Stehlé's approach, only one double-precision number is required for S h , C h and the corrective term, which fits on about 20 bits, to reach an accuracy of 66 bits. Again, S h and C h need to be multiplied by the two results of the polynomial evaluations, which can be considered stored as expansions of size 2. Thus the quick phase requires 2 + 1 double-precision memory accesses, 1 addition E 2 = E 1 + E 2 for the corrective term, 2 multiplications E 1 × E 2 , and 1 addition E 2 + E 2 . The cost of the quick phase with this table becomes 3 memory accesses and 35 FLOPs (63 without FMA).

To reach the 150-bit accuracy required by the accurate phase, it is necessary to get 2 extra floating-point numbers for S h and C h , which are not exact. The corrective term is then incorporated in the final result using an addition E 3 = E 1 + E 3 . The remaining operations need to be done using size-3 expansions. The total cost for the accurate phase becomes: 4 extra memory accesses, 2 multiplications E 3 × E 3 ,

1 addition E 1 + E 3 , and 1 addition E 3 + E 3 , that is, 7 memory accesses, 129 FLOPs (257 without FMA), and a 45 080-byte table.

Exact lookup tables solution

With our solution, as shown in Table 2.7, at most 39 bits are required to store S h and C h , that is, only one floating-point number per entry. Our corrective terms are inexact, so that the cost of the quick phase is the same as Gal's approach in Section 2.6.1, except that the addition for the corrective term is an addition E 2 = E 2 + E 2 , which accounts for 1 extra FLOP and 1 extra memory access. Hence, the quick phase for the exact lookup table solution is 36 FLOPs (64 without FMA), which is as much as Tang's solution when no FMA is available.

However, for the accurate phase, values S h and C h that were accessed during the quick phase are exact, and do not require any extra memory access. The corrective term is stored as an expansion of size 3 and it requires 1 extra memory access to reach the 150 bits of accuracy. The addition for the corrective term is performed using an addition with expansions of size 3. Multiplications correspond to E 3 = E 1 × E 3 as the results of the polynomial evaluations can be considered as expansions of size 3. The final addition is done using an E 3 = E 3 + E 3 operation. That is, the total cost of this step becomes 5 memory accesses and 92 FLOPs (148 without FMA), for a 32 200 byte table.

Comparison results

Table 2.11 synthesizes the estimated costs for those three range reduction algorithms based on tabulated values. This table reports the number of memory accesses and FLOPs for the quick and accurate phases, together with the size in bytes of each row of the precomputed table. First, it can be seen that the proposed table-based range reduction requires less memory per table row than other solutions. Tang's method requires 48 bytes per row, Gal's method 56 bytes and our exact lookup table 40 bytes. This is an improvement in memory usage of ≈ 17% and ≈ 29% over Tang's and Gal's methods, respectively. Second, regarding the number of operations, our solution requires two more FLOPs than Tang's solution for the quick phase (0 without FMA), and one more FLOP than Gal's solution. This represents an overhead of ≈ 6% (0% without FMA) and ≈ 2%, respectively. This small downside comes with at least three advantages:

1. The accurate phase shows greater gains. Indeed, for this phase, there is an improvement in favor of our approach of ≈ 19% (39% without FMA) and ≈ 29% (42% without FMA) over Tang's and Gal's methods, respectively.

2. Our method reduces the number of memory accesses during the accurate phase. Indeed, the number of accesses is reduced from 7 to 5 compared to Gal's approach. This is an improvement of ≈ 29%, which translates to ≈ 17% compared to Tang's solution.

3. Error bounds might be easier to determine, as half the terms in the reconstruction are exact. Following this, the maximum acceptable error for polynomials might be loosened and lead to polynomial approximations of lesser degree, which would improve the overall performance.

Costs for the hyperbolic functions

The costs for memory accesses and table sizes for the hyperbolic functions are the same as the ones for the trigonometric functions. Indeed, although the exact lookup tables are smaller for the hyperbolic functions because ln(2)/2 ≈ 0.347 is less than π/4 ≈ 0.785, they both contain as many words per row. We recall the reconstruction step for the hyperbolic functions:

sinh(y) = 2 n-1 -2 -n-1 • (C h • P C (x ) + S h • P S (x )) ± 2 n-1 + 2 -n-1 • (S h • P C (x ) + C h • P S (x )) . (2.13)
If we do not count the computation of the terms 2 n-1 ± 2 -n-1 , there are 6 multiplications, and 3 additions. Indeed these terms can be neglected in the cost estimations as they are easy to generate using shifts and adds in integer arithmetic.

Tang's and Gal's solutions

On the first hand, let us consider Tang's solution: Overall, the table requires 17 088 bytes. For the quick phase, we must use E 2 expansions for 6 multiplications and 3 additions, as we have seen in Equation (2.13). These extended-precisions operations account for 102 FLOPs (192 without FMA). During the accurate phase, all computations are done using E 3 expansions, which require a total of 339 FLOPs (723 without FMA). On the other hand, using Gal's method, the table takes 19 936 bytes and one addition E 2 = E 1 + E 2 for x -corr is added to the cost computations. The quick phase accounts now for 4 multiplications E 1 × E 2 between tabulated values and polynomial results, 2 multiplications E 2 × E 2 for the remaining terms, 1 addition E 2 = E 1 + E 2 , and

3 additions E 2 = E 2 + E 2 .
Hence, the cost of the quick phase with this table is 93 FLOPs (179 without FMA). For the accurate phase, all operations need to be done using size-3 expansions, except for the addition of the corrective term E 3 = E 1 + E 3 . The total cost for the accurate phase becomes: 355 floating-point operations (739 without FMA).

Our exact lookup tables solution

Again, as for the trigonometric functions, the cost of the quick phase is the same as Gal's approach plus one extra memory access and one extra FLOP for the addition of the inexact corrective term, which is 94 FLOPs (180 without FMA). However, the table is smaller as it takes only 14 240 bytes.

For the accurate phase, as it was the case for the trigonometric functions, values S h and C h are exact and do not require any extra memory access. Only the corrective term is stored as an expansion of size 3 and requires 1 extra memory accesses in order to reach an accuracy of 150 bits. Hence, the accurate phase corresponds to 4 multiplications E 1 × E 3 , 2 multiplications E 3 × E 3 , and 4 additions E 3 + E 3 . Thus the total cost of this step is 270 floating-point operations (510 without FMA).

Comparison results

As in Section 2.6.1, we compare the three methods using synthesized results from Table 2.11. Regarding memory usage, we have quantitatively the same benefits as for the trigonometric functions. Regarding the number of floating-point operations, our solution has the same small overhead of one FLOP for the quick phase compared to Gal's solution, but it requires 8% (6% without FMA) less FLOPs than Tang's solution. For the accurate phase, there is an improvement in favor of our approach of ≈ 20% (29% without FMA) and ≈ 24% (31% without FMA) over Tang and Gal, respectively. This is mainly due to the huge cost of multiplications E 3 × E 3 that our method avoids in 4 out of 6 multiplications.

example on the trigonometric sine

We propose to illustrate the use of the proposed error-free tabulated values compared to Tang's solution, which is probably the most pervasive, using Table 2.3 and Table 2.12. We consider the evaluation of the sine function for the input value y = 10 when targeting 8 bits of accuracy. When necessary, storage precision and intermediate computations will be done on 16 bits of precision in order to emulate 8-bit expansions of size 2, like what is commonly used for the quick phase. For clarity, bit fields will be represented in radix 2 (suffix " 2 "). First, the reduced argument x is computed on 16 bits: As explained earlier, x is essentially inexact and has to be rounded (it is rounded to 16 bits in this case). By construction, x belongs to the interval [-π/4, π/4], but in practice, properties of symmetry are used so that it finally lies in [0, π/4]. We have q = 6, which entails q ≡ 2 mod 4. This means that an additional rotation by π radians was made during the first range reduction, so that the sign changed and we now have sin(10) = -sin(x). Now let us compute x h and x , such that x = x h + x with x h = i • 2 -4 , i ∈ [0, 13] and |x | 2 -5 . This leads to

q = 10 • 2/π = 6 x = 10 -q • π/2 = 1.00100110100001 2 • 2 -1 • (1 + ε -16 ) ≈ 0.575225830078125.
x h = 9 • 2 -4 and x = x -x h = 1.10100001 2 • 2 -7 • (1 + ε -16 ).
For row index i = 9, Table 2.12 gives us

S h = 1.000100010000111 2 • 2 -1 • (1 + ε -16 ) C h = 1.101100010001111 2 • 2 -1 • (1 + ε -16 ), while Table 2.3 gives us S h = 2880 C h = 4715 corr = -1.110100000010101 2 • 2 -7 • (1 + ε -16 )
with k = 5525 as mentioned in Section 2.5.1. Now is when the different table-based methods diverge. They will be analyzed separately.

Tang's table lookup method

Using Tang's table, the evaluation of two polynomial approximations to sin(x) and cos(x) is performed with an output precision of 12 bits at input x = 1.10100001 2 • 2 -7 :

sin(x ) ≈ 1.10100001 2 • 2 -7 • (1 + ε -12 ) cos(x ) ≈ 1.0 2 • (1 + ε -12 ).
(2.14)

The final result is then computed on 8 bits using tabulated values S h and C h with the results from Equation (2.14): sin(10) = -sin(x)

≈ -S h • 1.0 2 -C h • 1.10100001 2 • 2 -7 ≈ -1.0001011 2 • 2 -1 • (1 + ε -8 ) ≈ -0.54296875.

Exact table lookup method

With the proposed solution, the evaluation of polynomial approximations to sin(x)/k and cos(x)/k are performed with an output precision of 12 bits at input x -corr ≈ 1.10111000100101 2 • 2 (2.15)

Finally, the result is reconstructed as follows, using the error-free tabulated values S h and C h from Table 2.3 and the approximations from Equation (2.15):

sin(10) = -sin(x) = -S h • cos(x -corr)/k -C h • sin(x -corr)/k ≈ -1.0001011 2 • 2 -1 ≈ -0.54296875.

Comparison between both methods

One can see that for the exact lookup table method, the error is quickly concentrated into the reduced argument x , then into xcorr and finally into the polynomial approximations. The main advantage of this method relies on the fact that S h and C h embed no rounding error at all, which allows for an easier, faster, and more elegant reconstruction step.

conclusions and perspectives

In this chapter, we have presented a new approach to address tablebased range reductions in the evaluation process of elementary functions. We have shown that this method allows for simpler and faster evaluation of these functions. For the trigonometric and hyperbolic functions, we used Pythagorean triples to build exact lookup tables that concentrate the error into the polynomial evaluations. Compared to other solutions, our exact lookup tables eliminate the rounding errors on certain tabulated values, and transfer these errors into the remaining reduced argument, which generally already contains rounding error. When targeting correct rounding in double precision, we have shown that our method allows to reduce the table sizes, the number of memory accesses and the number of floating-point operations involved in the reconstruction step by up to 29%, 29% and 42%, respectively, compared to Tang's and Gal's tables. This benefit could be even higher for greater targeted accuracies since all tabulated values except the corrective terms are exactly stored and do not depend on the targeted accuracy. Therefore we encourage integration and testing of our exact lookup tables into the efficient implementations of the trigonometric functions in the medium-precision range detailed in [START_REF] Johansson | Efficient Implementation of Elementary Functions in the Medium-Precision Range[END_REF].

Finally, in this chapter we have only focused on the trigonometric and hyperbolic functions, as they both benefit from the use of Pythagorean triples. However, the concept remains valid and worthwhile for other functions, provided at least two terms need to be tabulated per row and that exact values can be found for such terms. Hence it may be interesting to investigate whether and how other elementary functions could benefit from these exact lookup tables.

"Multitasking: A clever method of simultaneously slowing down the multitude of computer programs that insist on running too fast." Software Terms, https://www.gnu.org/fun/jokes/software.terms.html

As we have seen, the evaluation of elementary functions often makes use of polynomial approximations. Their theoretical performance is closely linked to the order in which multiplications and additions are performed. Such an order is commonly referred to as an evaluation scheme. In this chapter, we intend to show that the widespread Horner scheme is in fact rarely the most efficient. Other evaluation schemes can better take advantage of modern computer architectures by reducing data dependencies. Vectorization is a key factor to improve throughput performance for large, independent computations [START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF]. Hence, focusing on single-instruction multiple-data vector architectures (SIMD), we benchmarked three classic schemes as well as thousands of automatically generated ones, which we vectorized using auto-vectorizing capabilities of the GNU Compiler Collection (GCC). From our various measurements, we confirmed the important impact on performance of the chosen scheme as well as the input interval for the evaluation. These conclusions are used in Chapter 4, which targets the efficiency of automatically-generated source code.

introduction

Software implementations of elementary functions often use one or more polynomial evaluations of reasonable degree, on a small input interval [123, § 11.4]. The evaluation of those polynomials, i.e. the order in which the arithmetical operations are executed, is usually statically decided at the time the source code is written. This order is called a polynomial evaluation scheme. They may differ for custom implementations for specific architectures, but generic codes often use a single scheme for several architectures. 1Horner's scheme, or sometimes Horner's rule, is one of the most commonly used schemes [90, §4.6.4]. It is often chosen for its simplicity of implementation, its good accuracy properties, and because it minimizes operation count, hence having good potential throughput [START_REF] Boldo | A Simple Test Qualifying the Accuracy of Horner's Rule for Polynomials[END_REF]. Figure 3.1 gives a visual representation of this scheme, which shows that it is essentially sequential. Therefore its latency on architectures benefiting from unbounded parallelism is theoretically higher than the latency of more parallel schemes such as Estrin's [START_REF] Munro | Optimal algorithms for parallel polynomial evaluation[END_REF]. Regardless of potential precision loss, one may want to use a more parallel scheme to reduce the latency or increase the throughput. However, manually selecting a fast scheme on a given hardware can be a very long and complex task as the polynomial degree increases.

Our objective is to automatically generate polynomial schemes that benefit most from features of modern architectures. More particularly, we are interested in leveraging Instruction Level Parallelism (ILP) and data parallelism through the Single Instruction Multiple Data (SIMD) execution model. In this chapter, the throughput performances of different polynomial schemes are studied on an Intel Haswell architecture1 . We present two useful results:

1. High-ILP schemes can have both better latency and throughput than more sequential schemes like Horner, which goes against the usual recommendations [123, § 11.5].

2. We observe that the input interval for a given polynomial can be a key factor for the choice of a fast evaluation scheme, especially when this interval is around zero.

The chapter is organized as follows: Section 3.2 presents the classic polynomial schemes that are studied later on. Section 3.3 analyzes the measured throughputs of several schemes on an Intel Haswell CPU. Finally, Section 3.4 gives perspectives about efficient code generation for polynomial schemes with respect to vectorized implementations and ILP-leveraging environments.

reminders about polynomial evaluation

Let P ∈ F n [X] be a degree-n polynomial with nonzero representable floating-point coefficients. Suppose that we want to evaluate this polynomial at a certain point x in the same floating-point format, using a Simplified, Not Pipelined and Infinitely Parallel Architecture (SNPIP) whose addition and multiplication latencies are respectively 3 and 5 cycles, 2 , and that has no FMA. We do not focus on evaluation schemes with coefficient adaptation. The main reason of such a restriction is that coefficient adaptation can change polynomial coefficients and can lead to catastrophic cancellations near the roots of the polynomial. The interested reader can find out more about coefficient adaptation in [START_REF] Eve | The evaluation of polynomials[END_REF][START_REF] Paterson | On the Number of Nonscalar Multiplications Necessary to Evaluate Polynomials[END_REF]. In our case, the number of additions remains constant and increasing instruction parallelism can only increase the number of multiplications [START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors -Polynomial evaluation-based algorithms and certified code generation[END_REF]. One of the most common schemes is probably Horner's, which is often taught as a general and elegant method to evaluate a polynomial [90, §4.6.4]. It allows to evaluate P = n i=0 a i X i in n additions and n multiplications. For instance, for n = 7, one has:

P(x) = a 0 + x • (a 1 + x • (a 2 + x • (a 3 + x • (a 4 + x • (a 5 + x • (a 6 + x • a 7 )))))). (3.1)
Its interest is threefold:

1. It minimizes the number of multiplications that are needed;

2. There is no need to store intermediate results, which accounts for a lower pressure put on hardware registers and thus potential reductions of register spilling [90, p. 486];

3. In the context of elementary function evaluation, it is known to be numerically stable [14, ch.9].

Therefore, when the selection criterion is accuracy, Horner can be seen as a "good" scheme. But, as seen before, its latency on an architecture like SNPIP is higher than many schemes due to its intrinsic sequentiality. For example, the second-order Horner's scheme, which we denote Horner-2, is parenthesized this way for n = 7:

P(x) = a 0 + x 2 • a 2 + x 2 • a 4 + x 2 • a 6 + x • a 1 + x 2 • a 3 + x 2 • a 5 + x 2 • a 7 . ( 3.2)
This scheme shows some ILP, as illustrated in Figure 3.2: Odd-and even-degree monomials can be evaluated in parallel, at the cost of an additional multiplication. Its latency is thus asymptotically divided by two compared to Horner's rule [90, p. 488].

Estrin's scheme, illustrated in Figure 3.3, shows even more ILP by splitting the evaluation tree according to powers of x of the form 2 i [START_REF] Estrin | Organization of Computer Systems: The Fixed Plus Variable Structure Computer[END_REF]. For example, still for n = 7, Estrin's scheme is: In the case of a degree n = 2 k -1 polynomial, Estrin's evaluation tree is well-balanced, but a bit less in the general case. Its latency on unbounded parallelism is in O(log n).

P(x) = (a 0 + a 1 • x) + x 2 • (a 2 + a 3 • x) + x 4 (a 4 + a 5 • x) + x 2 (a 6 + a 7 • x) . ( 3
Using the SNPIP architecture, for n = 7 and with Equations (3.1), (3.2) and (3.3), we can theoretically expect the latencies given in Table 3.1 for all the aforementioned evaluation schemes. We can see that the choice of an evaluation scheme has a direct impact on its efficiency. Therefore it can be important to choose the right evaluation scheme when efficiency is required. However, the combinatorics of the number of existing schemes for a degree-n polynomial is incredibly huge, which makes exhaustive search impossible as soon as n 6. Fortunately, heuristics exist to allow non-exhaustive search to higher degrees up to n ≈ 12 [144, ch.6].

Other schemes as well as the asymptotic optimality problem on bounded or unbounded parallelism have been studied in [START_REF] Kiyoshi | On the Parallel Evaluation of Polynomials[END_REF][START_REF] Munro | Optimal algorithms for parallel polynomial evaluation[END_REF][START_REF] Muraoka | Parallelism Exposure and Exploitation in Programs[END_REF]. Here, in the context of elementary function evaluation, and especially their vectorization, we performed various performance measurements on scalar and autovectorized classic schemes as well as on many generated schemes, that are presented in Section 3.3.

benchmarking polynomial evaluation schemes

In this section, benchmark results are analyzed for different evaluation schemes that were compiled with the GNU Compiler Collection The reciprocal throughput is a common metric for benchmarks. It corresponds to the mean number of cycles between two independent results: the lower, the better. In the following, we implicitly convert this metric by dealing with throughput: the higher, the better. Unless otherwise stated, this is the metric we use in all our benchmarks. For a given degree, it a priori depends upon the chosen scheme and the architecture, for which we alternate between enabling and disabling the generation of FMA instructions. On the targeted Haswell architecture, as long as no operation involves or produces subnormal numbers, we observe that the more ILP a scheme exhibits, the better its throughput. In Section 3.3.1, we perform and analyze throughput benchmarks of the three aforementioned "classic" schemes, making the polynomial degree vary from 3 to 32, as well as enabling or disabling compiler autovectorization and FMA instruction generation. We run the same benchmarks for both the binary32 and the binary64 precisions, which proves to be useful when trying to explain huge performance drops in the binary32 precision. Then, in Section 3.3.2, we run latency micro-benchmarks to gain confidence in our hypothesis that the interval on which a polynomial is designed to be evaluated can have an important impact on performance. Finally, in Section 3.3.3, we perform several benchmarks on Horner, Horner-2 and Estrin's schemes against all schemes of theoretical lowest latency on our simplified architecture for a degree-12 polynomial, and for various evaluation intervals. They confirm that the choice of an efficient polynomial evaluation scheme can greatly vary with the evaluation interval. They also show that on architectures like Haswell, when no subnormal number is produced during the evaluation, the less the number of multiplications involved, the greater the throughput.

Benchmarking classic schemes

In this section, in order to better represent the context of elementary functions evaluation, we show benchmarks for polynomials whose degree ranges between 3 and 32, and we use the Taylor coefficients for the function x → log(1 + x) around 0, computed with Sollya 6.0 [START_REF] Sylvain Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF]. This corresponds to the following Sollya command: > P = taylor(log1p(x), degree, 0); which generates a Taylor approximation of degree degree in the working precision. Then the polynomial coefficients are rounded to the target precision, e.g. binary32, with the following command: Figures 3.4 and 3.5 present the reciprocal throughputs measured for binary32 precision Horner, Horner-2 and Estrin schemes, with or without vectorization enabled and with or without FMA enabled. As expected, the higher the polynomial degree, the lower the throughput. However, starting from degree 16, we can see a huge performance loss for Estrin's scheme only. In Section 3.3.2, we explain this by the fact that the intermediate result x 16 , computed by Estrin's scheme only, produces subnormal numbers in approximately 13% of calls. Apart from this problem, we can notice that, without FMA, Estrin has the best throughput until the degree 15: For this degree, it is more than twice as good as Horner's scheme and 50% better than Horner-2's scheme, for both the scalar and vectorized versions. For the degree 16, when the subnormal production issue arises, Estrin's throughput is from 27 to 67% less than Horner's and from 49 to 79% less than Horner-2's, for the scalar and vectorized versions, respectively. With FMA enabled, the three schemes have comparable throughputs until degree 15, although the scalar version of Horner's scheme has the least throughput. Starting from degree 16, only Horner-2 is better than Horner, by 40%. Until degree 32, none of the three schemes seem to suffer from the limited number of registers. Indeed, GCC's scheduler manages to keep intermediate results in registers during each evaluation. This is an incentive for more use of highly parallel schemes in the context of elementary functions evaluation, for which the polynomial approximations are not likely to be of degrees higher than 32.

By replicating the same experiment but in binary64 precision, we do not observe the performance loss we had for Estrin's scheme in binary32 precision. This is clearly shown on Figures 3.6 and 3.7. This reinforces the hypothesis that the throughput loss was caused by the production of subnormal results. The throughput decrease is then almost linear with respect to the polynomial degree. Furthermore, Estrin's scheme has a greater throughput than Horner's starting from degree 10 in all cases: with or without FMA, with or without vectorization. We conclude that a lower latency does not always imply a throughput decrease, which can make schemes like second-order Horner or Estrin more interesting than Horner if performance is more important than precision. This conclusion may be a little counter-intuitive, since bounded parallelism can make us ex- pect some kind of saturation point on arithmetic units for high ILP evaluation schemes.

Influence of the evaluation interval

In the following, we note RN(x) the rounding of a real number x to the nearest (ties to even) binary32 floating-point number. The polynomials used for the evaluation of elementary functions are usually precise enough on rather small intervals. Inputs that lie outside of these intervals are either treated separately, or reduced to one of them. This interval is generally a neighborhood of the point of approximation of the elementary function, which is often 0. It is the case for ln(1 + x), exp(x), or sin(x). As a result, small magnitude inputs, if raised to certain powers, can rapidly produce subnormal numbers. For instance, for x ∈ [2 -12 ; 2 -5 ], which was the interval used in the benchmarks in the previous section, x 16 is a binary32 subnormal for about (2 -126/16 -2 -12 )/2 -5 ≈ 13% cases.

Figures 3.8 and 3.9 show the measured latency of two schemes of raising a binary32 floating-point number x in the aforementioned interval to the 16 th power. The first one corresponds to a naive exponentiation composed of 15 multiplications, while the other one uses a binary or "divide-and-conquer" (D&C) approach. Their minimal latencies are a few cycles long, except when x is small enough so that x 16 is a subnormal. Then a threshold can be observed: the measured latencies switch from a few cycles to more than 100. According to Agner Fog, these high cycle counts come from micro-coded special treatments for all operations whose operands are normalized floating-point numbers and whose result is a subnormal, as well as for multiplications with one normal and one subnormal operands [66, The naive evaluation scheme shows other thresholds, which correspond to monomials of lower degree that also produce subnormal numbers. The D&C approach gets back to having a latency of a few cycles after RN(x 16 ) = 0. This may contradict Fog's manual [66, § 10.9], which states: "There is no penalty for overflow, underflow, infinity or not-a-number results."

This suggests that the micro-code gets executed only for a precise interval of outputs. Also, when RN(x 15 ) = 0, i.e. |x| 2 -10 , a decreasing threshold can be observed for the naive scheme. This is interpreted by the fact that a multiplication by zero is not subject to this micro-code treatment penalty. Therefore, on some architectures, the fact that an intermediate result produces a subnormal can make computation performances decrease a lot. This conclusion will justify some algorithmic choices that we make in Chapter 4, when designing a high-throughput vectorized logarithm. In the context of elementary function evaluation, this performance issue can also affect polynomial evaluations. They are often valid approximations on intervals of small magnitudes, therefore higher-degree monomials can quickly produce subnormals. When efficiency is required it is thus relevant, if not necessary, to ensure that, for a given evaluation scheme, the input interval will not produce subnormal numbers. This conclusion could allow for better automatic selection of polynomial evaluation schemes, which would then be efficient on their whole interval of possible inputs. lowest-latency schemes, because it is able to produce fewer subnormal intermediate results, especially when the polynomial degree is not a power of two, like 12 here. On the other hand, Figure 3.12 shows that Estrin's scheme is about 12 times as less efficient as Horner and Horner-2 schemes when the production of subnormal results occurs starting from x 8 . Although Horner and Horner-2 schemes are more efficient than most lowest-latency schemes, there are about 50 schemes that stand better than them, among those that involve few multiplications. This goes in favor of our conclusion drawn in Section 3.3.2, i.e. that although not trivial, it is important to ensure that the scheme we choose does not produce subnormal results on its evaluation interval if we want to use it on architectures that have huge performance penalties for such results.

A deeper study of the lowest-latency schemes generated by CGPE may allow us to spot common patterns, that would be reusable in the context of elementary functions approximations by univariate polynomials. Otherwise, autotuning by benchmarking each scheme on a given architecture could allow us to select a scheme that has a better throughput than the usual schemes. A posteriori, bounds on the accuracy of the final result could be computed with Gappa to verify if the selected evaluation scheme meets the requirements.

conclusion and perspectives

We draw several conclusions about code generation for polynomial evaluations from the observations made in this chapter. In the context of elementary functions evaluation, we have shown that the production of subnormal results can have a huge impact on performance, which makes it a key factor when selecting an evaluation scheme that is efficient over all its evaluation interval. Finally, we have seen that highly parallel schemes, i.e. high ILP schemes, like Estrin or those generated by CGPE have often a better throughput than the more sequential schemes like Horner.

We believe that this study could be developed in three directions. First, a deeper analysis of the best schemes generated by CGPE may allow us to directly build efficient schemes without having to perform an expensive search. Second, CGPE could use new criteria to discard evaluation schemes that may produce unwanted kinds of intermediate or final results for a given input interval and a given set of coefficients. For example, smart interval arithmetic using the Gappa library could be used to implement such a heuristic. Finally, it might be interesting to study if an efficient evaluation scheme has a noticeable impact on the performance of elementary functions evaluations.

In the long term, an efficient code generator for polynomial evaluation in the context of elementary functions approximation would be desirable. Then it would be interesting to deal with performance/accuracy tradeoffs, as well as formal guarantees. We have seen several steps of the design of elementary functions implementations in the previous chapters. First, we developed a generator of exact lookup tables for trigonometric and hyperbolic functions (Chapter 2). Then, based on multiple benchmarks, we defined guidelines for the choice of efficient polynomial evaluation schemes. In order to generate efficient code, code generation frameworks should be able to leverage any source of performance enhancement, such as SIMD instruction sets extensions, e.g. SSE, AVX, or AVX-512, which provide support for vector instructions. Such instructions enable to treat packed inputs -currently up to 16 binary32 or 8 binary64 words for AVX-512 -in parallel, at the same announced latency and throughput as the equivalent scalar instruction. Hence, they allow for better sharing out of latency and multiplication of throughput compared to the same scalar instructions. In this chapter, we focus on automatic and efficient code generation for vectorized elementary functions. Such functions leverage the aforementioned SIMD capabilities to receive a short vector of inputs and return an output vector of same length. With the democratization of parallel computing units such as SIMD instruction sets or GPUs, this kind of functions have been gaining in popularity for efficient libm implementations [START_REF] Hauth | Development and Evaluation of Vectorised and Multi-Core Event Reconstruction Algorithms within the CMS Software Framework[END_REF][START_REF] Quirin | A new open-source SIMD vector libm fully implemented with high-level scalar C[END_REF][START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF]. This chapter presents a joint collaboration with Nicolas Brunie from Kalray. 1 He created and has been maintaining the Metalibmlugdunum code generation framework, which has been described in Chapter 1 (p.28). In particular, we present:

• a classic range reduction in a branchless fashion so as to use at best recent micro-architecture features, like the fast hardware reciprocal instruction rcp, and to treat all inputs in the same flow;

the power of meta-implementations

• a detailed "meta-design" of faithfully-rounded implementations of the natural logarithm;

• a meta-implementation within the Metalibm-lugdunum framework, including contributions to the SSE/AVX and AVX2 backends;

• how our meta-implementation allows to achieve high throughput implementations for the binary32 and binary64 formats in a fully automated way.

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 in a single instruction, thus potentially increasing the performance of floatingpoint programs. By taking into account these vector instructions during the design process, it is possible to achieve efficient implementations of vectorized mathematical functions. These are widely used in scientific programs manipulating floating-point computations, for which performance improvement is crucial [START_REF] Hauth | Development and Evaluation of Vectorised and Multi-Core Event Reconstruction Algorithms within the CMS Software Framework[END_REF].

Here, we deal with vectorized implementations of the natural logarithm function log(x) in floating-point arithmetic, with a particular focus on their automatic generation 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]. This framework enables to describe implementations of a function using a meta-language and to generate C source codes optimized for, e.g. different target architectures. For performance reasons, we target faithfully-rounded implementations, that is, with a maximum error of 1 ulp. Indeed, even if, the IEEE 754-2008 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][END_REF]. Furthermore, faithful rounding is often sufficient in high-performance computing contexts, for example in physics [START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[END_REF].

We provide 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 generate various 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. For performance reasons, we do not consider input arguments possibly requiring special treatments, e.g. NaNs, ±∞, or nonpositive numbers. However, we do consider subnormal inputs 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.

This chapter is organized as follows: Section 4.2 details the range reduction that we use to achieve the targeted accuracy. Section 4. 3 gives some guidelines on how to build accurate-enough programs. Then Section 4.4 provides some implementation details. Section 4.5 presents our meta-implementation within the Metalibm framework. Finally, Section 4.6 discusses some experimental results compared to existing implementations, and highlights advantages and disadvantages of meta-implementations, before a conclusion in Section 4.7.

algorithm for log(x) with faithful rounding

We consider here that the input x is neither a special input (NaN, ±∞, or x 0), nor the particular input x = 1. All floating-point numbers manipulated, as well as the floating-point operations carried out, are in precision p 2. This section presents the range reduction we use to implement the natural logarithm.

Range reduction

Let x = 1 be a positive floating-point number as defined in [START_REF][END_REF] with

s = 0: x = m • 2 e , (4.1)
where m ∈ [1, 2 -2

(x) = (e + τ) • log(2) -log(2 τ • r i ) + log(1 + u), where u = r i • m -1, and τ = [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, which provides a fast reciprocal approximation, 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 [START_REF] Daramy-Loirat | A library of correctly rounded elementary functions in double-precision[END_REF][START_REF] Revy | Automated Design of Floating-Point Logarithm Functions on Integer Processors[END_REF]. Hence let x be the rescaled floating-point number defined as:

x = x • 2 λ

Determination of u

Since u = r i • m -1, we deduce that u = ε in (4.5). Therefore the quantity u in (4.3) is the error occurring when approximating 1/m by r i . Property 1. Let x = 1 be a positive floating-point number as in (4.1). The quantity u in (4.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 . (4.6)
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 , 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 (4.5), 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 4.1, for different combinations (e max , κ). As a consequence, in order to implement the sequence computing u, we choose to use an fma, returning exactly u in a single instruction, and available on recent micro-architectures. If it had not been available, an alternative would have been to use multi-word arithmetic, like double-double arithmetic in [START_REF] Florent De Dinechin | Fast and correctly rounded logarithms in doubleprecision[END_REF][START_REF] Quirin | Arrondi correct de fonctions mathématiques -Fonctions univariées et bivariées, certification et automatisation[END_REF].

Determination of τ andlog(2 τ • r i )

The valuelog(2 τ • r i ) is an approximation oflog(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 oflog(1/m ), while for the others (τ = 1), it contains approximations oflog(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 follows:

τ = 1 if idx (2/ √ 2 -1) • 2 i , 0 otherwise. (4.7)
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. Notice that τ in (4.7) may not correspond exactly to the definition of τ in Section 4.2.1: but it 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.

How to ensure faithful rounding?

Let r be an approximation of log(x), with x = 1 a positive floatingpoint 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 | . (4.8)
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)). Following [START_REF] Muller | On the definition of ulp(x)[END_REF], a sufficient condition is:

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

|r h -log(x)| |(r h + r ) -log(x)| + |r | .
some coefficients of this polynomial, so that the special input 1 can also be handled in the general flow.

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 (4.11). 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 (4.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) is an approximation of log(2) stored as the unevaluated sum of two floating-point numbers:

L (2) = L (2) h + L (2) with L (2) -log(2) θ 1 .
Then log_tbl(i) returnslog(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 (4.5) 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)| |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 4.2.3, in order to ensure that the condition in (4.11) holds, we must build a program so that this sum verifies:

e + τ -λ • θ 1 + θ 2 + θ 3 + θ 4 θ. (4.12)
This is detailed in Section 4.4.2.

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.4), we know that

rcp(m ) ∈ [1 -κ, 1 + κ].
Recall that the quantity r i equals rcp(m ) plus one half-ulp 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 -κ, (4.13)
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 4 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 in order to handle the special input 1 in the general flow, i must be chosen to satisfy (4.13), and the polynomial approximant a(u) in Section 4.3.1 must be built so that the first coefficient (a 0 ) is zero, to ensure that a(u) = 0 in this case.

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.

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. Then the quantity λ in (4.2) is determined in a similar way, using comparison and logical instructions.

Second, once λ is known, we must compute x = x • 2 λ . We could have performed this by using a floating-point multiplication. However, as observed in [START_REF] De | 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 [76, § 3.4]. Let the following bit string encode a subnormal number x in binary32 format: 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 e min . This sequence of operations relies on arithmetic and logical instructions, but it requires also a means to determine λ, and more particularly a routine that returns λ if x is a subnormal number, and 0 otherwise. The instruction set we target does not embrace the nlz 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] Warren | 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.

Otherwise if

x is a subnormal number:

x = 2 -126 × 0. 000 • • • 000 λ 1 • • • .
Then at Line 24, we have:

Z = 2 0 × 1. 000 • • • 000 λ 1 • • • .
At Line 25, 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 (X >> 23) being zero, we deduce that mask = 2 32 -1 = 1111 • • • 1111, and value = -λ. Hence the routine returns λ.

This is apparently a well-known technique to compute the number of leading zeros in the bit string of a floating-point number. But to our knowledge, it has never been published, and was brought to our knowledge during a discussion with Christoph Lauter.

Certified evaluation program

Recall that in this section, we target the binary32 arithmetic, that is, for p = 24. For accuracy purpose, we must build an evaluation program, so that the bound in (4.12) is satisfied. Note that in the design process, as mentioned in Section 4.2.3, we must distinguish three cases. Indeed, the error bound θ = 2 -48 , which is the tightest of the three cases in (4.11), 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 θ 3 0, we know from (4.11) and (4.12) that the approximation error bound must satisfy θ 3 2 -48 . Remark 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 Sollya [START_REF] Sylvain Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF]. This way, we obtain θ 1 ≈ 2 -53. 33 . In addition Table 4.2 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 size i. In the rest of the example, we choose to implement log(x) using a table indexed by 7 bits and a degree-5 polynomial.

The remaining part is to compute the program evaluation error bound θ 4 , and to check if the bound (4.12) 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 [START_REF] Mouilleron | Automatic Generation of Fast and Certified Code for Polynomial Evaluation[END_REF], which enables to build several schemes to evaluate a given polynomial. Then, bounding the evaluation error is carried out using the Gappa tool [START_REF] Melquiond | De l'arithmétique d'intervalles à la certification de programmes[END_REF]. 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 [102, § 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 -50.86 , which is actually less than 2 -48 . Then according to the case, we obtain:

θ 4        2 -49.86 in Case 1, 2 -45.45 in Case 2, 2 -36.21 in Case 3,
and in all cases, the condition in (4.12) holds.

towards automated implementations

This section presents insights on how to automate the implementation of the logarithm function. It starts with a description of the Metalibm framework, before some details on its extension to our context.

The Metalibm-lugdunum framework

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 function libraries optimized for different architectures. The second kind is non-expert 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 transparently support different accuracies (improved, normal, degraded or customized), producing automatically different code versions,

• To factor common code across software and hardware architectures, for various I/O precisions and micro-architectures, thus reducing the development time,

• To vectorize meta-code according to different micro-architectures, 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 in the Cephes library [START_REF] Hauth | Development and Evaluation of Vectorised and Multi-Core Event Reconstruction Algorithms within the CMS Software Framework[END_REF][START_REF] Piparo | Speeding up HEP experiment software with a library of fast and auto-vectorisable mathematical functions[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 the possibility to use CGPE to build efficient polynomial evaluation schemes, as detailed in Section 4.5.2 below.

The case of polynomial evaluation schemes

As shown in Section 4.2.1, our implementations of log(x) function rely on the evaluation of a polynomial, whose degree varies from 4 up to 15 (Table 4.2). 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 on customizable architectures. More precisely, multiplier and adder latencies are customizable, but the use of heuristics has been automatized 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 thesis.

numerical results

Using the process presented before, we have generated different versions of log(x), in order to evaluate their performance. This section presents some numerical results. It compares our implementations to existing solutions. Then it studies the impact of the table index size and the cost of handling subnormal in the general flow.

Experimental protocol

We measure the reciprocal throughput of handmade and generated routines, and evaluate them against the Libmvec routines. All benchmarks are run on an Intel ® Xeon ® Processor E5-2650v4, which features the AVX2 instruction set extensions. The operating system is CentOS Linux release 7.4.1708, running a Linux kernel version 3.10.0-693.2.2.el7.x86_64. Although this CPU has 12 physical cores, we make sure that no other compute-intensive jobs are running at the same time as our benchmarks to reduce potential noise. All source code is compiled using GCC 7.2.0, linked against glibc 2.26, both compiled from source. Shared compiler flags include -O3 -mtune=native. For Libmvec benchmarks, we also have 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 use the -march=core-avx2 flag, while we use -march=corei7 to restrict GCC to emitting SSE4 code.

To evaluate these performances, we use either automated microbenchmarks provided by Metalibm or custom micro-benchmarks for the handmade version and Libmvec routines. These custom microbenchmarks try to reproduce plausible workloads by using pseudorandom inputs, but also rare workloads by using only subnormal inputs. The benchmarks are designed to warm up the L1 data and instruction caches before the main measurement, so that cache-miss penalties may be kept to a minimum. Also, we choose 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 benchmark separately pseudorandom input vectors versus subnormals-only vectors, which may yield lower performance. We also benchmark constants-only vectors, which might yield better throughputs for vectorized memory gather. We observe that using our approach, we achieve a binary32 implementation with a reciprocal throughput of 20.7 CPE for vector size = 4 and 11.4 CPE for vector size = 8. In binary64, the measured performances are 88.2 CPE for vector size = 2 and 44.9 CPE for vector size = 4. This is greater than binary32 versions by a factor close to 4. Obviously, our implementations are slower than Libmvec implementations: for example, in binary32 their reciprocal throughputs vary from 0.3 up to 6.7 CPE with respect 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 (binary32), while in our case, polynomials are evaluated using double-single arithmetic. A direct consequence is the accuracy of the output. Indeed, our implementations are correct within 1 ulp (i.e. faithful), while Libmvec provides functions with as much as 4 ulp absolute error. 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 less efficient when dealing with a subnormalsonly 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.)

An interesting observation is that the multiple-indices vectorized accesses to tabulated values are not penalizing, although data is accessed at non-contiguous locations. This may be explained by the fact that all tables are small -typically 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 targeted processor.

Note finally that, as mentioned above, our routines are guaranteed faithfully rounded. However, to gain in confidence, the binary32 versions have been exhaustively verified and compared with the values returned by the MPFR library.

.6.3 Impact of the table size and subnormal handling

In this section, we measure the impact on performance of the table size and the cost of supporting subnormals in the main flow. To do this, we generate routines in the binary32 format with the Metalibm framework, for a parameter i ranging from 3 up to 9, for different vector sizes (typically 4 and 8), and with and without subnormal handling. To unplug subnormal handling, it suffices to remove λ from all the computations. Hence we simply adapt our meta-code in this sense. Table 4.4 shows the performance of these routines in CPE.

As expected, the greater the table size, the better the routines performance. Using this table we may conclude that, for performance purpose, a table indexed by 8 bits seems to be a good choice.

Furthermore, an interesting observation is that treating subnormal floating-point numbers in the main flow has no great impact on the performance of the whole function. Indeed, except for i = 5 and vector size = 4, the overhead due to this technique remains no greater than 2.3 CPE, which is acceptable. This reinforces our choice to avoid branches and fallback routines to treat subnormals separately.

conclusion and perspectives

In this chapter we have addressed the implementation of logarithm functions on vector micro-architectures, with a particular focus on its automation through the Metalibm framework. We achieve high throughput implementations with 1-ulp accuracy, optimized for SSE, AVX and AVX2 micro-architectures. For example, on AVX2 microarchitectures, we measure throughputs of 11.4 and 44.9 CPE for bi-nary32 and binary64, respectively, with a relatively small overhead due to handling subnormals in the main flow (≈ 2 CPE).

In Table 4.3, "n/a" indicates that Metalibm was not able to generate these functions. This is actually due to the lack of certain meta-instructions in the x86 backend. Ongoing work focuses on integrating these requirements. In a near future, this could enable automatic generation of binary64 implementations as well. More generally, the IEEE 754 standard defines two other formats, namely, the binary16 and binary128. Since our design does not depend upon the precision, generating implementations for these formats is reachable and of interest. This would require efforts to eventually emulate the underlying required missing arithmetic (SIMD support for binary128, for example), and to adapt the Metalibm backend in consequence.

Furthermore 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 highperformance computing.

"This is it. Let me give you one piece of advice: be honest. He knows more than you can imagine."

Trinity, in The Matrix, by the Wachowski brothers.

This thesis aimed at developing more automatic code generation algorithms and techniques for the efficient evaluation of elementary functions. In Chapter 2 we introduced a novel and elegant method to automatically generate exact lookup tables for trigonometric and hyperbolic functions, a goal that had never been reached for decades. This technique could overtake improvements such as those that were brought by Gal and Bachelis or Stehlé and Zimmermann at least for the trigonometric and hyperbolic functions.

In Chapter 3 we analyzed the impact on throughput brought by vectorizing polynomial evaluation schemes. First, we confirmed previous results showing that elementary arithmetic operations producing subnormal floating-point numbers could drastically reduce the throughput of applications such as polynomial evaluation schemes. Second, contrary to the generally assumed idea, we observed that higher-ILP schemes could increase the throughput on Intel Haswell architectures.

Finally, in Chapter 4, we advocated the concept of "metalibms", i.e. frameworks that allow to describe rather than to implement elementary functions or evaluation schemes for elementary functions, and that are then able to generate source codes for different optimization criteria. For example, such criteria include precision, accuracy, target architecture or vectorization. Through our meta-implementation of the natural logarithm, we were able to describe our contributions to Metalibm-lugdunum, a framework for developers and users of mathematical libraries that allows to describe "meta-evaluation schemes" and to generate optimized source code for various targets. Our contributions included improved support for the Intel x86 backend and its vectorized instructions, so that we were able to generate fullyvectorized logarithms in the binary32 and the binary64 floating-point precisions for portable C code and the Intel AVX2 architecture. One of the main advantages of this framework is that we only had to describe one meta-algorithm to be able to generate many optimized implementations, which can bring huge time-saving opportunities in the future.

The AVX2 versions of our meta-logarithm were measured to be competitive on several aspects compared to the GNU Libmvec implementations: First, although the automatically generated codes showed lower measured throughputs than GNU Libmvec, our meta-imple-mentation is proven faithfully-rounded while the GNU Libmvec is only accurate to a tested 4 ulps accuracy. Second, the GNU Libmvec consists in several routines written in x86 assembly, whereas our meta-implementation can generate multiple source codes for several languages and targets. For our example, we used the contributed x86 backend to generate C codes that use the Intel Intrinsics C API.

new horizons

In Chapter 2, we presented an exact table lookup generation heuristic using Fässler's algorithms, for the trigonometric functions only. When targeting correct rounding in double precision, we achieved to shrink lookup tables by up to 29%, compared to Tang's and Gal's methods. We mentioned that this benefit could even be higher for greater targeted accuracies since all tabulated values except the corrective terms are exactly stored and, thus, do not depend on the targeted accuracy. Also, we only focused on the trigonometric and hyperbolic functions, as they both benefit from the use of Pythagorean triples. However, the concept remains valid and worthwhile for other functions, provided at least two terms need to be tabulated per row and that exact values can be found for such terms. Therefore, shortto-medium term interesting perspectives could be:

• the implementation of a heuristic for the hyperbolic functions, by generating Pythagorean triples sharing a common leg. A first lead might be to follow the results presented in [START_REF] Rothbart | Pythagorean Triples: A New, Easy-to-Derive Formula With Some Geometric Applications[END_REF] by Rothbart and Paulsell;

• to investigate sufficient conditions for the optimality of such heuristics, as a first step towards a more formal validation of our approach;

• the integration and benchmarking of our exact lookup tables, e.g. into the high performance implementations of the trigonometric functions in the medium-precision range (4096 bits) detailed in [START_REF] Johansson | Efficient Implementation of Elementary Functions in the Medium-Precision Range[END_REF]. Indeed, an actual proof-of-concept using our method would probably contribute much to its spreading in real-world code bases.

Also, a short-term objective could be to integrate our tables to the Metalibm-lugdunum framework, so as to enable metaimplementations to use them for accurate evaluation schemes;

• to investigate whether other elementary functions could benefit from such exact lookup tables. The sinPi and cosPi routines might be a first lead, although we do not foresee much more than that. Indeed, we must insist that our method is only interesting whenever two or more values must be tabulated. Otherwise, Gal's accurate tables are probably a better option.

As far as evaluating polynomials is concerned, we believe that the study performed in Chapter 3 could be developed in several directions, from short-term to long-term perspectives:

• first, it might be interesting to study if an efficient evaluation scheme has a noticeable impact on the performance of elementary functions evaluations.

Early tests on classic function evaluation schemes tend to show negligible impacts. However, further investigations, e.g. within the Metalibm-lutetia open-ended generation framework, might demonstrate the contrary, because it relies heavily on a domainsplitting algorithm [START_REF] Brunie | Code Generators for Mathematical Functions[END_REF];

• next, CGPE could be developed to use new criteria to discard polynomial evaluation schemes. For instance, schemes that may produce subnormal floating-point results for a given input interval and a given set of coefficients. The Gappa tool could probably be used to implement such a heuristic quite efficiently. As a result, CGPE would be able to certify such property when targeting high performance polynomial evaluation;

• a medium to long-term goal would be to merge CGPE features into the different Metalibm projects so as to benefit from an efficient code generator for polynomial evaluations. Such a project has already started within the Metalibm-lugdunum framework, through the PythonCGPE bindings presented in Chapter 4;

• finally, a long-term research opportunity would be to do an in-depth analysis of the lowest-latency schemes generated by CGPE. It may allow us to directly build efficient schemes without having to perform an exponentially expensive search. This study should probably be combined with theoretical research on "polynomial chains" [START_REF] Knuth | The Art of Computer Programming[END_REF].

The paradigm shift towards more meta-and generative programming presented in Chapter 1, Section 1.3 and in Chapter 4, looks very promising to produce tailored mathematical libraries for high performance software. In addition, it could reduce libms maintenance costs and development time by pushing for generic code reuse. This is why we believe that, along with other complementary concepts of "metalibms" such as MetaLibm-lutetia, meta-implementations are the modern way to implement, optimize, maintain, and generate elementary functions evaluation schemes. Hence, we think that potentially interesting future research topics could include:

• merging similar meta-implementations, e.g. multiple logarithms into an efficient, generic, base-β logarithm. This should be a rather straightforward application of [START_REF] Revy | Automated Design of Floating-Point Logarithm Functions on Integer Processors[END_REF] in the context of Metalibm-lugdunum;

• providing support for other IEEE 754 formats such as the bi-nary16 and binary128 formats, or decimal floating-point for-mats. Indeed, these formats may be relevant choices for specific applications. Metalibm-lugdunum already supports several mixed-precision FMA operators for binary16 and low precision fixed-point implementations, both useful for e.g. machine learning applications [START_REF] Brunie | Modified Fused Multiply and Add for Exact Low Precision Product Accumulation[END_REF]. In this sense, it could be interesting to list and implement what is missing in the different MetaLibm projects to support code generation for these formats;

• making an inventory of generic, reusable meta-blocks. For instance, in Metalibm-lugdunum, there is a somewhat generic Payne and Hanek range reduction [START_REF] Payne | Radian Reduction for Trigonometric Functions[END_REF], that could be further abstracted.

Also, following the specific hardware support for the fast hardware reciprocal approximation that we implemented in Metalibm-lugdunum, we could provide new meta-blocks or optimization passes for specific hardware features, e.g. for NVIDIA Tensor Cores [START_REF] Nvidia Tesla | GPU architecture. the world's most advanced data center GPU[END_REF]; These units basically provide floating-point 4 × 4 × 4 matrix-matrix FMA in mixed-precision. 1 According to NVIDIA, Tensor Cores can increase throughput for GEMM 2 by a factor up to 12x compared to GP100 accelerators. Therefore, more support of such hardware specificities through higherlevel meta-code and optimization passes within the Metalibmlugdunum framework could bring decisive advantages over existing solutions like CUDA programming. Metalibm-lugdunum already has an OpenCL backend, hence it should not be difficult to develop a CUDA backend too. Then, BLAS 3 could also benefit from tailored code generation within a unified code generation framework.

• finally, Metalibm-lutetia has been supporting automatic blackbox detection of periodicity and symmetry properties [START_REF] Brunie | Code Generators for Mathematical Functions[END_REF]. A medium-term research problem could be to provide support for automated asymptotic analysis, too. First, a white-box approach seems more realistic as it could reason on algebraic properties. Such analysis might allow for an improved domainsplitting algorithm whenever clear asymptotic behaviors are detected.

On retrouve des fonctions mathématiques élémentaires dans de nombreux codes de calcul hautes performances. Or, bien que les bibliothèques de fonctions mathématiques (libm) auxquelles font appel ces codes proposent en général plusieurs variétés d'une même fonction, ces dernières sont figées lors de leur implémentation. Cette caractéristique monolithique des libms représente donc un frein à la performance des programmes qui les utilisent, car elles sont conçues pour être polyvalentes au détriment d'optimisations spécifiques. De plus, la duplication de modèles partagés dans le code rend la maintenance de ces libms plus difficile et sujette à l'introduction de bugs. Un défi actuel est de proposer des "méta-outils" visant la génération automatique de codes performants pour l'évaluation des fonctions élémentaires. Ces outils doivent permettre la réutilisation d'algorithmes efficaces et génériques pour différentes variétés de fonctions ou architectures matérielles. Il devient alors possible de générer des libms optimisées pour des besoins très spécifiques avec du code générateur factorisé, ce qui facilite sa maintenance. Dans un premier temps, nous proposons un algorithme original permettant de générer des tables de correspondances sans erreur d'arrondi pour les fonctions trigonométriques et hyperboliques. Puis, nous étudions les performances de schémas d'évaluation polynomiale vectorisés, premier pas vers la génération de fonctions vectorisées efficaces. Enfin, nous proposons une méta-implémentation d'un logarithme vectorisé, factorisant la génération de code pour différents formats et architectures. Nous montrons que ces contributions sont compétitives comparées à des solutions libres ou commerciales, justifiant le développement de ce nouveau paradigme.

mots-clés Fonctions élémentaires, performance, génération de code, évaluation polynomiale.

Elementary mathematical functions are pervasive in many high performance computing programs. However, although the mathematical libraries (libms), on which these programs rely, generally provide several flavors of the same function, these are fixed at implementation time. Hence this monolithic characteristic of libms is an obstacle for the performance of programs relying on them, because they are designed to be versatile at the expense of specific optimizations. Moreover, the duplication of shared patterns in the source code makes maintaining such code bases more error prone and difficult. A current challenge is to propose "meta-tools" targeting automated high performance code generation for the evaluation of elementary functions. These tools must allow reuse of generic and efficient algorithms for different flavours of functions or hardware architectures. Then, it becomes possible to generate optimized tailored libms with factorized generative code, which eases its maintenance. First, we propose an novel algorithm that allows to generate lookup tables that remove rounding errors for trigonometric and hyperbolic functions. The, we study the performance of vectorized polynomial evaluation schemes, a first step towards the generation of efficient vectorized elementary functions. Finally, we develop a meta-implementation of a vectorized logarithm, which factors code generation for different formats and architectures. Our contributions are shown competitive compared to free or commercial solutions, which is a strong incentive to push for developing this new paradigm.

keywords Elementary functions, performance, code generation, polynomial evaluation.

0. 1

 1 arithmétique flottante xi Gappa, pour formaliser et certifier le calcul d'erreur

A

  = (a/c) • k et B = (b/c) • k, avec k ∈ N * et k unique pour toute la

0. 3 . 5

 35 Vers des schémas plus performants L'outil CGPE 1 permet de calculer tous les schémas d'évaluation polynomiale de plus faible latence sur parallélisme infini, pour un degré et une architecture donnés. Nous l'avons utilisé pour générer ceux de degré 12, ce qui représente 11 944 schémas de latence 26 avec l'architecture considérée. Les débits de ces schémas compilés sans FMA et avec autovectorisation, ont ensuite été mesurés pour des entrées dont certaines puissances dénormalisent. Les Figure 3.10 (p.83), 3.11 (p.84), et 3.12 (p.84) représentent ces 11 944 mesures par des points. Les traits verticaux distinguent les schémas selon le nombre de multiplications (de 17 à 25) qu'ils impliquent.

1 . 2 .

 12 They are not versatile, i.e. they can only be used on fixed-length inputs [78, § 8.3.7]; They are not very accurate on their whole domain [78, § 8.3.10];

1 C

 1 O M P U TAT I O N O F E L E M E N TA R Y F U N C T I O N S "You take the red pill: you stay in Wonderland, and I show you how deep the rabbit-hole goes." Morpheus, in The Matrix, by the Wachowski brothers.

Figure 1 . 1 -

 11 Figure 1.1 -IEEE 754 Binary Interchange Floating-Point Format.

Figure 1 . 2 -

 12 Figure 1.2 -Effect of IEEE 754 rounding-direction attributes for positive results.

  function-specific identities These are often used for range reductions and may overlap above properties. For example: log(ab) = log a + log b; sin(a + b) = sin a cos b + sin b cos a; cos x + (2n + 1) π 2 = (-1) n sin x; cosh(a + b) = cosh a cosh b + sinh a sinh b;

@

  rnd = float<ieee _ 32, ne>; x = rnd(Mx); { Mx in [-0.1,0.1] -> rnd(x * x) in ? } 1 http://cgpe.gforge.inria.fr/ 2 Génération Automatique de Preuves de Propriétés Arithmétiques [automatic proof generation of arithmetic properties]. http://gappa.gforge.inria.fr/ computation of elementary functions Gappa then returns Results: float<24,-149,ne>(x * x) in [0, 10737419b-30 {0.01, 2^(-6.64385)}]

  carg = GET _ VEC _ FIELD _ ADDR(vec _ x); tmp = _ mm256 _ load _ ps(carg); y = _ mm256 _ mul _ ps(tmp, tmp);

Figure 2 . 2 -

 22 Figure 2.2 -Visual representation of the second range reduction with a regularly-spaced trigonometric table.

Figure 2 . 3 -

 23 Figure 2.3 -Visual representation of the second range reduction with an example of Gal's trigonometric table.

  Hence, for each Pythagorean triple (a, b, c), assuming b = 0, ∃θ ∈ [0, π/2[ , cos(θ) = b c and sin(θ) = a c ∃ϕ ∈ R 0 , cosh(ϕ) = c b and sinh(ϕ) = a b .

Figure 2 . 4 -

 24 Figure 2.4 -Visual representation of the Pythagorean triple (3, 4, 5).

4 Figure 2 . 5 -

 425 Figure 2.5 -Primitive Pythagorean triples with a hypotenuse c < 2 12 .

Figure 2 . 6 -

 26 Figure 2.6 -Visual representation of PPTs and their symmetric counterparts falling into a trigonometric exact lookup table.

1 .

 1 For the trigonometric functions, Figure 2.6 illustrates the following advantages: (a) If arcsin(b/c) > π/4, one has arcsin(a/c) < π/4, which falls into the range of the exact lookup table.

  only the symmetric PPT (b, a, c) falls into the range of the lookup table.

Figure 2 . 8 -

 28 Figure 2.8 -Number of PPTs per entry of a trigonometric table for p = 7 with a hypotenuse less than 2 14 . Only the first entry does not contain any PPT.

4 Figure 2 . 9 -

 429 Figure 2.9 -Visual representation of the trigonometric exact lookup table from table 2.3 in the Euclidean plane.

- 6 :

 6 sin(x -corr)/k ≈ 1.01000110101 2 • 2 -18 cos(x -corr)/k ≈ 1.011110111 2 • 2 -13 .

Figure 3 . 1 -

 31 Figure 3.1 -Horner's scheme for a degree-3 polynomial.

Figure 3 . 2 -Figure 3 . 3 -

 3233 Figure 3.2 -Second-order Horner's scheme for a degree-3 polynomial.

Figure 3 . 4 -

 34 Figure 3.4 -Means of reciprocal throughputs of three usual schemes in binary32 precision, without FMA.

>

  P _ binary32 = roundcoefficients(P, [|single...|]); Several evaluations are benchmarked on large vectors of 2 14 floatingpoint numbers uniformly distributed in the interval [2 -12 ; 2 -5 ].

Figure 3 . 5 -

 35 Figure 3.5 -Means of reciprocal throughputs of three usual schemes in binary32 precision, with FMA.

Figure 3 . 6 -

 36 Figure 3.6 -Means of reciprocal throughputs of three usual schemes in binary64 precision, without FMA.

Figure 3 . 7 -

 37 Figure 3.7 -Means of reciprocal throughputs of three usual schemes in binary64 precision, with FMA.

Figure 3 . 8 -

 38 Figure 3.8 -Latencies of the computation of x 16 in binary32 by naive exponentiations.

Figure 3 . 9 -

 39 Figure 3.9 -Latencies of the computation of x 16 in binary32 by binary exponentiations.

Figure 3 . 11 -

 311 Schemes sorted by increasing number of multiplications

  Schemes sorted by increasing number of multiplications

Figure 3 . 12 -

 312 Figure 3.12 -Mean reciprocal throughputs of lowest-latency schemes of degree 12, sorted by increasing number of multiplications. x n is subnormal when n = 8.

T

  H E P O W E R O F M E TA -I M P L E M E N TAT I O N S "I've done this because, philosophically, I am sympathetic to your aim. I can tell you with no ego, this is my finest sword. If on your journey, you should encounter God, God will be cut." Hattori Hanzo, in Kill Bill: Vol.1 by Quentin Tarantino.

(4. 2 )

 2 where λ ∈ {0, • • • , p -1} is the number of leading zeros in the binary representation of m. It follows that x = 2 e • m 2 e min and log(x) = (e + τ -λ) • log(2) -log(2 τ • r i ) + log(1 + u), (4.3) with m ∈ [1, 2 -2 1-p ], represented exactly with only one precisionp floating-point number. The quantity r i is thus an approximation of 1/m obtained with the rcp instruction truncated on i bits, and

Listing 4 . 1 -

 41 Vectorizable Algorithm to Compute λ.

1 .

 1 typedef union { uint32 _ t i; float f; } cfloat32; uint32 _ t lambda _ or _ zero(uint32 _ t X) { cfloat32 Z;Z.i = X | 0x3f800000; // Exponent mask of 1.0f Z.f -= 1.f; uint32 _ t mask = 0xffffffff + ((X >> 23) != 0); uint32 _ t value = ((Z.i >> 23) -127) & mask; return -value; } If x is a normal number, its exponent field (X >> 23) is nonzero,and mask = 0. Hence the routine returns 0.

1

 1

4

 4

  et 2 -150 . On choisit p = 10 pour notre table, ce qui correspond à π/4 × 2 10 = 805 entrées. Pour faciliter la comparaison, on ne considère que le nombre d'accès mémoire (MA) requis par la seconde réduction d'argument et le nombre d'opérations flottantes (FLOP) effectuées pendant l'étape de reconstruction. Aussi, on considère que des algorithmes de calcul en expansions sont utilisés quand une grande précision est requise, comme c'est le cas par exemple dans la bibliothèque CR-Libm. La Table 2.10 (p.65), extraite de [34, § 2.3, 101, § 3.2, 102], sera utilisée comme référence pour calculer le coût de ces algorithmes quand aucune instruction FMA (Fused-Multiply and Add) n'est disponible dans le jeu d'instructions. La notation E n désigne une expansion de taille n, c'est-à-dire, un nombre flottant représenté en machine comme la somme non évaluée de n nombres flottants IEEE 754. Avec ce formalisme, E 1 représente un nombre flottant habituel. Le coût de ces opérations élémentaires en précision étendue est égal au nombre d'instructions flottantes nécessaires pour les réaliser. Nos estimations sont résumées dans la Table 2.11 (p.67), qui contient le nombre d'opérations flottantes et d'accès mémoire pour la phase rapide et la phase précise ainsi que l'occupation mémoire de chaque table en octets/ligne. On remarque tout d'abord que la réduction d'argument proposée requiert moins de mémoire par entrée dans la table que les autres solutions. Il faut en effet 48 octets par entrée pour la méthode de Tang et 56 octets par entrée pour celle de Gal contre seulement 40 octets par entrée pour la table de points exacts. Cela représente un gain en mémoire d'environ 17% et 29% par rapport à Tang et Gal, respectivement.

	xix
	Comparaison des résultats

  minimise le nombre de multiplications, n'a pas besoin de stocker de résultat intermédiaire[90, p. 486] et il est stable numériquement[14, ch.9]. Quand le critère de choix est la précision, Horner est donc un bon schéma. Mais, on l'a vu, sa latence est élevée. Or, il existe des schémas plus rapides. Par exemple, le schéma de Horner d'ordre 2, noté Horner-2, a le parenthésage suivant au degré 7 :

1 https://software.intel.com/haswell 0.3 vectorisation de schémas polynomiaux xxi Son intérêt est triple : il

  3.2 (p.76) : les mônomes de degrés pair et impair peuvent être évalués en parallèle. Sa latence est donc environ divisée par deux [90, p. 488]. Le schéma d'Estrin, illustré en Figure 3.3 (p.76), expose encore plus d'ILP en découpant l'arbre d'évaluation selon les puissances de x de la forme 2 i [61]. Par exemple, pour n = 7, on aura :

  1 . Nous réalisons plusieurs évaluations sur des vecteurs de 214 nombres flottants uniformément répartis dans l'intervalle [2 -12 ; 2 -5 ]. Les Figures 3.4 (p.78) et 3.5 (p.79) présentent les débits inverses mesurés en simple précision pour Horner, Horner-2 et Estrin, avec ou sans vectorisation, avec ou sans FMA, respectivement. Sans surprise, le débit diminue avec le degré pour tous les schémas. En revanche, on observe une perte de performance drastique à partir du degré 16, uniquement pour le schéma d'Estrin. Dans la section suivante, nous expliquons cette perte de performance par le fait que le résultat intermédiaire x 16 , calculé par Estrin, dénormalise dans 13% des cas. Sans FMA, Estrin a le meilleur débit jusqu'au degré 15 : pour ce degré, il est plus de deux fois meilleur que Horner et 50% meilleur que Horner-2, pour les deux versions (scalaire et vectorisée). Au degré 16, il est de 27 à 67% inférieur à Horner et de 49 à 79% inférieur à Horner-2, pour les versions scalaire et vectorisée, respectivement. Avec FMA, les trois schémas se valent jusqu'au degré 15, avec un désavantage pour Horner non vectorisé. À partir du degré 16, seul Horner-2 est meilleur que Horner, d'environ 40%. Globalement, jusqu'au degré 32,

  Ces résultats sont présentés dans laTable 4.3 (p.101). Également, nous avons mesuré l'impact du traitement des dénormalisés et de la taille de notre table sur le débit de 28 implémentations générées automatiquement. Ces mesures sont reportées dans la Table 4.4 (p.103).

	xxviii	résumé substantiel en français
	le schéma d'Estrin. Afin de générer d'autres schémas potentiellement
	plus efficaces automatiquement, nous avons développé et intégré une
	interface entre CGPE et Metalibm-lugdunum nommée PythonCGPE.
	Comme CGPE est écrit en C++ et Metalibm en Python, notre choix
	s'est porté sur Cython pour l'implémentation de cette interface.
	0.4.4 Mesures de performances du code généré
	Nous avons mesuré le débit inverse de plusieurs routines vectorisées,
	soit écrites manuellement, soit générées automatiquement, ainsi que
	plusieurs routines de la GNU Libmvec présente dans glibc 2.26.
	où κ est une borne d'erreur relative sur l'approximation de 1/m
	par l'instruction rcp. C'est pourquoi nous privilégions la génération
	d'une instruction fma pour calculer u exactement, disponible sur la
	plupart des architectures récentes. Sans cette possibilité, une solution
	est de générer une séquence de calcul en arithmétique multi-mots
	comme par exemple la précision double-double [53, 101].
	Nous montrons finalement comment garantir l'arrondi fidèle du
	résultat r h en distinguant trois cas qui nous fournissent trois bornes
	d'erreur sur |log(x) -r| où r = r h + r l est le résultat intermédiaire
	sur deux mots de précision p. Ces bornes nous ont ensuite permis de
	certifier que notre méta-implémentation renvoie l'arrondi fidèle du
	logarithme de x.
	0.4.3 Méta-implémentation dans Metalibm-lugdunum
	Dans cette section, nous montrons comment calculer x en utilisant
	uniquement des instructions en arithmétique entière, permettant ainsi
	d'éviter les surcouts liés au traitement des nombres dénormalisés en
	arithmétique flottante. Le Listing 4.1 (p.96) présente ainsi un algo-
	rithme efficace et vectorisable qui calcule λ avec uniquement des opé-
	rations entières. Le calcul de x = x • 2 λ revient ensuite simplement
	à décaler x de λ bits vers la gauche, qui est une opération rapide et
	vectorisable.
	Un tel algorithme a été implémenté sous la forme d'un "méta-
	bloc" dans Metalibm-lugdunum. Metalibm-lugdunum est un envi-
	ronnement de développement permettant de générer rapidement du
	code logiciel ou matériel efficace, principalement pour les fonctions
	élémentaires et le calcul numérique. Il repose sur un langage spéci-
	fique capable d'abstraire des concepts comme la cible matérielle ou la
	précision, entre autres. Sur ce principe, nous avons ajouté du support
	pour le générateur de code ciblant les architectures Intel AVX2, afin
	que Metalibm puisse générer du code vectorisé utilisant les instruc-
	tions vrcpps ou vgatherdps, par exemple.
	En ce qui concerne l'évaluation polynomiale, Metalibm-lugdunum
	propose un générateur de schéma de Horner ainsi qu'un autre pour

Table 2 . 6

 26 Tables Exactes pour les Fonctions Élémentaires . . . . . . . xii 0.2.1 Évaluation à base de tables . . . . . . . . . . . . . . . . . xii 0.2.2 Tabulation de termes exacts . . . . . . . . . . . . . . . . . xiv 0.2.3 Implantation et résultats numériques . . . . . . . . . . . xvi 0.2.4 Comparaisons avec les méthodes existantes . . . . . . . xviii 0.2.5 Conclusions et perspectives . . . . . . . . . . . . . . . . . xix 0.3 Vectorisation de Schémas Polynomiaux . . . . . . . . . . . xix 0.3.1 Contexte et Motivations . . . . . . . . . . . . . . . . . . . xx 0.3.2 Rappels sur l'évaluation polynomiale . . . . . . . . . . Mesures de performances du code généré . . . . . . . .xxviii 0.4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . .xxviii Conclusion Générale . . . . . . . . . . . . . . . . . . . . . . . . .xxviii Evaluation of Elementary Functions . . . . . . . . . . . . . 14 1.2.1 The Table Maker's Dilemma . . . . . . . . . . . . . . . . 16 1.2.2 Algorithms for the evaluation of elementary functions . 18 1.2.3 SIMD vectorization of elementary functions . . . . . . . 23 1.3 Code Generation for Elementary Functions . . . . . . . . . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.1.1 Overview of this chapter . . . . . . . . . . . . . . . . . . . 32 2.2 Table-Based Evaluation Schemes . . . . . . . . . . . . . . . 33 2.2.1 Usual framework . . . . . . . . . . . . . . . . . . . . . . . 33 2.2.2 Focus on table-lookup range reductions . . . . . . . . . . 35 2.3 Design of an Exact Lookup Table . . . . . . . . . . . . . . . 40 2.4 Application to Sine and Cosine Functions . . . . . . . . . . 42 2.4.1 Pythagorean triples and their properties . . . . . . . . . 42 2.4.2 Construction of primitive Pythagorean triples . . . . . . 43 2.4.3 Selection of primitive Pythagorean triples . . . . . . . . . 47 2.5 Implementations and Numeric Results . . . . . . . . . . . 51 2.5.1 Exhaustive search . . . . . . . . . . . . . . . . . . . . . . . 52 2.5.2 Characterizing exhaustive results: a first heuristic . . . . 57 2.5.3 A much faster heuristic using Fässler's algorithms . . . 59 2.6 Comparisons with Other Methods . . . . . . . . . . . . . . Range reduction . . . . . . . . . . . . . . . . . . . . . . . 89 4.2.2 How to determine r i , u, τ, andlog(2 τ • r i )? . . . . . . 90 4.2.3 How to ensure faithful rounding? . . . . . . . . . . . . . 92 4.3 Guidelines to Build Accurate Enough Programs . . . . . . 93 Program to handle the general flow . . . . . . . . . . . . 94 4.3.2 How to handle the special input 1 in the general flow? . 95 4.4 Implementation Details . . . . . . . . . . . . . . . . . . . . 95 4.4.1 Branchless program . . . . . . . . . . . . . . . . . . . . . 95 4.4.2 Certified evaluation program . . . . . . . . . . . . . . . . 97 4.5 Towards Automated Implementations . . . . . . . . . . . . 98 4.5.1 The Metalibm-lugdunum framework . . . . . . . . . . . 98 4.5.2 The case of polynomial evaluation schemes . . . . . . . . 99 4.6 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . 100 4.6.1 Experimental protocol . . . . . . . . . . . . . . . . . . . . 100 4.6.2 Performance of the log(x) function . . . . . . . . . . . . . 102 4.6.3 Impact of the table size and subnormal handling . . . . 102 4.7 Conclusion and Perspectives . . . . . . . . . . . . . . . . . 103 Prime

	C O N T E N T S
	Dédicace

Factorization of Found Common Multiples for sinh and cosh. . . . . . . . . . . . . . . Table 2.7 Heuristic Search Results for sin and cos. . . . . Table 2.8 Heuristic Search Results for sinh and cosh. . . Table 2.9 "Fässler's" Heuristic Search Results for sin and cos. . . . . . . . . . . . . . . . . . . . . . . . . . Table 2.10 Costs of Addition and Multiplication of Expansions. . . . . . . . . . . . . . . . . . . . . . . . . Table 2.11 Computation and Memory Costs for Three Table-Based Range Reductions for Trigonometric and Hyperbolic Functions. . . . . . . . . . . . . . . Table 2.12 Tang's Table for p = 4 With 16-Bit Precision. . Table 3.1 Theoretical latencies for common schemes on a simplified Haswell architecture. . . . . . . . . Table 4.1 Suitable Values for i (Satisfying Property 1) for Different Combinations of (e max , κ). . . . . . . Table 4.2 Implementation Parameters and Table and Polynomial Coefficient Sizes (# Bytes) for Various Values of i When p = 24. . . . . . . . . . . . .

Table 4 . 3

 43 Measured Performances of our Implementations Compared to the GNU Libmvec. . . . . . . . .

	Table 4.4	Measured Performances (CPE) on AVX2 of our
		Generated Implementations, in Binary32 Arith-
		metic and for Vectors of Size 4 (v4) or 8 (v8),
		With (sub) or Without (no-sub) Subnormal Han-

dling. . . . . . . . . . . . . . . . . . . . . . . . . xxxvii xxxviii LIST OF TABLES I N T R O D U C T I O N "I have never been so wrong, in all my life."

Table 1 .

 1 1 -Basic Binary Floating-Point Formats Parameters.Because the storage of a floating-point number is only possible on a finite amount of memory, the data s, m, e must be stored on a finite number of bits. For the sign s, one bit is enough, but for the fields m and e, the IEEE 754-2008 standard specifies different lengths or ranges, which partially define floating-point formats. For binary

	Parameter	binary32 binary64 binary128
	k, storage width in bits	32	64	128
	p, precision in bits	24	53	113
	emax	127	1023	16383
	1.1.1 IEEE 754-2008 binary floating-point formats

floating-point arithmetic, i.e. when β = 2, a floating-point format is almost fully characterized by three integers [123, § 2.1]:

Table 1 .

 1 2 -Selection of IEEE 754-2008 Recommended Correctly Rounded Functions.

	Operation	Function	Domain
	exp	e x	[-∞, +∞]
	log log2 log10	log e (x) log 2 (x) log 10 (x)	[0, +∞]

  Table Maker's Dilemma (TMD). The Table Maker's Dilemma Most functions of Table 1.2 are transcendental, i.e. for most finite precision floating-point inputs, they return a transcendental number.Transcendental numbers do not have a repetitive sequence of digits in their fractional part, so they can be neither a finite precision floating-point number nor a midpoint between two finite precision floating-point numbers, but they may be extremely close to such values. As can be seen in Figure1.2, rounding values can change at these boundaries, which is why they are sometimes called rounding

	1.2.1

breakpoints

[123, p. 406

].

Table 1 .

 1 3 -Register Size and Availability of Intel's SIMD Extensions.

	Extension name Register size (bits)	ISA
	MMX	64	Pentium P55C
	SSE SSE2	128	Pentium III Pentium 4
	AVX AVX2	256	Sandy Bridge Haswell
	AVX512	512	Knights Landing
	four double precision (binary64) or eight single precision (binary32)
	floating-point numbers.

  table entry, this means that a naive exhaustive search would have to compute the LCM for every set made of 101 elements, each of which having 413 different possible values in average. Think of it as finding the right combination of a padlock that has 101 rotating discs, each of which having approximately -Number of PPTs per entry of a trigonometric table for p = 7 with a hypotenuse less than 2 18 . Each table entry contains at least one PPT.

	450 500							Mean
	400
	350
	300
	250
	200
	150
	100
	50
	0	0	10	20	30	40 Entry number 50 60	70	80	90 100
	Figure 2.7

Table 2 .

 2 1 -Exhaustive Search Results for sin and cos.

	p	k	n Time (s)	PPTs Hypotenuses
	3	425	9	1	87	33
	4	5525	13	1	1405	428
	5	160 225	18	0.2	42 329	11 765
	6 1 698 385	21	7	335 345	87 633
	7 6 569 225	23	31	1 347 954	335 645
	8	> 2 27 > 27	> 6700 > 21 407 993	> 4 976 110

Table 2 .

 2 2 -Exhaustive Search Results for sinh and cosh.

	p	k	n time (s)	PPTs Bigger legs
	3	144	8	1	23	12
	4	840	10	1	86	43
	5	10 080	14	1	1202	610
	6	171 360	18	9	18 674	9312
	7 1 081 080	21	328	147 748	72 476
	8	> 2 24 > 24 > 60 000 > 1 188 585	> 574 800

Table 2 .

 2 3 -A Trigonometric Exact Lookup Table Computed for p = 4.

	Index	S h	C h corr
	0	0	5525 +0x0.0000000000000p+0
	1 235 5520 -0x1.46e9e7603049fp-6
	2 612 5491 -0x1.cad996fe25a24p-7
	3 1036 5427 +0x1.27ac440de0a8cp-10
	4 1360 5355 -0x1.522b2a9e8491dp-10
	5 1547 5304 -0x1.d6513b89c7237p-6
	6 2044 5133 +0x1.038b12ae4eba1p-8
	7 2340 5005 -0x1.53f734851f48bp-13
	8 2600 4875 -0x1.49140da6fe454p-7
	9 2880 4715 -0x1.d02973d03a1f6p-7
	10 3315 4420 +0x1.2f1f464d3dc25p-6
	11 3500 4275 -0x1.7caa112f287aep-10
	12 3720 4085 -0x1.735972faced77p-7
	13 3952 3861 -0x1.fa6ed9240ab1ap-7

Table 2 .

 2 4 -A Hyperbolic Exact Lookup Table Computed for p = 5.

	Index	S h	C h	corr
	0	0	10 080 +0x0.0000000000000p+0
	1 284 10 084 -0x1.93963974f0cb6p-9
	2 651 10 101 +0x1.0b316b3c740d1p-9
	3 1064 10 136 +0x1.7c74108520aebp-7
	4 1190 10 150 -0x1.d8f891d50d1a1p-8
	5 1560 10 200 -0x1.13297ef8b55bbp-9
	6 1848 10 248 -0x1.535fdc36d3139p-8
	7 2222 10 322 -0x1.fe04ef1053a97p-15
	8 2560 10 400 +0x1.5891c9eaef76ap-10
	9 2940 10 500 +0x1.a58844d36e49ep-8
	10 3237 10 587 +0x1.b77a5031ebc86p-9
	11 3456 10 656 -0x1.dcf49bb32dc17p-8

table from table 2

 from2

	4000	Exact points Row boundaries ϕ = ln(2)/2
	3000
	a
	2000
	1000
	0	0	10080
		b

.3 in the Euclidean plane. Figure 2.10 -Visual representation of the hyperbolic exact lookup table from table 2.4 in the Euclidean plane.

Table 2 .

 2 5 -Prime Factorization of Found Common Multiples for sin and cos.

	k Prime factorization
	425	5 2 • 17
	5525	5 2 • 13 • 17
	160 225	5 2 • 13 • 17 • 29
	1 698 385	5 • 13 • 17 • 29 • 53
	6 569 225	5

2 • 13 • 17 • 29 • 41 Table 2.6 -Prime Factorization of Found Common Multiples for sinh and cosh.

Table 2 .

 2 7 -Heuristic Search Results for sin and cos.

	p	k n time (s)	PPTs Hypotenuses
	3	425	9	1	42	8
	4	5525 13	1	211	17
	5	160 225 18	1	996	40
	6	1 698 385 21	0.1	2172	66
	7	6 569 225 23	0.4	3453	69
	8	314 201 225 29	9.5 10 468	100
	9	12 882 250 225 34	294 20 312	109
	10 279 827 610 985 39	9393 33 057	110

Table 2 .

 2 8 -Heuristic Search Results for sinh and cosh.

	p	k n time (s) PPTs Bigger legs
	3	144	8	1	23	12
	4	840 10	1	65	24
	5	10 080 14	1	247	79
	6	180 180 18	1	917	193
	7	1 081 080 21	0.3 1743	248
	8	17 907 120 25	3.2 3909	388
	9	147 026 880 28	23 5802	400
	10 2 793 510 720 32	350 9012	502
	where P

π is the set of Pythagorean primes: P π = {5, 13, 17, 29, 37, 41, 53, 61, 73, . . . } .

  pp. 505-511]. Fässler use of several results from number theory: First, he brings back an antique way of representing any PPT (a, b, c), which dates back to at least the time of Euclid:

	Exhaustive search results Heuristic search results OLS regression (r ≈ 0.994) binary64 53-bit limit
	2 (k)
	log
	n =
	0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
	p
	Exhaustive search results Heuristic search results OLS regression (r ≈ 0.999) binary64 53-bit limit
	2 (k)
	log
	n =
	0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
	p

Figure 2.11 -Bit-length n of k with respect to the table index size p for the trigonometric functions.

Figure 2.12 -Bit-length n of k with respect to the table index size p for the hyperbolic functions.

makes

Table 2 .

 2 9 -"Fässler's" Heuristic Search Results for sin and cos.Description of the algorithmNow that we have improved our knowledge about PTs sharing a common hypotenuse and the lower bound for k, let us explain how we can use Equation (2.12) in a much faster heuristic to generate trigonometric tables. First, we must determine the decompositions of the first Pythagorean primes up to 73 as sums of two squares. This can be done by a precomputed simple exhaustive search: Assuming

	p		k	n time (s)	PPTs Hypotenuses
	3	425	9	1	43	8
	4	5525	13	1	225	16
	5	160 225	18	1	1071	34
	6	1 698 385	21	0.1	2317	49
	7	6 569 225	23	0.2	3757	50
	8	314 201 225	29	0.8 10 861	44
	9	12 882 250 225	34	2 19 243	27
	10	279 827 610 985	39	6 26 947	11
	11	3 929 086 318 625	42	6 29 434	3
	12 286 823 301 259 625	49	15 32 581	1
	13	not found > 53	1 32 576	0
	c > b > a > 0, we find
		c/2 < b	√ c -1 and 1 a	√ c -1 -1.

For example, if c = 73, we can search for (a, b) ∈ [1 . . 8] × [7 . . 8], which is almost instantaneous: 73 = 3 2 + 8 2 . Second, if the trigonometric table needs N triples, we initialize the lower bound for k with

Table 2 .

 2 10 -Costs of Addition and Multiplication of Expansions.

	Operation	Floating-point operations Without FMA With FMA

Table 2 .

 2 11 -Computation and Memory Costs for Three Table-Based Range Reductions for Trigonometric and Hyperbolic Functions. The number of memory accesses and the number of floating-point operations (without FMA in parentheses) are reported.

		Table Loads	Flops for the reconstruction (without FMA)
	Method	Quick Accurate	Trigonometric Quick Accurate	Hyperbolic Quick Accurate	Size (Bytes/row)
	Tang	4	6	34 (64) 113 (241) 102 (192)	339 (723)	48
	Gal	3	7	35 (63) 129 (257) 93 (179)	355 (739)	56
	Proposed	4	5	36 (64) 92 (148)	94 (180)	270 (510)

Table 2 .

 2 12 -Tang's Table for p = 4 With 16-Bit Precision.

	Index	S h	C h
	0 0x0.0000p+0 0x1.0000p+0
	1 0x1.ffaap-5 0x1.ff00p-1
	2 0x1.feaap-4 0x1.fc02p-1
	3 0x1.7dc2p-3 0x1.f706p-1
	4 0x1.faaep-3 0x1.f016p-1
	5 0x1.3ad2p-2 0x1.e734p-1
	6 0x1.7710p-2 0x1.dc6cp-1
	7 0x1.b1d8p-2 0x1.cfc6p-1
	8 0x1.eaeep-2 0x1.c152p-1
	9 0x1.110ep-1 0x1.b11ep-1
	10 0x1.2b92p-1 0x1.9f36p-1
	11 0x1.44ecp-1 0x1.8bb2p-1
	12 0x1.5d00p-1 0x1.76a0p-1
	13 0x1.73b8p-1 0x1.6018p-1

Table 3 .

 3 1 -Theoretical latencies for common schemes on a simplified Haswell architecture. -O2 optimization flag enabled, for an Intel Haswell CPU featuring the AVX2 and FMA instruction set extensions. These extensions can operate on 16 registers of 256 bits, either on a single word or on multiple words contained in the same register. The Haswell architecture is deeply pipelined, and has one floating-point adder and two multipliers/FMA, which allows for a throughput of four instructions per cycle[66, § 10]. Therefore it is able to leverage ILP on independent multiplications and independent additions/multiplications, although this is clearly bounded parallelism. GCC 5.2 is able to generate vector instructions via an autovectorization pass, which we used to produce vectorized binary code using scalar C source code.

	Scheme	Latency (cycles) without FMA with FMA
	Horner	56	35
	Horner-2	37	25
	Estrin	24	15
	(GCC) 5.2,

Table 4 .

 4 1 -Suitable Values for i (Satisfying Property 1) for Different Combinations of (e max , κ). • 2 -12 + 2 -24 + 1.5 • 2 -36[START_REF] Aad | Measurements of the Higgs boson production and decay rates and coupling strengths using pp collision data at √ s = 7 and 8 TeV in the ATLAS experiment[END_REF][START_REF] Berggren | Pytagoreiska trianglar [Pythagorean triangles[END_REF]

	e max	κ	i ∈ • • •
	127	1.5 • 2 -12	[1, 10]
	1023	2 -14	[1, 12]
	1023 1.5

Table 4 .

 4 2 -Implementation Parameters and Table and Polynomial Coefficient Sizes (# Bytes) for Various Values of i When p = 24.

	i d	θ 2	θ 3	table coefficient
	1 15 2 -52.69 2 -50.40	16	116
	2 11 2 -52.69 2 -50.27	32	84
	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

Table 4 .

 4 3 -Measured Performances of our Implementations Compared to the GNU Libmvec, in Cycles per Element. Performance of the log(x) functionIn this section, we measure the reciprocal throughput of the generated routines, and compare 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 tables indexed by i = 7 bits, and we use Metalibm and custom micro-benchmarks, for generated and handmade codes, respectively. Performance results are presented in Table4.3, in cycles per element (CPE).

	Version	ISA	fast-math	Workload	Rcpr. throughput Rcpr. throughput
					in binary32 (CPE) in binary64 (CPE)
	Libmvec 128 bits	SSE4	yes	PRNG	1.2	5.9
	Libmvec 256 bits	AVX2	yes	PRNG	0.3	1.3
	Libmvec 128 bits	SSE4	yes	Subnormals	6.7	40.5
	Libmvec 256 bits	AVX2	yes	Subnormals	2.6	15.2
	Handmade 128 bits	AVX2	no	Any	20.7	88.2
	Handmade 256 bits	AVX2	no	Any	11.4	44.9
	ML-generated 128 bits AVX2	no	Any	24.0	n/a
	ML-generated 256 bits AVX2	no	Any	16.0	n/a

Table 4 .

 4 4 -Measured Performances (CPE) on AVX2 of our Generated Implementations, in Binary32 Arithmetic and for Vectors of Size 4 (v4) or 8 (v8), With (sub) or Without (no-sub) Subnormal Handling. v4 / no-sub 58.4 49.5 24.5 22.2 22.0 20.6 21.3 v8 / no-sub 31.2 26.8 15.9 14.9 14.7 19.6 19.7

	i	3	4	5	6	7	8	9
	v4 / sub	60.7 51.0 46.0 24.0 24.0 21.6 22.5
	v8 / sub	32.3 27.5 16.8 15.7 16.0 14.1 14.0

See http://metalibm.org.

Latences minimales documentées pour l'architecture Haswell.

Voir http://cgpe.gforge.inria.fr/ et[START_REF] Mouilleron | Automatic Generation of Fast and Certified Code for Polynomial Evaluation[END_REF][START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors -Polynomial evaluation-based algorithms and certified code generation[END_REF].

Kalray, https://kalray.eu.

https://github.com/kalray/metalibm

CERN (Conseil Européen pour la Recherche Nucléaire [European Organization for Nuclear Research]) studies the most basic constituent of matter: the fundamental particles. See https://home.cern/about.

Sollya is a software tool used in this thesis and described in Chapter 1.

L'Agence ationale de la recherche, the French National Research Agency.

See http://metalibm.org.

The author humbly counts himself as part of such civilizations, with his two pentadactyl hands, although he also does his best to be in the set of people counting on their 10 hands, having 101 fingers each.

The elementary charge, usually noted e, is the electric charge carried by a proton. Its value is about 1.602 × 10 -19 C.

https://gcc.gnu.org/onlinedocs/gcc/Floating-Types.html

Not to be confused with Indian naans. These are delicious.

Except for some special cases of specific elementary functions, with which this thesis does not deal, fortunately.

http://www.netlib.org/cephes/

http://gruntthepeon.free.fr/ssemath/

http://software-lisc.fbk.eu/avx _ mathfun/

Commercial off-the-shelf.

"Full-custom application-specific integrated circuits" is a self-defined expression. They potentially maximize the hardware performance.

http://www.metalibm.org/

http://www.metalibm.org/ANRMetaLibm

https://github.com/kalray/metalibm

For instance, see the GNU libm in glibc

2.25.

OEIS Foundation Inc. (2018), The On-Line Encyclopedia of Integer Sequences, http://oeis.org/A002144.

See http://gnuplot.info/ and [177, p. 173].

See for instance the GNU libm implementations for generic architectures.

See https://software.intel.com/haswell.

Latencies for Haswell architectures[START_REF][END_REF].

See http://cgpe.gforge.inria.fr/.

https://kalray.eu

http://www.mpfr.org/ and[START_REF] Fousse | MPFR: A Multiple-precision Binary Floating-point Library with Correct Rounding[END_REF].

(cit. on p. 7).

Towards more efficient evaluation schemes

CGPE, 1 which stands for Code Generation for Polynomial Evaluation, is a software tool that can, among other features, compute low-latency polynomial evaluation schemes on unbounded parallelism for a given degree and architecture [START_REF] Mouilleron | Automatic Generation of Fast and Certified Code for Polynomial Evaluation[END_REF][START_REF] Revy | Implementation of binary floating-point arithmetic on embedded integer processors -Polynomial evaluation-based algorithms and certified code generation[END_REF]. Using our simplified SNPIP architecture defined in Section 3.2, we used CGPE to generate low-latency schemes of degree 12. This represents a total of 11 944 schemes, each having a latency of 26 cycles. The reciprocal throughputs of these schemes were measured for different inputs. Some of those inputs were chosen because they produce subnormal numbers when raised to a certain power. The schemes were compiled with FMA generation disabled and using GCC's autovectorizing capabilities. Figures 3.10, 3.11 and 3.12 show these 11 944 measurements as dots. Vertical lines split the set of polynomial evaluation schemes into subsets having the same number of multiplications. For example, evaluation schemes in the first subset (on the left) involve 17 multiplications, while in the last subset (on the right) they involve 25 multiplications.

In Figure 3.10, we note that all lowest-latency schemes have a greater throughput than Horner and Horner-2 schemes, but only about ten are at least as efficient as Estrin's scheme. These are the schemes that use the fewest multiplications. Quite rationally, thresholds indicating lower throughputs can be observed as the number of multiplications for a scheme increases, although they have the same theoretical latency on unbounded parallelism. Figure 3.11 confirms that Estrin's evaluation scheme can be a lot more efficient than most the power of meta-implementations u = r i • m -1. In addition the valuelog(2 τ • r i ) is retrieved from a lookup table using the i > 0 most significant fraction bits of the significand of r i .

Note that instead of extracting m and approximating 1/m , we could have preferred computing directly an approximation of 1/x . However when the result of 1/x falls into the subnormal range, that is, x > 2 e max -1 , the rcp instruction flushes the result to 0. To overcome this issue, a rescaling similar to the one used here for handling subnormal numbers should have been performed: but it would have been as costly as the extraction of m , or even more.

Note also 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 4.6.3.

4.2.2

How to determine r i , u, τ, andlog(2 τ • r i )?

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

Determination of r i

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

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

Hence:

where

Using the Intel Intrinsics Guide [START_REF][END_REF], 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 . On more recent micro-architectures though, the AVX-512 instruction set provides also a binary64 rcp instruction with ε rcp 2 -14 . Note that on SSE/AVX-supporting micro-architectures, the binary32 instruction can also be used to approximate 1/m for a binary64 input: m must first be cast to the binary32 format, thus slightly increasing the approximation error:

Hence from (4.8) and assuming ulp(log(x)) = ulp(r h ), we deduce that if

then (4.9) holds. When log

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. Now let us distinguish three cases:

1. When e + τ -λ = 0, and x ∈ X with It follows that if

then (4.10) holds, and r h is a faithful rounding of log(x). For example, for (p, i) = (24, 7), we have: min ≈ 3.5 • 10 -3 and ulp(min) = 2 -32 .

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 4.4.2.

guidelines to build accurate enough programs

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 constrain