Quels algorithmes d’exponentiation les langages de programmation / CPU utilisent-ils?

J’ai appris des algorithmes d’exponentiation plus rapides (k-aire, porte coulissante, etc.), et je me demandais lesquels sont utilisés dans les processeurs / langages de programmation? (Je ne sais pas si cela se produit dans le processeur ou via le compilateur)

Et juste pour les coups de pied, lequel est le plus rapide?

Modifier en ce qui concerne la largeur: C’est intentionnellement large parce que je sais qu’il existe un tas de techniques différentes pour faire cela. La réponse cochée avait ce que je cherchais.

Je suppose que votre intérêt est de mettre en œuvre les fonctions d’exponentiation que l’on trouve dans les bibliothèques de mathématiques standard pour HLL, en particulier C / C ++. Ceux-ci incluent les fonctions exp() , exp2() , exp10() et pow() , ainsi que les contreparties à simple précision expf() , exp2f() , exp10f() et powf() .

Les méthodes d’exponentiation que vous mentionnez (telles que k-aire, fenêtre glissante) sont généralement utilisées dans les algorithmes de chiffrement, tels que RSA, qui repose sur l’exponentiation. Ils ne sont généralement pas utilisés pour les fonctions d’exponentiation fournies par math.h ou cmath . Les détails d’implémentation des fonctions mathématiques standard telles que exp() diffèrent, mais un schéma commun suit un processus en trois étapes:

  1. réduction de l’argument de la fonction à un intervalle d’approximation primaire
  2. approximation d’une fonction de base appropriée sur l’intervalle d’approximation primaire
  3. mappage du résultat de l’intervalle primaire sur toute la plage de la fonction

Une étape auxiliaire est souvent le traitement de cas particuliers. Celles-ci peuvent concerner des situations mathématiques particulières telles que log(0.0) ou des opérandes spéciaux à virgule flottante tels que NaN (Pas un nombre).

Le code C99 pour expf(float) ci-dessous montre, à titre d’exemple, à quoi ces étapes ressemblent pour un exemple concret. L’argument a est d’abord divisé de telle sorte que exp(a) = e r * 2 i , où i est un entier et r est dans [log (sqrt (0.5), log (sqrt (2.0)]), l’intervalle d’approximation principal. Dans la deuxième étape, nous approchons maintenant avec un polynôme, qui peut être conçu en fonction de divers critères de conception, tels que la minimisation des erreurs absolues ou relatives, et évalué de différentes manières, notamment les schémas de Horner et d’Essortingn.

Le code ci-dessous utilise une approche très commune en utilisant une approximation minimax, qui minimise l’erreur maximale sur tout l’intervalle d’approximation. Un algorithme standard pour calculer de telles approximations est l’algorithme de Remez. L’évaluation se fait via le schéma de Horner; la précision numérique de cette évaluation est améliorée par l’utilisation de fmaf() .

Cette fonction mathématique standard implémente ce que l’on appelle un FLD multiply-add fusionné. Ceci calcule a*b+c utilisant le produit complet a*b lors de l’addition et en appliquant un arrondi unique à la fin. Sur la plupart des matériels modernes, tels que les GPU, les processeurs IBM Power, les processeurs x86 récents (Haswell, par exemple) et les processeurs ARM récents (en tant qu’extension facultative), cela correspond directement à une instruction matérielle. Sur les plates-formes dépourvues d’une telle instruction, fmaf() correspondra à un code d’émulation assez lent, auquel cas nous ne voudrions pas l’utiliser si la performance nous intéresse.

Le calcul final est la multiplication par 2 i , pour laquelle C et C ++ fournissent la fonction ldexp() . Dans le code de bibliothèque “à valeur indussortingelle”, on utilise généralement un idiome spécifique à une machine, qui tire parti de l’utilisation de l’arithmétique binary IEEE-754 pour float . Enfin, le code élimine les cas de débordement et de débordement.

La FPU x87 dans les processeurs x86 a une instruction F2XM1 qui calcule 2 x -1 sur [-1,1]. Ceci peut être utilisé pour la deuxième étape du calcul de exp() et exp2() . Il existe une instruction FSCALE qui est utilisée pour multiplier par i dans la troisième étape. Une façon courante d’implémenter F2XM1 est d’utiliser un microcode utilisant une approximation rationnelle ou polynomiale. Notez que la FPU x87 est maintenue principalement pour la prise en charge existante de nos jours. Sur la plate-forme x86 moderne, les bibliothèques utilisent généralement des implémentations logicielles pures basées sur SSE et des algorithmes similaires à celui présenté ci-dessous. Certains combinent des petites tables avec des approximations polynomiales.

pow(x,y) peut être implémenté de manière conceptuelle sous la forme exp(y*log(x)) , mais présente une perte de précision significative lorsque x est proche de l’unité et y de grande ampleur, ainsi qu’un traitement incorrect de nombreux cas spéciaux. spécifié dans les normes C / C ++. Une façon de contourner le problème de précision consiste à calculer log(x) et le produit y*log(x)) sous une forme de précision étendue. Les détails rempliraient une réponse complète, longue et séparée, et je n’ai pas de code à scope de main pour le démontrer. Dans diverses bibliothèques mathématiques C / C ++, pow(double,int) et powf(float, int) sont calculés par un chemin de code distinct qui applique la méthode “square-and-multiply” avec un balayage binary de la représentation binary de l’exposant entier.

 #include  /* import fmaf(), ldexpf() */ /* Like rintf(), but -0.0f -> +0.0f, and |a| must be < 2**22 */ float quick_and_dirty_rintf (float a) { float cvt_magic = 0x1.800000p+23f; return (a + cvt_magic) - cvt_magic; } /* Approximate exp(a) on the interval [log(sqrt(0.5)), log(sqrt(2.0))]. */ float expf_poly (float a) { float r; r = 0x1.694000p-10f; // 1.37805939e-3 r = fmaf (r, a, 0x1.125edcp-07f); // 8.37312452e-3 r = fmaf (r, a, 0x1.555b5ap-05f); // 4.16695364e-2 r = fmaf (r, a, 0x1.555450p-03f); // 1.66664720e-1 r = fmaf (r, a, 0x1.fffff6p-02f); // 4.99999851e-1 r = fmaf (r, a, 0x1.000000p+00f); // 1.00000000e+0 r = fmaf (r, a, 0x1.000000p+00f); // 1.00000000e+0 return r; } /* Approximate exp2() on interval [-0.5,+0.5] */ float exp2f_poly (float a) { float r; r = 0x1.418000p-13f; // 1.53303146e-4 r = fmaf (r, a, 0x1.5efa94p-10f); // 1.33887795e-3 r = fmaf (r, a, 0x1.3b2c6cp-07f); // 9.61833261e-3 r = fmaf (r, a, 0x1.c6af8ep-05f); // 5.55036329e-2 r = fmaf (r, a, 0x1.ebfbe0p-03f); // 2.40226507e-1 r = fmaf (r, a, 0x1.62e430p-01f); // 6.93147182e-1 r = fmaf (r, a, 0x1.000000p+00f); // 1.00000000e+0 return r; } /* Approximate exp10(a) on [log(sqrt(0.5))/log(10), log(sqrt(2.0))/log(10)] */ float exp10f_poly (float a) { float r; r = 0x1.a5a000p-3f; // 0.20587158 r = fmaf (r, a, 0x1.155dcap-1f); // 0.54173118 r = fmaf (r, a, 0x1.2bda68p+0f); // 1.17130136 r = fmaf (r, a, 0x1.046fa8p+1f); // 2.03465748 r = fmaf (r, a, 0x1.53524ap+1f); // 2.65094876 r = fmaf (r, a, 0x1.26bb1cp+1f); // 2.30258512 r = fmaf (r, a, 0x1.000000p+0f); // 1.00000000 return r; } /* Compute exponential base e. Maximum ulp error = 0.86565 */ float my_expf (float a) { float t, r; int i; t = a * 0x1.715476p+0f; // 1/log(2); 1.442695 t = quick_and_dirty_rintf (t); i = (int)t; r = fmaf (t, -0x1.62e400p-01f, a); // log_2_hi; -6.93145752e-1 r = fmaf (t, -0x1.7f7d1cp-20f, r); // log_2_lo; -1.42860677e-6 t = expf_poly (r); r = ldexpf (t, i); if (a < -105.0f) r = 0.0f; if (a > 105.0f) r = 1.0f/0.0f; // +INF return r; } /* Compute exponential base 2. Maximum ulp error = 0.86770 */ float my_exp2f (float a) { float t, r; int i; t = quick_and_dirty_rintf (a); i = (int)t; r = a - t; t = exp2f_poly (r); r = ldexpf (t, i); if (a < -152.0f) r = 0.0f; if (a > 152.0f) r = 1.0f/0.0f; // +INF return r; } /* Compute exponential base 10. Maximum ulp error = 0.95588 */ float my_exp10f (float a) { float r, t; int i; t = a * 0x1.a934f0p+1f; // log2(10); 3.321928 t = quick_and_dirty_rintf (t); i = (int)t; r = fmaf (t, -0x1.344136p-2f, a); // log10(2)_hi; 3.01030010e-1 r = fmaf (t, 0x1.ec10c0p-27f, r); // log10(2)_lo; 1.43209888e-8 t = exp10f_poly (r); r = ldexpf (t, i); if (a < -46.0f) r = 0.0f; if (a > 46.0f) r = 1.0f/0.0f; // +INF return r; }