De bonnes premières estimations pour la division Goldschmidt

Je calcule des inverses de point fixe en Q22.10 avec la division Goldschmidt pour les utiliser dans mon logiciel de tramage logiciel sur ARM.

Cela se fait en réglant simplement le numérateur sur 1, c’est-à-dire que le numérateur devient le scalaire à la première itération. Pour être honnête, je suis en quelque sorte de suivre aveuglément l’algorithme de Wikipedia. L’article indique que si le dénominateur est mis à l’échelle dans la plage semi-ouverte (0.5, 1.0], une bonne première estimation peut être basée sur le seul dénominateur: Soit F le scalaire estimé et D le dénominateur, alors F = 2 – RÉ.

Mais en faisant cela, je perds beaucoup de précision. Dites si je veux trouver l’inverse de 512.00002f. Afin de réduire le nombre, je perds 10 bits de précision dans la partie fraction qui est décalée. Donc, mes questions sont:

  • Y at-il un moyen de choisir une meilleure estimation qui ne nécessite pas de normalisation? Pourquoi? Pourquoi pas? Une preuve mathématique expliquant pourquoi cela est possible ou non serait formidable.
  • Aussi, est-il possible de pré-calculer les premières estimations pour que la série converge plus rapidement? Pour le moment, il converge en moyenne après la 4ème itération. Sur ARM, il s’agit d’environ 50 cycles dans le cas le plus défavorable, ce qui ne tient pas compte de l’émulation de clz / bsr, ni des recherches dans la mémoire. Si c’est possible, j’aimerais savoir si cela augmente l’erreur et de combien.

Voici mon test. Remarque: la mise en œuvre logicielle de clz sur la ligne 13 est tirée de mon article ici . Vous pouvez le remplacer par un insortingnsèque si vous le souhaitez. clz doit renvoyer le nombre de zéros non clz et 32 ​​pour la valeur 0.

 #include  #include  const unsigned int BASE = 22ULL; static unsigned int divfp(unsigned int val, int* iter) { /* Numerator, denominator, estimate scalar and previous denominator */ unsigned long long N,D,F, DPREV; int bitpos; *iter = 1; D = val; /* Get the shift amount + is right-shift, - is left-shift. */ bitpos = 31 - clz(val) - BASE; /* Normalize into the half-range (0.5, 1.0] */ if(0 >= bitpos; else D <<= (-bitpos); /* (FNi / FDi) == (FN(i+1) / FD(i+1)) */ /* F = 2 - D */ F = (2ULL<>BASE; while(1){ DPREV = D; F = (2<>BASE; /* Bail when we get the same value for two denominators in a row. This means that the error is too small to make any further progress. */ if(D == DPREV) break; N = ((unsigned long long)N*F)>>BASE; *iter = *iter + 1; } if(0 >= bitpos; else N <<= (-bitpos); return N; } int main(int argc, char* argv[]) { double fv, fa; int iter; unsigned int D, result; sscanf(argv[1], "%lf", &fv); D = fv*(double)(1<<BASE); result = divfp(D, &iter); fa = (double)result / (double)(1UL << BASE); printf("Value: %8.8lf 1/value: %8.8lf FP value: 0x%.8X\n", fv, fa, result); printf("iteration: %d\n",iter); return 0; } 

    Je ne pouvais pas résister à passer une heure sur ton problème …

    Cet algorithme est décrit dans la section 5.5.2 de “Arithmetique des ordinateurs” de Jean-Michel Muller. Il s’agit en fait d’un cas particulier d’itérations de Newton avec 1 comme sharepoint départ. Le livre donne une formulation simple de l’algorithme pour calculer N / D, avec D normalisé dans l’intervalle [1 / 2,1 [:

     e = 1 - D Q = N repeat K times: Q = Q * (1+e) e = e*e 

    Le nombre de bits corrects double à chaque itération. Dans le cas de 32 bits, 4 itérations suffiront. Vous pouvez également itérer jusqu’à ce que e devienne trop petit pour modifier Q

    La normalisation est utilisée car elle fournit le nombre maximal de bits significatifs dans le résultat. Il est également plus facile de calculer l’erreur et le nombre d’itérations nécessaires lorsque les entrées sont dans une plage connue.

    Une fois que votre valeur d’entrée est normalisée, vous n’avez pas besoin de vous préoccuper de la valeur de BASE avant d’avoir l’inverse. Vous avez simplement un nombre 32 bits X normalisé dans la plage 0x80000000 à 0xFFFFFFFF et calculez une approximation de Y = 2 ^ 64 / X (Y est au maximum 2 ^ 33).

    Cet algorithme simplifié peut être implémenté pour votre représentation Q22.10 comme suit:

     // Fixed point inversion // EB Apr 2010 #include  #include  // Number X is represented by integer I: X = I/2^BASE. // We have (32-BASE) bits in integral part, and BASE bits in fractional part #define BASE 22 typedef unsigned int uint32; typedef unsigned long long int uint64; // Convert FP to/from double (debug) double toDouble(uint32 fp) { return fp/(double)(1<>(uint64)32; e = (e*e)>>(uint64)32; printf("Q=0x%llx E=0x%llx\n",q,e); } // Here, (Q/2^32) is the inverse of (NFP/2^32). // We have 2^31<=NFP<2^32 and 2^32>(64-2*BASE-shl)); } int main() { double x = 1.234567; uint32 xx = toFP(x); uint32 yy = inverse(xx); double y = toDouble(yy); printf("X=%f Y=%f X*Y=%f\n",x,y,x*y); printf("XX=0x%08x YY=0x%08x XX*YY=0x%016llx\n",xx,yy,(uint64)xx*(uint64)yy); } 

    Comme indiqué dans le code, les multiplications ne sont pas complètes 32×32-> 64 bits. E deviendra de plus en plus petit et s’adapte initialement sur 32 bits. Q sera toujours sur 34 bits. Nous prenons uniquement les 32 bits les plus élevés des produits.

    La dérivation de 64-2*BASE-shl est laissée comme exercice au lecteur :-). S’il devient 0 ou négatif, le résultat n’est pas représentable (la valeur d’entrée est trop petite).

    MODIFIER. Pour faire suite à mon commentaire, voici une deuxième version avec un 32ème bit implicite sur Q. E et Q sont maintenant stockés sur 32 bits:

     uint32 inverse2(uint32 fp) { if (fp == 0) return (uint32)-1; // invalid // Shift FP to have the most significant bit set int shl = 0; // normalization shift for FP uint32 nfp = fp; // normalized FP while ( (nfp & 0x80000000) == 0 ) { nfp <<= 1; shl++; } // use "clz" instead int shr = 64-2*BASE-shl; // normalization shift for Q if (shr <= 0) return (uint32)-1; // overflow uint64 e = 1 + (0xFFFFFFFF ^ nfp); // 2^32-NFP, max value is 2^31 uint64 q = e; // 2^32 implicit bit, and implicit first iteration int i; for (i=0;i<3;i++) // iterate { e = (e*e)>>(uint64)32; q += e + ((q*e)>>(uint64)32); } return (uint32)(q>>shr) + (1<<(32-shr)); // insert implicit bit } 

    Quelques idées pour vous, mais aucune ne résout directement votre problème, comme indiqué.

    1. Pourquoi cet algo de division? La plupart des fractures que j’ai vu dans ARM utilisent une variante de
       adcs hi, den, hi, lsl #1 subcc hi, hi, den adcs lo, lo, lo 

    répété n bits fois avec une recherche binary hors du clz pour déterminer par où commencer. C’est assez rapide.

    1. Si la précision est un gros problème, vous n’êtes pas limité à 32/64 bits pour votre représentation en virgule fixe. Ce sera un peu plus lent, mais vous pouvez append / adc ou sub / sbc pour déplacer des valeurs d’un registre à l’autre. Les mul / mla sont également conçus pour ce type de travail.

    Encore une fois, pas de réponses directes pour vous, mais peut-être quelques idées pour aller de l’avant. Voir le code ARM m’aiderait probablement aussi un peu.

    Mads, vous ne perdez aucune précision. Lorsque vous divisez 512.00002f par 2 ^ 10, vous réduisez simplement l’exposant de votre nombre à virgule flottante de 10. Mantissa rest identique. Bien sûr, à moins que l’exposant n’atteigne sa valeur minimale, mais cela ne devrait pas arriver puisque vous passez à (0.5, 1].

    EDIT: Ok, vous utilisez donc un point décimal fixe. Dans ce cas, vous devez autoriser une représentation différente du dénominateur dans votre algorithme. La valeur de D est comprise entre (0.5 et 1) non seulement au début, mais tout au long du calcul (il est facile de prouver que x * (2-x) <1 pour x <1). Vous devez donc représenter le dénominateur avec un nombre décimal. point at base = 32. De cette façon, vous aurez tout le temps 32 bits de précision.

    EDIT: Pour implémenter ceci, vous devrez changer les lignes suivantes de votre code:

      //bitpos = 31 - clz(val) - BASE; bitpos = 31 - clz(val) - 31; ... //F = (2ULL<>BASE; F = -D; N = F >> (31 - BASE); D = ((unsigned long long)D*F)>>31; ... //F = (2<<(BASE)) - D; //D = ((unsigned long long)D*F)>>BASE; F = -D; D = ((unsigned long long)D*F)>>31; ... //N = ((unsigned long long)N*F)>>BASE; N = ((unsigned long long)N*F)>>31; 

    Aussi, à la fin, vous devrez déplacer N non pas par bitpos mais par une valeur différente que je suis trop paresseux pour comprendre maintenant :).