Division efficace en virgule flottante avec diviseurs de nombres entiers constants

Question

Une récente question, si les compilateurs sont autorisés à remplacer division flottante avec multiplication en virgule flottante, m'a poussé à poser cette question.

En vertu de l'exigence stricte, que les résultats obtenus après transformation du code est binaire identique à l'actuelle opération de division, il est trivial de voir que pour les binaires IEEE-754 arithmétique, il est possible pour les diviseurs d'une puissance de deux. Tant que la réciproque de le diviseur est représentable, en multipliant par l'inverse du diviseur donne des résultats identiques à la division. Par exemple, la multiplication par 0.5 peut remplacer la division par 2.0.

On se demande alors pour quelle autre diviseurs de tels remplacements de travail, en supposant que nous permettent de tout de courtes séquences d'instructions qui remplace la division, mais est significativement plus rapide, tout en offrant des bits des résultats identiques. En particulier permettre fused multiply-add opérations en plus de la simple multiplication. Dans les commentaires je l'ai fait à la suite des papier:

Nicolas Brisebarre, Jean-Michel Muller, et Saurabh Kumar Raina. L'accélération arrondie division flottante lorsque le diviseur est connu à l'avance. IEEE Transactions sur les Ordinateurs, Vol. 53, n ° 8, août 2004, pp. 1069-1072.

La technique préconisée par les auteurs de l'article précalcule la réciproque du diviseur de y comme un normalisée tête-queue paire z_h:z_l comme suit: z_h = 1 / y, z_l = fma (y, z,_h, 1) / y. Plus tard, la division de q = x / y est calculée comme q = fma (z,_h, x, z_l * x). Le papier provient de diverses conditions qui diviseur de y doit satisfaire pour que cet algorithme fonctionne. Que l'on observe facilement, cet algorithme a des problèmes avec les infinis et zéro lorsque les premiers signes de la tête et la queue diffèrent. Plus important encore, il ne parvient pas à fournir des résultats corrects pour les dividendes x qui sont de très faible ampleur, car le calcul du quotient de la queue, z_l * x, souffre de dépassement de capacité.

Le document fait également le passage d'une référence à un autre FMA-division basée sur l'algorithme, mis au point par Peter Markstein, quand il était à IBM. La référence pertinente est:

P. W. Markstein. Calcul de fonctions élémentaires sur les IBM RISC System/6000 processeur. IBM Journal de Recherche & Développement, Vol. 34, N ° 1, janvier 1990, pp. 111-119

Dans Markstein de l'algorithme, on calcule une réciproque rc, à partir de laquelle un premier quotient q = x * rc est formé. Ensuite, le reste de la division est calculée avec précision avec un FMA comme r = fma (-y, q, x), et une amélioration de plus justes, plus le quotient est enfin calculée comme q = fma (r, rc, q).

Cet algorithme a aussi des problèmes de x qui sont des zéros ou des infinis (facilement avec approprié de l'exécution conditionnelle), mais des tests exhaustifs à l'aide de la norme IEEE-754 simple précision float des données montre qu'il fournit la bonne quotient dans tous les possibe dividendes x pour de nombreux diviseurs y, parmi ces nombreux petits nombres. Ce code C la met en œuvre:

/* precompute reciprocal */
rc = 1.0f / y;

/* compute quotient q=x/y */
q = x * rc;
if ((x != 0) && (!isinf(x))) {
    r = fmaf (-y, q, x);
    q = fmaf (r, rc, q);
}

Sur la plupart des architectures de processeur, cela devrait se traduire par une dépourvu de branches séquence d'instructions, en utilisant soit la prédication, à condition de déplacements, ou sélectionnez-le type d'instructions. Pour donner un exemple concret: Pour la division par 3.0f, l' nvcc compilateur CUDA 7.5 génère les éléments suivants du code machine pour un Kepler GPU:

    LDG.E R5, [R2];                        // load x
    FSETP.NEU.AND P0, PT, |R5|, +INF , PT; // pred0 = fabsf(x) != INF
    FMUL32I R2, R5, 0.3333333432674408;    // q = x * (1.0f/3.0f)
    FSETP.NEU.AND P0, PT, R5, RZ, P0;      // pred0 = (x != 0.0f) && (fabsf(x) != INF)
    FMA R5, R2, -3, R5;                    // r = fmaf (q, -3.0f, x);
    MOV R4, R2                             // q
@P0 FFMA R4, R5, c[0x2][0x0], R2;          // if (pred0) q = fmaf (r, (1.0f/3.0f), q)
    ST.E [R6], R4;                         // store q

Pour mes expériences, j'ai écrit le petit C test de programme ci-dessous les étapes à travers entier diviseurs dans l'ordre croissant et pour chacun d'eux de manière exhaustive les tests ci-dessus séquence de code à l'encontre de la bonne division. Il imprime une liste des diviseurs que passé ce test exhaustif. Sortie partielle se présente comme suit:

PASS: 1, 2, 3, 4, 5, 7, 8, 9, 11, 13, 15, 16, 17, 19, 21, 23, 25, 27, 29, 31, 32, 33, 35, 37, 39, 41, 43, 45, 47, 49, 51, 53, 55, 57, 59, 61, 63, 64, 65, 67, 69,

Pour intégrer l'algorithme de remplacement dans un compilateur comme une optimisation, une liste de diviseurs à laquelle le code ci-dessus de transformation peuvent être utilisées sans danger n'est pas pratique. La sortie du programme jusqu'à présent (à un taux d'environ un résultat par minute) suggère que le jeûne code fonctionne correctement dans tous les encodages possibles de l' x pour ceux diviseurs y qui sont des entiers impairs ou sont des puissances de deux. Les preuves anecdotiques, pas une preuve, bien sûr.

Quel ensemble de conditions mathématiques peut déterminer a priori si la transformation de la division dans le code ci-dessus, la séquence est sûr? Réponses peut supposer que toutes les opérations à virgule flottante sont effectuées dans le mode d'arrondi par défaut de "tour la plus proche, ou encore".

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

int main (void)
{
    float r, q, x, y, rc;
    volatile union {
        float f;
        unsigned int i;
    } arg, res, ref;
    int err;

    y = 1.0f;
    printf ("PASS: ");
    while (1) {
        /* precompute reciprocal */
        rc = 1.0f / y;

        arg.i = 0x80000000;
        err = 0;
        do {
            /* do the division, fast */
            x = arg.f;
            q = x * rc;
            if ((x != 0) && (!isinf(x))) {
                r = fmaf (-y, q, x);
                q = fmaf (r, rc, q);
            }
            res.f = q;
            /* compute the reference, slowly */
            ref.f = x / y;

            if (res.i != ref.i) {
                err = 1;
                break;
            }
            arg.i--;
        } while (arg.i != 0x80000000);

        if (!err) printf ("%g, ", y);
        y += 1.0f;
    }
    return EXIT_SUCCESS;
}

Answer 1

Cette question demande un moyen d'identifier les valeurs de la constante Y qui font qu'il est sécuritaire de le transformer x / Y en un moins cher calcul à l'aide de FMA pour toutes les valeurs possibles de l' x. Une autre approche consiste à utiliser l'analyse statique pour déterminer une sur-approximation des valeurs de x peut prendre, de sorte que le généralement dénuée de transformation peut être appliquée dans la connaissance que les valeurs pour lesquelles la transformation du code diffère de l'original de la division ne se produisent pas.

À l'aide de représentations des ensembles de valeurs en virgule flottante qui sont bien adaptées aux problèmes de l'flottantes, même avant l'analyse à partir du début de la fonction peut produire de l'information utile. Par exemple:

float f(float z) {
  float x = 1.0f + z;
  float r = x / Y;
  return r;
}

En supposant que le défaut d'arrondi au plus proche de mode(*), dans la fonction ci-dessus x ne peut être NaN (si l'entrée est NaN), +0.0 f, ou un nombre plus grand que 2^-24 en ampleur, mais pas -0.0 f ou quelque chose de plus proche de zéro que 2^-24. Ce qui justifie la transformation en l'un des deux formes de la question pour de nombreuses valeurs de la constante Y.

(*) hypothèse sans laquelle de nombreuses optimisations sont impossibles et que les compilateurs C déjà moins que le programme utilise explicitement #pragma STDC FENV_ACCESS ON

---

Avant l'analyse statique, qui prédit l'information pour x ci-dessus peut être basée sur une représentation des ensembles de valeurs à virgule flottante, une expression peut prendre comme un n-uplet de:

• une représentation pour les ensembles de possibles valeurs NaN (Depuis comportements de NaN sont underspecified, un choix est à utiliser uniquement une valeur booléenne, true sens certains NaNs peut être présent, et false indiquant l'absence de NaN est présent.),
• quatre boolean drapeaux qui indiquent respectivement la présence de +inf, -inf, +0.0, -0.0,
• un intervalle inclusif de négatif finis les valeurs à virgule flottante, et
• un intervalle inclusif de positif finis les valeurs à virgule flottante.

Afin de suivre cette approche, toutes les opérations à virgule flottante qui peut se produire dans un programme C doit être compris par l'analyseur statique. Pour illustrer, l'addition entre des ensembles de valeurs de U et de V, afin d'être utilisée pour manipuler + dans le code analysé, peut être mis en œuvre:

• Si NaN est présent dans l'un des opérandes, ou si les opérandes peuvent être infinis de signes opposés, NaN est présent dans le résultat.
• Si la valeur 0 ne peut pas être un résultat de l'addition d'une valeur de U et une valeur de V, utilisation standard de l'arithmétique d'intervalle. La limite supérieure de la résultat est obtenu pour l'arrondi au plus proche plus de la plus grande valeur de U et la plus grande valeur de V, de sorte que ces limites doivent être calculées avec arrondi au plus proche.
• Si la valeur 0 peut être un résultat de l'addition d'une valeur positive de U et une valeur négative de V, alors soit M la plus petite valeur positive dans U telle que M est présent dans V.
  • si le succ(M) est présent dans U, alors cette paire de valeurs contribue succ(M) - M pour les valeurs positives de la suite.
  • si -succ(M) est présent dans V, alors cette paire de valeurs contribue à la valeur négative de M - succ(M) pour les valeurs négatives de la suite.
  • si pred(M) est présent dans U, alors cette paire de valeurs contribue à la valeur négative pred(M) - M pour les valeurs négatives de la suite.
  • si -pred(M) est présent dans V, alors cette paire de valeurs contribue à la valeur de M - pred(M) pour les valeurs positives de la suite.
• Faire le même travail si la valeur 0 peut être le résultat de l'addition de la valeur négative de U et une valeur positive de V.

Remerciements: le ci-dessus emprunte des idées à partir de "l'Amélioration de la virgule Flottante d'Addition et de Soustraction des Contraintes", Bruno Marre & Claude Michel

---

Exemple: la compilation de la fonction f ci-dessous:

float f(float z, float t) {
  float x = 1.0f + z;
  if (x + t == 0.0f) {
    float r = x / 6.0f;
    return r;
  }
  return 0.0f;
}

L'approche en question refuse de transformer la division dans la fonction f dans une autre forme, parce que 6 n'est pas celui de la valeur pour laquelle la division peut être inconditionnellement transformé. Au lieu de cela, ce que je suggère, c'est l'application d'une simple analyse de la valeur à partir du début de la fonction qui, dans ce cas, détermine qu' x est un corps fini de flotter +0.0f ou d'au moins 2^-24 en ampleur, et à utiliser cette information à appliquer Brisebarre et al transformation du, confiant dans la connaissance que l' x * C2 n'a pas de dépassement de capacité.

Pour être explicite, je suggère d'utiliser un algorithme tel que celui ci-dessous pour décider si ou de ne pas transformer la division en quelque chose de plus simple:

1. Est - Y l'une des valeurs qui peuvent être transformés à l'aide Brisebarre et al de la méthode en fonction de leur algorithme?
2. N'C1 et C2 de la méthode qu'ils ont le même signe, ou est-il possible d'exclure la possibilité que le dividende est infini?
3. N'C1 et C2 de la méthode qu'ils ont le même signe, ou peut - x ne prendre que l'une des deux représentations de 0? Si dans le cas où C1 et C2 ont des signes différents et x ne peut être qu'une représentation de zéro, n'oubliez pas de violon(**) avec les signes de la FMA-base de calcul à faire produire le bon de zéro lors de l' x est égal à zéro.
4. Peut l'ampleur du dividende être garanti d'être assez grand pour exclure la possibilité qu' x * C2 underflows?

Si la réponse à ces quatre questions est "oui", puis la division peut être transformé en une multiplication et une FMA dans le contexte de la fonction en cours d'élaboration. L'analyse statique décrit ci-dessus sert à répondre à des questions 2., 3. et 4.

(**) "jouer avec les signes" signifie l'utilisation de -FMA(-C1, x, (-C2)*x) à la place de FMA(C1, x, C2*x) lorsque cela est nécessaire pour rendre le résultat correctement lorsque x ne peut être l'un des deux signé des zéros

Answer 2

Permettez-moi de redémarrer pour la troisième fois. Nous essayons d'accélérer

    q = x / y

où y est une constante entière, et q, x, et y sont tous IEEE 754-2008 binary32 des valeurs à virgule flottante. Ci-dessous, fmaf(a,b,c) indique une fused multiply-add a * b + c à l'aide de binary32 valeurs.

L'algorithme naïf est par l'intermédiaire d'un précalculées réciproque,

    C = 1.0f / y

de sorte que lors de l'exécution d'une (beaucoup plus rapide) multiplication suffit:

    q = x * C

Le Brisebarre-Muller-Raina accélération utilise deux précalculées constantes,

    zh = 1.0f / y
    zl = -fmaf(zh, y, -1.0f) / y

de sorte qu'au moment de l'exécution, une multiplication et une fused multiply-add suffit:

    q = fmaf(x, zh, x * zl)

Le Markstein algorithme combine l'approche naïve avec deux fused multiply-ajoute que donne le résultat correct si l'approche naïve donne un résultat dans un délai de 1 unité dans le moins significatif de la place, par précalculer

    C1 = 1.0f / y
    C2 = -y

de sorte que la division peut être approchée à l'aide de

    t1 = x * C1
    t2 = fmaf(C1, t1, x)
    q  = fmaf(C2, t2, t1)

---

L'approche naïve fonctionne pour toutes les puissances de deux y, mais sinon c'est assez mauvais. Par exemple, pour les diviseurs, les 7, 14, 15, 28, et 30, il donne un résultat incorrect pour plus de la moitié de tous les possibles x.

Le Brisebarre-Muller-Raina approche de la même façon échoue pour presque tous les non-puissance de deux y, mais beaucoup moins x de rendement le résultat incorrect (moins de la moitié de un pour cent de tous les possibles x, varie en fonction des y).

Le Brisebarre-Muller-Raina l'article montre que l'erreur maximale dans l'approche naïve est de ±1,5 ULPs.

Le Markstein approche donne des résultats corrects pour des puissances de deux y, et aussi pour entier impair y. (Je n'ai pas trouvé un défaut entier impair diviseur pour le Markstein approche.)

---

Pour le Markstein approche, j'ai analysé les diviseurs 1 - 19700 (données brutes ici).

Traçant le nombre de cas d'insuffisance (diviseur dans l'axe horizontal, le nombre de valeurs de x où Markstein approche échoue pour dit diviseur), on peut voir un modèle simple de se produire:

Notez que ces parcelles ont à la fois horizontale et verticale des axes logarithmiques. Il n'y a pas de points pour les impair de diviseurs, comme l'approche donne des résultats corrects pour tous les impair de diviseurs que j'ai testé.

Si nous changeons l'axe des x pour le bit inverse (chiffres binaires dans l'ordre inverse, c'est à dire 0b11101101 → 0b10110111, données) des diviseurs, nous avons un schéma clair:

Si nous traçons une ligne droite passant par le centre des ensembles de points, nous obtenons la courbe 4194304/x. (Rappelez-vous, l'intrigue ne tient compte que de la moitié de la possible flotte, de sorte que lorsque l'on considère tous les possible de flotteurs, de le doubler.) 8388608/x et 2097152/x support de l'ensemble du motif d'erreur complètement.

Ainsi, si nous utilisons rev(y) pour calculer le bit inverse du diviseur y, alors 8388608/rev(y) est une bonne approximation du premier ordre de nombre de cas (de tous les possibles float) où le Markstein approche donne un résultat incorrect pour une même, la non-puissance de deux diviseur y. (Ou, 16777216/rev(x) pour la limite supérieure.)

Ajouté 2016-02-28: j'ai trouvé un rapprochement pour le nombre de cas d'erreur à l'aide de la Markstein approche, étant donné un entier (binary32) diviseur. Ici, c'est que les pseudo-code:

function markstein_failure_estimate(divisor):
    if (divisor is zero)
        return no estimate
    if (divisor is not an integer)
        return no estimate

    if (divisor is negative)
        negate divisor

    # Consider, for avoiding underflow cases,
    if (divisor is very large, say 1e+30 or larger)
        return no estimate - do as division

    while (divisor > 16777216)
        divisor = divisor / 2

    if (divisor is a power of two)
        return 0

    if (divisor is odd)
        return 0

    while (divisor is not odd)
        divisor = divisor / 2

    # Use return (1 + 83833608 / divisor) / 2
    # if only nonnegative finite float divisors are counted!
    return 1 + 8388608 / divisor

Cela permet de corriger l'erreur d'estimation à l'intérieur de ±1 sur le Markstein cas j'ai testé (mais je n'ai pas encore suffisamment testé diviseurs de plus de 8388608). La finale de la division doit être telle qu'il n'y a pas de faux zéros, mais je ne peux pas garantir qu'il (encore). Il ne prend pas en compte les très grands diviseurs (dire 0x1p100, ou 1e+30, et de plus grande ampleur), qui ont des problèmes de dépassement de capacité, je n'irais certainement exclure de tels diviseurs de l'accélération de toute façon.

Dans les tests préliminaires, l'estimation semble étrangement précis. Je n'ai pas de dessiner un tracé de comparer les estimations et les erreurs réelles pour les diviseurs de 1 à 20000, parce que les points coïncident exactement dans les parcelles. (Dans cette gamme, l'estimation est exacte, ou un trop grand.) Essentiellement, les estimations de reproduire la première parcelle dans cette réponse, exactement.

---

Le modèle d'échecs pour le Markstein approche est régulière, et très intéressant. L'approche fonctionne pour toutes les puissances de deux diviseurs, et tout entier impair diviseurs.

Pour les diviseurs de plus de 16777216, j'ai toujours voir les mêmes erreurs que pour un diviseur qui est divisée par la plus petite puissance de deux pour donner une valeur de moins de 16777216. Par exemple, 0x1.3cdfa4p+23 et 0x1.3cdfa4p+41, 0x1.d8874p+23 et 0x1.d8874p+32, 0x1.cf84f8p+23 et 0x1.cf84f8p+34, 0x1.e4a7fp+23 et 0x1.e4a7fp+37. (Au sein de chaque paire, la mantisse est la même, et que la puissance de deux varie.)

En supposant que mon banc d'essai n'est pas dans l'erreur, cela signifie que le Markstein approche fonctionne aussi diviseurs de plus de 16777216 ampleur (mais plus petit que, disons, la 1e+30), si le diviseur est telle que lorsqu'il est divisé par la plus petite puissance de deux qui donne un quotient de moins de 16777216 en ampleur, et le quotient est impair.

Answer 3

J'aime @Pascal's réponse, mais en optimisation, il est souvent préférable d'avoir un moyen simple et bien compris le sous-ensemble de transformations plutôt que d'une solution parfaite.

Tous les actuels et historique commun à virgule flottante formats ont une chose en commun: un binaire de la mantisse.

Par conséquent, toutes les fractions ont été les nombres rationnels de la forme:

x / 2ⁿ

Ceci est en contraste à l'une des constantes dans le programme (et de tous les possibles en base 10 fractions) qui sont des nombres rationnels de la forme:

x / (2ⁿ * 5^m)

Ainsi, une optimisation serait tout simplement le test de l'entrée et réciproque pour m == 0, puisque ces nombres sont représentés exactement dans la FP de format et opérations avec eux devrait produire des données exactes dans le format.

Ainsi, par exemple, au sein de l' (virgule 2 chiffres) plage de .01 de 0.99 de diviser ou de multiplier par le nombre serait optimisé:

.25 .50 .75

Et tout le reste ne serait pas. (Je pense, ne tester d'abord, lol.)

Answer 4

Le résultat d'une division à virgule flottante est:

• un signe drapeau
• un significande
• un exposant
• un ensemble d'indicateurs (dépassement de, dépassement de capacité, inexactes, etc - voir fenv())

L'obtention de la première 3 pièces correctes (mais le jeu de drapeaux incorrect) n'est pas suffisant. Sans plus de connaissances (p. ex. pièces de pièces de la suite fait de l'importance, les valeurs possibles du dividende, etc) je suppose que le remplacement de la division par une constante avec la multiplication par une constante (et/ou un alambiqué FMA mess) est presque jamais à l'abri.

En outre, pour les Processeurs modernes, je n'ai pas supposer que le remplacement d'une division avec 2 Scfg est toujours une amélioration. Par exemple, si le goulot d'étranglement est l'instruction fetch/décoder, puis cette "optimisation" ferait rendement pire. Pour un autre exemple, si des instructions ne dépend pas du résultat (le CPU peut faire beaucoup d'autres instructions en parallèle pendant l'attente des résultats) FMA version peut présenter de multiples dépendance des stands et l'exécution est de pire en pire. Pour un troisième exemple, si tous les registres sont utilisés alors que le FMA version (ce qui exige plus de "live" variables) peut augmenter le "renversement" et de faire une performance pire.

Il faut noter que, dans de nombreux cas mais pas tous) de la division ou la multiplication par une constante multiple de 2 peut être fait avec plus seul (en particulier, l'ajout d'un décalage de l'exposant).