Comment mettre en œuvre la multiplication en virgule fixe

Question

Je veux implémenter la multiplication en C et en virgule fixe avec une largeur de bit de fraction évolutive (c'est à dire de 1 bit à 30 bits), je sais que la façon la plus simple est la suivante :

typedef int32_t fixedpt;
typedef int64_t fixedptd;
fixedpt fixedpt_mul(fixedpt A, fixedpt B)
{
    return (((fixedptd)A * (fixedptd)B) >> FIXEDPT_FBITS);
}

mais supposons que je ne puisse pas utiliser int64_t, cette méthode est donc susceptible de déborder lorsque la largeur de bit fractionnaire est importante. J'ai trouvé un repo existant : libfixmath qui sépare les parties entières et les parties fractionnaires avant la multiplication, cette méthode permet de résoudre le problème ci-dessus, mais elle n'implémente que la largeur de bit de 16 fractions, c'est pourquoi je la modifie pour qu'elle s'adapte à des cas plus généraux :

#include <stdint.h>
#include <stdio.h>

typedef int32_t fix_t;
// w means fraction bitwidth
#define fix_ONE(w)      ((fix_t)(0x00000001 << w))
#define MASK_M(w)       (((fix_t)1 << w) - 1) // mask of mantissa, all LSB set, all MSB clear
#define fix_rconst(x, w) ((fix_t)((x) * fix_ONE(w) + ((x) >= 0 ? 0.5 : -0.5)))  // convert a const to rounded fix_t
static const fix_t fix_overflow = 0x80000000; // 1000 0000 0000 0000 0000 0000 0000 0000
static inline float fix_to_float(fix_t a, int w) { return ((float)a / fix_ONE(w)); }

fix_t fix_mul(fix_t inArg0, fix_t inArg1, int w) {
    // w is fractional bitwidth
    // separate interger part and fraciton part
    int32_t A = (inArg0 >> w), C = (inArg1 >> w);
    uint32_t B = (inArg0 & MASK_M(w)), D = (inArg1 & MASK_M(w));

    int32_t AC = A*C;
    int32_t AD_CB = A*D + C*B;
    uint32_t BD = B*D;

    int32_t product_hi = AC + (AD_CB >> w); // product_hi stands for the interger part of final result

    uint32_t ad_cb_temp = AD_CB << (32-w); // get the fraction part of AD_CB
    uint32_t product_lo = BD + ad_cb_temp;  //product_lo stands for the fraction part of final result

    if (product_lo < BD || product_lo < ad_cb_temp) { // check if product_lo overflow
        product_hi++;
    } 

    // The upper part bits should all be the same (including the sign).
    if (product_hi >> 31 != product_hi >> (31-w)) {
        printf("Overflow in fix_mul(), please use other bitwidth \n");
        return fix_overflow;
    }
    // combine interger part and fraction part
    return (product_hi << (w)) | (product_lo >> (32-w));
}

int test_mul(void)
{
    // test cases
    float a = 0.267f; //0.50f;//1.267f; //-1.267f;//-1.267f;
    float b = 0.849f; //0.25f;//1.849f; //1.849f; //-1.849f;
    for (int w = 1; w < 28; w++) {
        fix_t aa = fix_rconst(a, w);
        fix_t bb = fix_rconst(b, w);
        fix_t result = fix_mul(aa, bb, w);
        float fresult = fix_to_float(result, w);
        printf("fix_rconst(%f, %i) = %i, fix_rconst(%f, %i) = %i, result = %i, fresult=%f \n", a, w, aa, b, w, bb, result, fresult);
    }
    return 0;
}

int main()
{
    test_mul();
    // system("pause"); 
    return 0;
}

Vous pouvez obtenir le code en ligne aquí . Mais les résultats du test ne sont pas corrects sauf pour la largeur de 16 bits, le résultat attendu est d'environ 0,267*0,849=0,226683, il est acceptable qu'il y ait une erreur mineure pour une petite largeur de bit fractionnaire, c'est-à-dire une précision plus faible :

fix_rconst(0.267000, 1) = 1, fix_rconst(0.849000, 1) = 2, result = 1, fresult=0.500000 
fix_rconst(0.267000, 2) = 1, fix_rconst(0.849000, 2) = 3, result = 0, fresult=0.000000 
fix_rconst(0.267000, 3) = 2, fix_rconst(0.849000, 3) = 7, result = 0, fresult=0.000000 
fix_rconst(0.267000, 4) = 4, fix_rconst(0.849000, 4) = 14, result = 0, fresult=0.000000 
fix_rconst(0.267000, 5) = 9, fix_rconst(0.849000, 5) = 27, result = 0, fresult=0.000000 
fix_rconst(0.267000, 6) = 17, fix_rconst(0.849000, 6) = 54, result = 0, fresult=0.000000 
fix_rconst(0.267000, 7) = 34, fix_rconst(0.849000, 7) = 109, result = 0, fresult=0.000000 
fix_rconst(0.267000, 8) = 68, fix_rconst(0.849000, 8) = 217, result = 0, fresult=0.000000 
fix_rconst(0.267000, 9) = 137, fix_rconst(0.849000, 9) = 435, result = 0, fresult=0.000000 
fix_rconst(0.267000, 10) = 273, fix_rconst(0.849000, 10) = 869, result = 0, fresult=0.000000 
fix_rconst(0.267000, 11) = 547, fix_rconst(0.849000, 11) = 1739, result = 0, fresult=0.000000 
fix_rconst(0.267000, 12) = 1094, fix_rconst(0.849000, 12) = 3478, result = 3, fresult=0.000732 
fix_rconst(0.267000, 13) = 2187, fix_rconst(0.849000, 13) = 6955, result = 29, fresult=0.003540 
fix_rconst(0.267000, 14) = 4375, fix_rconst(0.849000, 14) = 13910, result = 232, fresult=0.014160 
fix_rconst(0.267000, 15) = 8749, fix_rconst(0.849000, 15) = 27820, result = 1856, fresult=0.056641 
fix_rconst(0.267000, 16) = 17498, fix_rconst(0.849000, 16) = 55640, result = 14855, fresult=0.226669 
fix_rconst(0.267000, 17) = 34996, fix_rconst(0.849000, 17) = 111280, result = 118846, fresult=0.906723 
fix_rconst(0.267000, 18) = 69992, fix_rconst(0.849000, 18) = 222560, result = 164338, fresult=0.626900 
fix_rconst(0.267000, 19) = 139985, fix_rconst(0.849000, 19) = 445121, result = 266201, fresult=0.507738 
fix_rconst(0.267000, 20) = 279970, fix_rconst(0.849000, 20) = 890241, result = 32390, fresult=0.030890 
fix_rconst(0.267000, 21) = 559940, fix_rconst(0.849000, 21) = 1780482, result = 259120, fresult=0.123558 
fix_rconst(0.267000, 22) = 1119879, fix_rconst(0.849000, 22) = 3560964, result = 2069485, fresult=0.493404 
fix_rconst(0.267000, 23) = 2239758, fix_rconst(0.849000, 23) = 7121928, result = 8167272, fresult=0.973615 
fix_rconst(0.267000, 24) = 4479517, fix_rconst(0.849000, 24) = 14243856, result = 15062169, fresult=0.897775 
fix_rconst(0.267000, 25) = 8959033, fix_rconst(0.849000, 25) = 28487712, result = 19611502, fresult=0.584468 
fix_rconst(0.267000, 26) = 17918066, fix_rconst(0.849000, 26) = 56975424, result = 22674290, fresult=0.337873 
fix_rconst(0.267000, 27) = 35836132, fix_rconst(0.849000, 27) = 113950848, result = 47176592, fresult=0.351493 

---------------------- MISE À JOUR --------------------------

J'ai trouvé une solution provisoire. L'erreur vient de la partie fractionnaire, c'est-à-dire B et D. La solution consiste à insérer le code suivant :

    uint32_t B = (inArg0 & MASK_M(w)), D = (inArg1 & MASK_M(w));

    // Inserted code here
    if(w < 16){
        B <<= (16-w);
        D <<= (16-w);
    }else{
        B >>= (w-16);
        D >>= (w-16);
    }

    int32_t AC = A*C;

mais cette solution de contournement ne fonctionne pas si a = 5.567, b = 2.7835 .

Answer 1

// w is fractional bitwidth
// separate interger part and fraciton part
int32_t A = (inArg0 >> w), C = (inArg1 >> w);
uint32_t B = (inArg0 & MASK_M(w)), D = (inArg1 & MASK_M(w));

Le but de ce code est de diviser les valeurs de 32 bits en deux parties. Ce n'est qu'une coïncidence si, dans libfixmath, le nombre de bits fractionnaires est également de 16. Même si vous utilisez un emplacement de point fixe différent, la multiplication réelle peut rester identique.

Voici une image adaptée d'un de mes travaux antérieurs. Le montant du décalage w n'affecte que l'utilisation du résultat final, et non les valeurs A/B/C/D. Il convient de noter que lorsque w n'est pas 16, le résultat ne sera pas proprement aligné sur les limites du mot et devra être décalé.

Answer 2

Si ce type fixedpt y fixedptd sont respectivement int32_t y int64_t La mise en œuvre simpliste fonctionne pour toutes les largeurs de bits des fractions, sauf peut-être pour la méthode d'arrondi. Vous pouvez essayer cette variante :

fixedpt fixedpt_mul(fixedpt A, fixedpt B)
{
    return (((fixedptd)A * (fixedptd)B + (1 << (FIXEDPT_FBITS - 1))) >> FIXEDPT_FBITS);
}

La norme C impose le type long long doit avoir au moins 63 bits de valeur, de sorte que les éléments ci-dessus devraient être pris en charge par la cible si elle est conforme à une version, même ancienne, de la norme.

Si vous ne souhaitez pas utiliser la multiplication sur 64 bits, vous pouvez utiliser la multiplication 16x16 -> 32 du code affiché pour l'émuler, puis arrondir et décaler le résultat final vers la droite de FIXEDPT_FBITS bits, ce qui nécessitera quelques ajustements en fonction de la présence ou non de FIXEDPT_FBITS >= 16 ou non.

Answer 3

Afin d'implémenter des calculs décimaux en virgule fixe, disons de N chiffres après la virgule, vous pouvez considérer que les nombres sont simplement des nombres entiers (avec l'unité au plus petit nombre décimal subinteger) lorsque vous additionnez/soustrayez deux nombres, vous n'avez aucun ajustement à faire, lorsque vous multipliez deux nombres, vous devez multiplier chaque résultat par 10^N (comme les deux nombres contribuent, chacun, avec une puissance de 10^N, vous devez diviser par l'un d'entre eux, pour que la virgule décimale soit à la bonne place). Cela vous obligera à multiplier le résultat par 10^N lors de la division (les deux puissances de dix s'annulent lors de la division, vous devez donc fournir un facteur pour placer le point décimal au bon endroit) et, lors de l'impression, il vous suffit de placer un point décimal à N places avant le dernier chiffre décimal (même s'il s'agit d'un zéro).

Conclusion, utilisez des entiers simples, en prenant soin, pour les multiplications, de diviser le résultat (après la multiplication) par 10^N, et de multiplier le quotient d'une division par 10^N. Bien sûr, comme l'unité de mesure est maintenant 10^n fois plus petite, vous aurez plus facilement des débordements, donc utilisez toujours des entiers complets de 64 bits.

Vous pouvez, au cas où vous auriez besoin de valeurs supérieures à ce que l'arithmétique des entiers 64 bits peut vous donner, utiliser une bibliothèque multiprécision comme gmp ou similaire (toutes les bibliothèques multiprécision gèrent les types décimaux à virgule fixe, normalement).