Sélectionner de bonnes premières estimations pour Goldschmidt division

StackOverflow https://stackoverflow.com/questions/2661541

  •  27-09-2019
  •  | 
  •  

Question

Je suis le calcul de fixedpoint inverses dans Q22.10 avec Goldschmidt division pour l'utiliser dans mon logiciel rasterizer sur les BRAS.

Cela se fait simplement en paramètre le numérateur 1, je.e le numérateur devient le scalaire sur la première itération.Pour être honnête, je suis un peu à la suite de la wikipédia algorithme à l'aveuglette ici.L'article dit que si le dénominateur est réduite dans la demi-intervalle ouvert (0.5, 1.0], une première estimation peut être fondée sur le seul dénominateur:Soit F l'estimation des scalaires et D être le dénominateur, alors F = 2 - D.

Mais en faisant cela, je perds beaucoup de précision.Disons que si je veux trouver la réciproque de 512.00002 f.Pour le numéro, je perds 10 bits de précision dans la fraction de la partie, qui est décalée en dehors.Donc, mes questions sont les suivantes:

  • Est-il un moyen de choisir une meilleure estimation qui ne nécessite pas de normalisation?Pourquoi?Pourquoi pas?Une preuve mathématique de ce qui est ou n'est pas possible serait génial.
  • Aussi, est-il possible de pré-calculer les premières estimations de sorte que la série converge plus rapidement?Pour l'instant, elle converge après la 4ème itération en moyenne.Sur les BRAS c'est environ ~50 cycles pire des cas, et qui n'est pas en prenant de l'émulation de la clozapine/bsr en compte, ni de la mémoire, les recherches.Si c'est possible, je voudrais savoir si vous limiterez ainsi l'erreur, et de combien.

Voici mon cas de test.Note:La mise en œuvre des logiciels de clz sur la ligne 13 est de mon post ici.Vous pouvez le remplacer par une valeur intrinsèque si vous le souhaitez. clz doit retourner le nombre de zéros, et 32 pour la valeur 0.

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

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)
    D >>= bitpos;
  else
    D <<= (-bitpos);

  /* (FNi / FDi) == (FN(i+1) / FD(i+1)) */
  /* F = 2 - D */
  F = (2ULL<<BASE) - D;
  /* N = F for the first iteration, because the numerator is simply 1.
     So don't waste a 64-bit UMULL on a multiply with 1 */
  N = F;
  D = ((unsigned long long)D*F)>>BASE;

  while(1){
    DPREV = D;
    F = (2<<(BASE)) - D;
    D = ((unsigned long long)D*F)>>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)
    N >>= 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;
}
Était-ce utile?

La solution

Je ne pouvais pas résister à passer une heure à votre problème...

Cet algorithme est décrit dans la section 5.5.2 de "Arithmetique des ordinateurs" par Jean-Michel Muller (en français).Il est en fait un cas particulier de Newton itérations avec 1 comme point de départ.Le livre donne une formulation simple de l'algorithme pour calculer N/D, avec D normalisée 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 sera suffisant.Vous pouvez également faire une itération jusqu'à ce que e devient trop petite pour modifier Q.

La normalisation est utilisé, car il représente le nombre maximum de bits significatifs dans le résultat.Il est aussi plus facile de calculer l'erreur et le nombre d'itérations nécessaires quand les entrées sont dans une aire de répartition connue.

Une fois votre valeur d'entrée est normalisée, vous n'avez pas besoin de s'embêter avec la valeur de BASE jusqu'à ce que vous avez l'inverse.Il vous suffit de disposer d'un nombre de 32 bits X normalisé dans la gamme 0x80000000 à 0xFFFFFFFF, et de calculer une approximation de Y=2^64/X (Y est au plus 2^33).

Cette simplifiée de l'algorithme peut être mis en œuvre pour votre Q22.10 représentation comme suit:

// Fixed point inversion
// EB Apr 2010

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

// 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<<BASE); }
uint32 toFP(double x) { return (int)floor(0.5+x*(1<<BASE)); }

// Return inverse of FP
uint32 inverse(uint32 fp)
{
  if (fp == 0) return (uint32)-1; // invalid

  // Shift FP to have the most significant bit set
  int shl = 0; // normalization shift
  uint32 nfp = fp; // normalized FP
  while ( (nfp & 0x80000000) == 0 ) { nfp <<= 1; shl++; } // use "clz" instead

  uint64 q = 0x100000000ULL; // 2^32
  uint64 e = 0x100000000ULL - (uint64)nfp; // 2^32-NFP
  int i;
  for (i=0;i<4;i++) // iterate
    {
      // Both multiplications are actually
      // 32x32 bits truncated to the 32 high bits
      q += (q*e)>>(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<Q<=2^33
  return (uint32)(q>>(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 de 32x32->64 bits.E va devenir de plus en plus petites et s'inscrit d'abord sur 32 bits.Q sera toujours sur 34 bits.Nous ne prenons que le haut de 32 bits des produits.

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

EDIT.Comme suite à mon commentaire, voici une deuxième version avec un implicite de la 32-ème bit sur Q.Les deux 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
}

Autres conseils

Quelques idées pour vous, mais aucun qui permettent de résoudre votre problème directement comme indiqué.

  1. Pourquoi cette algo de division?La plupart divise que j'ai vu dans les BRAS d'utiliser certains variante de 
    
      adcs hi, den, hi, lsl #1
      subcc hi, hi, den
      adcs lo, lo, lo

répété n bits fois avec un fichier binaire de recherche hors de la clozapine à déterminer par où commencer.C'est assez dang 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.Ça va être un peu plus lent, mais vous pouvez ajouter/adc ou sub/sbc pour déplacer les valeurs dans les registres.mul/mla sont également conçus pour ce genre de travail.

Encore une fois, pas de réponse directe pour vous, mais peut-être quelques idées pour aller de l'avant cette.Voir le BRAS de code serait probablement m'aider un peu aussi.

Mads, vous n'êtes pas perdre de précision à tous.Lorsque vous divisez 512.00002 f par 2^10, il vous suffit de diminuer la puissance de votre nombre à virgule flottante par 10.Mantisse reste le même.Bien sûr, à moins que l'exposant atteint sa valeur minimale, mais qui ne devrait pas arriver puisque vous êtes mise à l'échelle (0.5, 1].

EDIT:Ok, donc vous êtes à l'utilisation d'un point décimal.Dans ce cas, vous devriez permettre à une représentation différente du dénominateur dans votre algorithme.La valeur de D est de (0.5, 1], non seulement au début, mais tout au long de l'ensemble du calcul (il est facile de prouver que x * (2-x) < 1 pour x < 1).De sorte que vous devrait représenter le dénominateur avec point décimal à la base = 32.De cette façon, vous aurez 32 bits de précision de tous les temps.

EDIT:Pour mettre en œuvre ce que vous aurez à modifier les lignes suivantes dans votre code:

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

Également à la fin vous aurez à déplacer N pas par bitpos mais certains différents de la valeur que je suis trop paresseux pour comprendre tout de suite :).

Licencié sous: CC-BY-SA avec attribution
Non affilié à StackOverflow
scroll top