Scegliere buone prime stime per la divisione Goldschmidt

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

  •  27-09-2019
  •  | 
  •  

Domanda

Sto calcolando i reciproci in virgola fissa in Q22.10 con Divisione Goldschmidt da utilizzare nel mio rasterizzatore software su ARM.

Questo viene fatto semplicemente impostando il numeratore su 1, ovvero il numeratore diventa lo scalare alla prima iterazione.Ad essere onesti, sto seguendo ciecamente l'algoritmo di Wikipedia qui.L'articolo afferma che se il denominatore viene ridimensionato nell'intervallo semiaperto (0,5, 1,0], una buona prima stima può essere basata solo sul denominatore:Sia F lo scalare stimato e D il denominatore, quindi F = 2 - D.

Ma così facendo perdo molta precisione.Diciamo se voglio trovare il reciproco di 512.00002f.Per ridurre il numero, perdo 10 bit di precisione nella parte frazionaria, che viene spostata.Quindi, le mie domande sono:

  • Esiste un modo per scegliere una stima migliore che non richieda la normalizzazione?Perché?Perché no?Una prova matematica del perché ciò sia o meno possibile sarebbe grandiosa.
  • Inoltre, è possibile precalcolare le prime stime in modo che la serie converga più velocemente?In questo momento, converge in media dopo la quarta iterazione.Su ARM si tratta del caso peggiore di circa ~ 50 cicli e non tiene conto dell'emulazione di clz/bsr né delle ricerche di memoria.Se è possibile, vorrei sapere se così facendo aumenta l'errore e di quanto.

Ecco il mio caso di prova.Nota:L'implementazione del software di clz sulla linea 13 è dal mio post Qui.Se lo desideri, puoi sostituirlo con un intrinseco. clz dovrebbe restituire il numero di zeri iniziali e 32 per il valore 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;
}
È stato utile?

Soluzione

non ho potuto resistere spendere un'ora sul vostro problema ...

Questo algoritmo è descritta nella sezione 5.5.2 di "arithmétique des ordinateurs" di Jean-Michel Muller (in francese). In realtà è un caso speciale di iterazioni Newton con 1 come punto di partenza. Il libro fornisce una semplice formulazione dell'algoritmo per calcolare N / D, con D normalizzato nell'intervallo [1 / 2,1 [: 

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

Il numero di bit corretti raddoppia ad ogni iterazione. Nel caso di 32 bit, 4 iterazioni sarà sufficiente. È anche possibile iterare fino e diventa troppo piccolo per modificare Q.

La normalizzazione viene usato perché fornisce il numero massimo di bit significativi nel risultato. È anche più facile calcolare l'errore e il numero di iterazioni necessarie quando gli ingressi sono in un intervallo noto.

Una volta che il valore di ingresso viene normalizzato, non c'è bisogno di perdere tempo con il valore di base fino a quando si ha l'inverso. È sufficiente un numero a 32 bit X normalizzata nell'intervallo 0x80000000 su 0xFFFFFFFF, e calcolare un'approssimazione di Y = 2 ^ 64 / X (Y è al massimo 2 ^ 33).

Questo algoritmo semplificato può essere implementato per la rappresentazione Q22.10 come segue: 

// 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);
}

Come indicato nel codice, le moltiplicazioni non sono piene 32x32-> 64 bit. E diventa sempre più piccoli e si adatta inizialmente su 32 bit. Q sarà sempre di 34 bit. Prendiamo soltanto gli alti 32 bit dei prodotti.

La derivazione 64-2*BASE-shl è lasciata come esercizio per il lettore :-). Se diventa 0 o negativo, il risultato non è rappresentabile (il valore di ingresso è troppo piccolo).

EDIT. Come un follow-up al mio commento, ecco una seconda versione con un implicito a 32-esimo bit su Q. Sia E e Q sono ora memorizzati su 32 bit: 

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
}

Altri suggerimenti

Un paio di idee per voi, anche se nessuno che risolve il problema direttamente, come indicato.

  1. Perchè questo algo per la divisione? La maggior parte delle divisioni che ho visto in ARM usano una certa variante di 
    
      adcs hi, den, hi, lsl #1
      subcc hi, hi, den
      adcs lo, lo, lo

ripetuti n bit volte con un binario Cerca off del CLZ per determinare dove iniziare. Questo è veloce piuttosto dang.

  1. Se la precisione è un grosso problema, non è limitato a 32/64 bit per la rappresentazione in virgola fissa. Sarà un po 'più lento, ma si può fare add / ADC o sub / SBC ai valori di movimento attraverso i registri. mul / mla sono inoltre progettati per questo tipo di lavoro.

Di nuovo, non risposte dirette per voi, ma forse un paio di idee per andare avanti questo. Vedendo il codice ARM reale sarebbe probabilmente mi aiuterà un po 'pure.

Mads, non stai perdendo assolutamente la precisione.Quando dividi 512.00002f per 2^10, diminuisci semplicemente l'esponente del tuo numero in virgola mobile di 10.Mantissa rimane la stessa.Ovviamente, a meno che l'esponente non raggiunga il suo valore minimo, ma ciò non dovrebbe accadere poiché stai scalando a (0,5, 1].

MODIFICARE:Ok, quindi stai usando un punto decimale fisso.In tal caso dovresti consentire una diversa rappresentazione del denominatore nel tuo algoritmo.Il valore di D è da (0,5, 1] ​​non solo all'inizio ma durante l'intero calcolo (è facile dimostrare che x * (2-x) < 1 per x < 1).Quindi dovresti rappresentare il denominatore con punto decimale in base = 32.In questo modo avrai sempre 32 bit di precisione.

MODIFICARE:Per implementarlo dovrai modificare le seguenti righe del tuo codice:

  //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;

Inoltre alla fine dovrai spostare N non tramite bitpos ma con un valore diverso che sono troppo pigro per capire in questo momento :).

Autorizzato sotto: CC-BY-SA insieme a attribuzione
Non affiliato a StackOverflow
scroll top