Point flottant : Puis-je me passer de la partie collante si je ne fais qu'arrondir les égalités à l'unité ?

Question

J'écris des fonctions en virgule flottante pour un processeur 8 bits très limité. Lorsque je décale le significand vers la droite afin d'aligner les points binaires pour une addition ou une soustraction, je décale les bits qui tombent du côté droit du champ du significand dans un octet d'arrondi. Si je n'applique que la méthode d'arrondi "round ties to away", dois-je conserver un sticky bit ou prêter attention à autre chose que le bit le plus significatif de l'octet d'arrondi ? Je pense qu'avec cette méthode d'arrondi, si le MSB est 1, j'arrondis à l'écart de zéro (vers le haut pour les nombres positifs), sinon vers zéro.

Answer 1

Ce qui suit suppose l'utilisation de l'un des types binaires à virgule flottante de la norme IEEE-754, par exemple binary32 , binary64 .

Si l'on considère le significand juste avant l'arrondi, alors après une éventuelle renormalisation (lorsque le significand est en dehors de l'intervalle [1, 2] après l'addition ou la multiplication du significand) ou une dénormalisation (pour les résultats provisoirement inférieurs à la normale), l'arrondi nécessite trois éléments d'information lors de l'arrondi au plus proche, égalité des chances : La partie ronde R Il s'agit en fait d'un bit de signification supplémentaire dont le poids numérique est égal à la moitié du bit de signification stocké le moins significatif, le "bit" collant S qui est le OU implicite de tous les bits de poids significatifs numériquement inférieurs au bit rond, et du bit de poids significatif le plus faible du significand, L .

En dehors des cas de cravate, R décide seul si le significatif doit être incrémenté (lorsque R==1 ) ou non. Dans un cas d'égalité, où le résultat intermédiaire infiniment précis est exactement à mi-chemin entre les deux valeurs significatives représentables les plus proches, nous avons R==1 y S==0 et doivent examiner L pour satisfaire à la disposition relative à l'égalité des chances. En résumé, le significand doit être incrémenté lorsque (R==1) && ((S==1) || (L==1)) .

Comme indiqué dans la question, cette logique d'arrondi est simplifiée lorsque l'on utilise le mode d'arrondi au plus proche, liens à l'écart . Le significand doit être incrémenté lorsque (R==1) && ((S==X) || (L==X)) (où X signifie : sans importance), car il n'est pas nécessaire de distinguer les cas où les bits rejetés représentent une valeur supérieure à 0,5 ulp ou exactement égale à 0,5 ulp, et il n'est pas nécessaire d'examiner le significand pour voir s'il est pair. Nous nous retrouvons avec R==1 comme seule condition d'incrémentation du significand, et il n'est pas nécessaire d'accumuler des informations "bit" collantes.

J'ai utilisé le terme "bit collant" parce que dans les émulations logicielles, il est d'usage de consacrer un registre aux informations rondes et collantes, le bit le plus significatif du registre représentant le bit rond et tous les bits de rang inférieur collectivement représentant (par OU implicite) l'information collante. Cela permet d'éviter de condenser l'information collante en un seul bit à un endroit particulier. Lorsque l'on utilise un octet d'arrondi comme indiqué dans la question, l'incrémentation conditionnelle du significand est basée sur la comparaison de l'octet d'arrondi avec 0x80 ou en ajoutant 0x80 o 0x7F (en fonction du type d'arrondi le plus proche utilisé) et de propager la retenue.

Parce que l'arrondi au plus proche, l'égalisation met en œuvre impartial Je recommanderais fortement de mettre en œuvre le premier mode, malgré son coût légèrement plus élevé. Cela est également utile pour maintenir la compatibilité avec d'autres plates-formes.

Ci-dessous, je montre le code d'émulation pour IEEE-754 binary32 qui démontre que le maintien de l'information collante et la mise en œuvre des égalités ajoutent moins de 10 % de surcharge, ce qui devrait être le cas pour l'émulation sur un processeur orienté octet. Trois algorithmes différents pour l'arrondi au plus proche avec égalité des chances sont indiquées. L'approche la plus efficace dépend des spécificités de la plate-forme ; j'ai utilisé les trois variantes dans le passé. Je m'excuse pour les incohérences : Les deux morceaux de code ont été écrits à des moments différents, dans des contextes différents et dans des styles quelque peu différents. Notez que ce code utilise un mot d'arrondi 32 bits plutôt qu'un octet d'arrondi puisque le code a été écrit pour une plate-forme 32 bits.

Code pour binary32 En outre, avec le mot d'arrondi dans la variable rdnstk :

#define RND_VARIANT (3)  // 1 ... 3
#define CLZ_VARIANT (4)  // 1 ... 5

int clz (uint32_t a);

uint32_t fp32_add_rn_kernel (uint32_t a, uint32_t b)
{
    /* sort arguments my magnitude. Larger one becomes a, smaller becomes b */
    if ((b << 1) > (a << 1)) {
        uint32_t t = a;
        a = b;
        b = t;
    }

    uint32_t expo_a_m1 = ((a >> 23) & 0xff) - 1;
    uint32_t expo_b_m1 = ((b >> 23) & 0xff) - 1;

    if ((expo_a_m1 >= 0xfe) || (expo_b_m1 >= 0xfe)) {
        uint32_t am = a << 1;
        uint32_t bm = b << 1;
        if (bm == 0) {                // b is zero
            if (am == 0) a = a & b;   // a is also zero
            if (am > 0xff000000) a = a | 0x00400000; // a is NaN, quieten it
            return a;
        } else if ((expo_b_m1 > 0xfe) &&  // b is subnormal
                   (expo_a_m1 != 0xfe)) { // a is not NaN or INF
            int shift = clz (bm) - 8;
            bm = bm << shift;
            expo_b_m1 = ~shift;
            b = (b & 0x80000000) | bm;
            if (expo_a_m1 > 0xfe) { // a is subnormal
                int shift = clz (am) - 8;
                am = am << shift;
                expo_a_m1 = ~shift;
                a = (a & 0x80000000) | am;
            }
        } else if (am > 0xff000000) {     // a is a NaN
            return a | 0x00400000;
        } else if (bm > 0xff000000) {     // b is a NaN
            return b | 0x00400000;
        } else if ((am & bm) == 0xff000000) { // a is INF and b is INF
            return a | ((a ^ b) ? 0xffc00000 : 0);
        } else if ((am == 0) || (bm == 0xff000000)) { // a is zero or b is INF
            return b;
        } else if ((bm == 0) || (am == 0xff000000)) { // b is zero or a is INF
            return a;
        } 
    }
    uint32_t expo_diff = expo_a_m1 - expo_b_m1;
    uint32_t sign_a = a & 0x80000000;
    uint32_t r;
    if (expo_diff <= 25) {  
        uint32_t mant_a = (a & 0x00ffffff) | 0x00800000; // add in implict bit
        uint32_t mant_b = (b & 0x00ffffff) | 0x00800000; // add in implict bit
        uint32_t rndstk = expo_diff ? (mant_b << (32 - expo_diff)) : 0;
        if ((int32_t)(a ^ b) < 0) { // effective subtraction
            rndstk = 0 - rndstk;
            mant_a = mant_a - (mant_b >> expo_diff) - (rndstk != 0);
            if (!(mant_a & 0x00800000)) { /* cancellation of leading bits */
                if ((mant_a | rndstk) == 0) {
                    return mant_a; /* addends cancel completely */
                }
                /* renormalize */
                uint32_t shift = clz (mant_a) - 8;
                mant_a = (mant_a << shift) | (rndstk >> (32 - shift));
                rndstk <<= shift;
                expo_a_m1 -= shift;
            }
        } else {
            mant_a = mant_a + (mant_b >> expo_diff);
            if (mant_a & 0x01000000) { // significand overflow, renormalize
                rndstk = (rndstk >> 1) | (mant_a << 31);
                mant_a = mant_a >> 1;
                expo_a_m1 += 1;
            }
        }
        if (expo_a_m1 < 0xfe) {
            /* normal: round to nearest or even */
#if RND_VARIANT == 1
            r = (expo_a_m1 << 23) + mant_a + 
                ((rndstk | (mant_a & 1)) > 0x80000000);
#elif RND_VARIANT == 2
            r = (expo_a_m1 << 23) + mant_a + (rndstk >= 0x80000000);
            if (rndstk == 0x80000000) r = r & ~1; // tie case
#elif RND_VARIANT == 3
            r = (expo_a_m1 << 23) + mant_a + 
                ((rndstk == 0x80000000) ? (mant_a & 1) : (rndstk >> 31));
#else // RND_VARIANT
#error unsupported RND_VARIANT
#endif // RND_VARIANT 
        } else if (expo_a_m1 > 0xfe) { // underflow
            int denorm_shift = 0 - expo_a_m1;
            r = mant_a >> denorm_shift;
        } else { // overflow
            r = 0x7f800000;
        }
    } else {
        r = a;
    }
    r = sign_a | r;
    return r;
}

int clz (uint32_t a)
{
#if CLZ_VARIANT == 1
    static const uint8_t clz_tab[32] = {
        31, 22, 30, 21, 18, 10, 29,  2, 20, 17, 15, 13, 9,  6, 28, 1,
        23, 19, 11,  3, 16, 14,  7, 24, 12,  4,  8, 25, 5, 26, 27, 0
    };
    a |= a >> 16;
    a |= a >> 8;
    a |= a >> 4;
    a |= a >> 2;
    a |= a >> 1;
    return clz_tab [0x07c4acddu * a >> 27] + (!a);
#elif CLZ_VARIANT == 2
    uint32_t b;
    int n = 31 + (!a);
    if ((b = (a & 0xffff0000u))) { n -= 16;  a = b; }
    if ((b = (a & 0xff00ff00u))) { n -=  8;  a = b; }
    if ((b = (a & 0xf0f0f0f0u))) { n -=  4;  a = b; }
    if ((b = (a & 0xccccccccu))) { n -=  2;  a = b; }
    if ((    (a & 0xaaaaaaaau))) { n -=  1;         }
    return n;
#elif CLZ_VARIANT == 3
    int n = 0;
    if (!(a & 0xffff0000u)) { n |= 16;  a <<= 16; }
    if (!(a & 0xff000000u)) { n |=  8;  a <<=  8; }
    if (!(a & 0xf0000000u)) { n |=  4;  a <<=  4; }
    if (!(a & 0xc0000000u)) { n |=  2;  a <<=  2; }
    if ((int32_t)a >= 0) n++;
    if ((int32_t)a == 0) n++;
    return n;
#elif CLZ_VARIANT == 4
    uint32_t b;
    int n = 32;
    if ((b = (a >> 16))) { n = n - 16;  a = b; }
    if ((b = (a >>  8))) { n = n -  8;  a = b; }
    if ((b = (a >>  4))) { n = n -  4;  a = b; }
    if ((b = (a >>  2))) { n = n -  2;  a = b; }
    if ((b = (a >>  1))) return n - 2;
    return n - a;
#elif CLZ_VARIANT == 5
    return __builtin_clz (a);
#else
#error unknown CLZ_VARIANT
#endif
}

uint32_t float_as_uint32 (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

float uint32_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

float fp32_add_rn (float a, float b)
{
    uint32_t ai, bi;

    ai = float_as_uint32 (a);
    bi = float_as_uint32 (b);
    return uint32_as_float (fp32_add_rn_kernel (ai, bi));
}

Il s'agit de l'émulation de binary32 multiplication, avec le mot d'arrondi en variable mantr_lo :

#define RND_VARIANT (3)     // 1 ... 3

#define FLOAT_MANT_BITS    (23)
#define FLOAT_EXPO_BITS    (8)
#define FLOAT_EXPO_BIAS    (127)
#define FLOAT_MANT_MASK    (~((~0u) << (FLOAT_MANT_BITS+1))) /* incl. integer bit */
#define EXPO_ADJUST        (1)   /* adjustment for performance reasons */
#define MIN_NORM_EXPO      (1)   /* minimum biased exponent of normals */
#define MAX_NORM_EXPO      (254) /* maximum biased exponent of normals */
#define INF_EXPO           (255) /* biased exponent of infinities */
#define EXPO_MASK          (~((~0u) << FLOAT_EXPO_BITS))
#define FLOAT_SIGN_MASK    (0x80000000u)
#define FLOAT_IMPLICIT_BIT (1 << FLOAT_MANT_BITS)
#define RND_BIT_SHIFT      (31)
#define RND_BIT_MASK       (1 << RND_BIT_SHIFT)
#define FLOAT_INFINITY     (0x7f800000)
#define FLOAT_INDEFINITE   (0xffc00000u)
#define MANT_LSB           (0x00000001)
#define FLOAT_QNAN_BIT     (0x00400000)
#define MAX_SHIFT          (FLOAT_MANT_BITS + 2)

uint32_t fp32_mul_core (uint32_t a, uint32_t b)
{
    uint64_t prod;
    uint32_t expoa, expob, manta, mantb, shift;
    uint32_t r, signr, expor, mantr_hi, mantr_lo;

    /* split arguments into sign, exponent, significand */
    expoa = ((a >> FLOAT_MANT_BITS) & EXPO_MASK) - EXPO_ADJUST;
    expob = ((b >> FLOAT_MANT_BITS) & EXPO_MASK) - EXPO_ADJUST;
    manta = (a | FLOAT_IMPLICIT_BIT) & FLOAT_MANT_MASK;
    mantb = (b | FLOAT_IMPLICIT_BIT) & FLOAT_MANT_MASK;
    /* result sign bit: XOR sign argument signs */
    signr = (a ^ b) & FLOAT_SIGN_MASK;
    if ((expoa >= (MAX_NORM_EXPO - EXPO_ADJUST)) || /* at least one argument is special */
        (expob >= (MAX_NORM_EXPO - EXPO_ADJUST))) { 
        if ((a & ~FLOAT_SIGN_MASK) > FLOAT_INFINITY) { /* a is NaN */
            /* return quietened NaN */
            return a | FLOAT_QNAN_BIT;
        }
        if ((b & ~FLOAT_SIGN_MASK) > FLOAT_INFINITY) { /* b is NaN */
            /* return quietened NaN */
            return b | FLOAT_QNAN_BIT;
        }
        if ((a & ~FLOAT_SIGN_MASK) == 0) { /* a is zero */
            /* return NaN if b is infinity, else zero */
            return (expob != (INF_EXPO - EXPO_ADJUST)) ? signr : FLOAT_INDEFINITE;
        }
        if ((b & ~FLOAT_SIGN_MASK) == 0) { /* b is zero */
            /* return NaN if a is infinity, else zero */
            return (expoa != (INF_EXPO - EXPO_ADJUST)) ? signr : FLOAT_INDEFINITE;
        }
        if (((a & ~FLOAT_SIGN_MASK) == FLOAT_INFINITY) || /* a or b infinity */
            ((b & ~FLOAT_SIGN_MASK) == FLOAT_INFINITY)) {
            return signr | FLOAT_INFINITY;
        }
        if ((int32_t)expoa < (MIN_NORM_EXPO - EXPO_ADJUST)) { /* a is subnormal */
            /* normalize significand of a */
            manta = a & FLOAT_MANT_MASK;
            expoa++;
            do {
                manta = 2 * manta;
                expoa--;
            } while (manta < FLOAT_IMPLICIT_BIT);
        } else if ((int32_t)expob < (MIN_NORM_EXPO - EXPO_ADJUST)) { /* b is subnormal */
            /* normalize significand of b */
            mantb = b & FLOAT_MANT_MASK;
            expob++;
            do {
                mantb = 2 * mantb;
                expob--;
            } while (mantb < FLOAT_IMPLICIT_BIT);
        }
    }
    /* result exponent: add argument exponents and adjust for biasing */
    expor = expoa + expob - FLOAT_EXPO_BIAS + 2 * EXPO_ADJUST;
    mantb = mantb << FLOAT_EXPO_BITS; /* preshift to align result signficand */
    /* result significand: multiply argument significands */
    prod = (uint64_t)manta * mantb;
    mantr_hi = (uint32_t)(prod >> 32);
    mantr_lo = (uint32_t)(prod >>  0);
    /* normalize significand */
    if (mantr_hi < FLOAT_IMPLICIT_BIT) {
        mantr_hi = (mantr_hi << 1) | (mantr_lo >> (32 - 1));
        mantr_lo = (mantr_lo << 1);
        expor--;
    }
    if (expor <= (MAX_NORM_EXPO - EXPO_ADJUST)) { /* normal, may overflow to infinity during rounding */
        /* combine biased exponent, sign and signficand */
        r = (expor << FLOAT_MANT_BITS) + signr + mantr_hi;
        /* round result to nearest or even; overflow to infinity possible */
#if RND_VARIANT == 1
        r = r + ((mantr_lo | (mantr_hi & MANT_LSB)) > RND_BIT_MASK);
#elif RND_VARIANT == 2
        r = r + (mantr_lo >= RND_BIT_MASK);
        if (mantr_lo == 0x80000000) r = r & ~MANT_LSB; // tie case
#elif RND_VARIANT == 3
        r = r + ((mantr_lo == RND_BIT_MASK) ? (mantr_hi & MANT_LSB) : (mantr_lo >> RND_BIT_SHIFT));
#else // RND_VARIANT
#error unsupported RND_VARIANT
#endif // RND_VARIANT 
    } else if ((int32_t)expor > (MAX_NORM_EXPO - EXPO_ADJUST)) { /* overflow */
        /* return infinity */
        r = signr | FLOAT_INFINITY;
    } else { /* underflow */
        /* return zero, normal, or smallest subnormal */
        shift = 0 - expor;
        if (shift > MAX_SHIFT) shift = MAX_SHIFT;
        /* denormalize significand */
        mantr_lo = mantr_hi << (32 - shift) | (mantr_lo ? 1 : 0);
        mantr_hi = mantr_hi >> shift;
        /* combine sign and signficand; biased exponent known to be zero */
        r = mantr_hi + signr;
        /* round result to nearest or even */
#if RND_VARIANT == 1
        r = r + ((mantr_lo | (mantr_hi & MANT_LSB)) > RND_BIT_MASK);
#elif RND_VARIANT == 2
        r = r + (mantr_lo >= RND_BIT_MASK);
        if (mantr_lo == 0x80000000) r = r & ~MANT_LSB; // tie case
#elif RND_VARIANT == 3
        r = r + ((mantr_lo == RND_BIT_MASK) ? (mantr_hi & MANT_LSB) : (mantr_lo >> RND_BIT_SHIFT));
#else // RND_VARIANT
#error unsupported RND_VARIANT
#endif // RND_VARIANT 
     }
    return r;
}

uint32_t float_as_uint (float a)
{
    uint32_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

float uint_as_float (uint32_t a)
{
    float r;
    memcpy (&r, &a, sizeof r);
    return r;
}

float fp32_mul (float a, float b)
{
    return uint_as_float (fp32_mul_core (float_as_uint (a), float_as_uint (b)));
}