Picking gute erste Schätzungen für die Teilung der Firma Goldschmidt

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

  •  27-09-2019
  •  | 
  •  

Frage

Ich Berechnung FixedPoint reziproke in Q22.10 mit der Firma Goldschmidt Division für die Verwendung in meiner Software Rasterizer auf ARM.

Dies ist nur durch die Einstellung der Zähler auf 1 durchgeführt, das heißt der Zähler die skalare auf der ersten Iteration wird. Um ehrlich zu sein, bin ich irgendwie nach der wikipedia Algorithmus blind hier. Der Artikel sagt, dass, wenn der Nenner in dem halboffenen Bereich skaliert wird (0,5, 1,0], eine gute erste Schätzung kann allein auf dem Nenner basieren: Sei F der geschätzte Skalar und D der Nenner, dann F = 2 sein werden - D.

Aber wenn dies zu tun, ich habe eine Menge Präzision verlieren. Sprich, wenn ich die reziproke 512.00002f finden möchten. Um die Zahl nach unten zu skalieren, verlieren I 10 Bit Genauigkeit in dem Bruchteil, der aus verschoben wird. Also, meine Fragen sind:

  • Gibt es eine Möglichkeit, eine bessere Schätzung zu holen, die nicht Normalisierung erfordert? Warum? Warum nicht? Ein mathematischer Beweis, warum dies oder nicht möglich wäre toll.
  • Auch ist es möglich, die ersten Schätzungen, so dass die Serie konvergiert schneller vorab berechnen? Gerade jetzt, sie konvergiert nach der 4. Iteration im Durchschnitt. Auf ARM ist dies etwa ~ 50 Zyklen schlimmster Fall, und das ist nicht Emulation von clz / bsr berücksichtigt, noch Speicher-Lookups. Wenn es möglich ist, würde ich gerne wissen, ob dabei den Fehler erhöht, und um wie viel.

Hier ist mein Testfall. Hinweis: Die Software-Implementierung von clz auf Linie 13 ist von meinem Posten hier . Sie können es mit einer intrinsischen ersetzen, wenn Sie wollen. clz sollte die Anzahl der führenden Nullen zurück und 32 für den Wert 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;
}
War es hilfreich?

Lösung

Ich konnte nicht wieder eine Stunde auf Ihrem Problem zu verbringen ...

Dieser Algorithmus ist in Abschnitt 5.5.2 von "arithmétique des ordinateurs" von Jean-Michel Muller (auf Französisch) beschrieben. Es ist eigentlich ein Sonderfall von Newton Iterationen mit 1 als Ausgangspunkt. Das Buch gibt eine einfache Formulierung des Algorithmus N / D, zu berechnen, mit D in Bereich normalisierten [1 / 2,1 [: 

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

Die Anzahl der korrekten Bits verdoppelt bei jeder Iteration. Im Fall von 32 Bit werden 4 Iterationen ausreichen. Sie können auch Iterierte bis e zu klein wird Q zu ändern.

Normalisierungs verwendet, da es die maximale Anzahl von Bits in dem Ergebnis liefert. Es ist auch einfacher, die Fehler und die Anzahl der Iterationen benötigt, um zu berechnen, wenn die Eingänge in einem bekannten Bereich liegen.

Wenn Sie Ihr Eingabewert normalisiert wird, Sie brauchen nicht mit dem Wert von BASE zu stören, bis Sie die inverse haben. Sie müssen lediglich eine 32-Bit-Zahl X im Bereich 0x80000000 auf 0xFFFFFFFF normalisiert und eine Annäherung an Y berechnen = 2 ^ 64 / X (Y ist höchstens 2 ^ 33).

Dieser vereinfachte Algorithmus für Ihre Q22.10 Darstellung umgesetzt werden können wie folgt: 

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

Wie im Code angegeben, sind die Multiplikationen nicht voll 32x32-> 64 Bit. E wird kleiner und kleiner werden und passt sich zunächst auf 32 Bit. Q wird immer auf 34 Bit betragen. Wir nehmen nur die hohen 32 Bits der Produkte.

Die Ableitung von 64-2*BASE-shl wird für den Leser als Übung :-). Wenn es 0 oder negativ wird, ist das Ergebnis nicht darstellbare (der Eingangswert zu klein ist).

EDIT. Als Follow-up zu meinem Kommentar, hier ist eine zweite Version mit einem impliziten 32-te Bit auf Q. Sowohl E und Q werden nun auf 32 Bit gespeichert: 

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
}

Andere Tipps

Ein paar Ideen für Sie, aber keine, die Ihr Problem lösen direkt, wie angegeben.

  1. Warum diese algo für die Division? Die meisten dividieren ich in ARM gesehen habe verwenden, um einige varient von 
    
      adcs hi, den, hi, lsl #1
      subcc hi, hi, den
      adcs lo, lo, lo

wiederholt n Bits mal mit einem binären dem clz sucht aus, um zu bestimmen, wo zu beginnen. Das ist ziemlich dingt schnell.

  1. Wenn Präzision ein großes Problem ist, sind Sie mit 32/64 Bits für Ihre Festpunktdarstellung beschränkt. Es wird etwas langsamer sein, aber Sie können zu bewegen Werte über Register Add / adc oder Sub / sbc tun. mul / mla ist auch für diese Art von Arbeit entwickelt.

Auch nicht direkte Antworten für Sie, aber vielleicht ein paar Ideen, diese vorwärts zu gehen. den tatsächlichen ARM-Code zu sehen, würde mir wahrscheinlich auch ein bisschen helfen.

Mads, Sie verlieren keine Präzision überhaupt. Wenn Sie teilen 512.00002f von 2 ^ 10, verringern Sie nur den Exponenten Ihrer Gleitkommazahl um 10. Mantissa gleich bleibt. Natürlich, es sei denn der Exponent seinen Minimalwert trifft aber, dass da Sie nicht passieren sollte ist Skalierung (0,5, 1].

EDIT: Ok, so dass Sie verwenden eine Festkomma. In diesem Fall sollten Sie eine andere Darstellung des Nenners in Ihrem Algorithmus erlauben. Der Wert von D ist aus (0,5, 1] ??nicht nur am Anfang, sondern während der gesamten Berechnung (es ist einfach, dass x * (2-x) <1 für x <1) zu beweisen. So ermitteln Sie den Nenner mit Nachkommastellen darstellen sollten Punkt an der Basis = 32. auf diese Weise können 32 Bit Genauigkeit die ganze Zeit haben wird.

EDIT: Um dies zu implementieren, werden Sie die folgenden Zeilen des Codes ändern müssen: 

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

Auch am Ende werden Sie N haben verschieben nicht durch bitpos aber einigen anderen Wert, den ich zu faul bin jetzt, um herauszufinden.)

Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top