Effiziente Berechnung von 2 ** 64/Divisor über schnellen Fließkomma-Reziprok

Question

Ich untersuche derzeit Wege, um die schnelle Gleitkomma-Reziprokfähigkeit verschiedener moderner Prozessoren zu verwenden, um eine Startapproximation für ein 64-Bit zu berechnen vorzeichenlose ganzzahlige Division basierend auf Festkomma-Newton-Raphson-Iterationen. Sie erfordert eine möglichst genaue Berechnung des Divisors, wobei die anfängliche Approximation kleiner oder gleich dem mathematischen Ergebnis sein muss, basierend auf den Anforderungen der folgenden Festkomma-Iterationen. Dies bedeutet, dass diese Berechnung zu unterschätzen ist. Im Moment habe ich den folgenden Code, der gut funktioniert, basierend auf umfangreichen Tests:

#include <stdint.h> // import uint64_t 
#include <math.h> // import nextafterf() 

uint64_t divisor, recip; 
float r, s, t; 

t = uint64_to_float_ru (divisor); // ensure t >= divisor 
r = 1.0f/t; 
s = 0x1.0p64f * nextafterf (r, 0.0f); 
recip = (uint64_t)s; // underestimate of 2**64/divisor 

Während dieser Code funktioniert, es ist nicht gerade schnell auf den meisten Plattformen. Eine offensichtliche Verbesserung, die ein wenig maschinenspezifischen Code erfordert, besteht darin, die Division r = 1.0f/t durch einen Code zu ersetzen, der von der Hardware bereitgestellte schnelle Fließkomma-Reziprokwerte verwendet. Dies kann durch Iteration erhöht werden, um ein Ergebnis zu erzeugen, das innerhalb von 1 ul des mathematischen Ergebnisses liegt, so dass im Kontext des bestehenden Codes eine Unterschätzung erzeugt wird. Eine Beispielimplementierung für x86_64 wäre:

#include <xmmintrin.h> 
/* Compute 1.0f/a almost correctly rounded. Halley iteration with cubic convergence */ 
inline float fast_recip_f32 (float a) 
{ 
    __m128 t; 
    float e, r; 
    t = _mm_set_ss (a); 
    t = _mm_rcp_ss (t); 
    _mm_store_ss (&r, t); 
    e = fmaf (r, -a, 1.0f); 
    e = fmaf (e, e, e); 
    r = fmaf (e, r, r); 
    return r; 
} 

Implementationen von nextafterf() sind in der Regel nicht die Leistung optimiert. Auf Plattformen, wo es Mittel gibt, um schnell eine IEEE 754 binary32 in eine int32 und umgekehrt über intrinsics float_as_int() und int_as_float(), können wir den Einsatz von nextafterf() und Skalierung kombinieren neu interpretiert wie folgt: Unter der Annahme,

s = int_as_float (float_as_int (r) + 0x1fffffff); 

diese Ansätze Auf einer gegebenen Plattform möglich, lassen uns die Umbauten zwischen float und uint64_t als Haupthindernisse. Die meisten Plattformen bieten keine Anweisung, die eine Konvertierung von uint64_t zu float mit statischem Rundungsmodus (hier: in Richtung positive Unendlichkeit = aufwärts) durchführt, und einige bieten keine Anweisungen zum Konvertieren zwischen uint64_t und Gleitkommatypen ein Leistungsengpass.

t = uint64_to_float_ru (divisor); 
r = fast_recip_f32 (t); 
s = int_as_float (float_as_int (r) + 0x1fffffff); 
recip = (uint64_t)s; /* underestimate of 2**64/divisor */ 

Eine tragbare, aber langsam, Implementierung von uint64_to_float_ru verwendet dynamische Änderungen FPU-Rundungsmodus:

#include <fenv.h> 
#pragma STDC FENV_ACCESS ON 

float uint64_to_float_ru (uint64_t a) 
{ 
    float res; 
    int curr_mode = fegetround(); 
    fesetround (FE_UPWARD); 
    res = (float)a; 
    fesetround (curr_mode); 
    return res; 
} 

ich verschiedene Splitting und Bit-Fummel Ansätze ausgesehen habe in dem Conversions zu beschäftigen (zB tun die Rundung auf der Ganzzahlseite, dann verwenden Sie eine normale Umwandlung zu float, die den Rundungsmodus IEEE 754 Round-to-Nearest-or-even verwendet, aber der dabei entstehende Overhead macht diese Berechnung über schnelle reziproke Gleitkommazahl aus einer Performance Perspektive. So wie es aussieht, scheint es mir besser zu sein, eine Anfangsapproximation zu erzeugen, indem ich eine klassische LUT mit Interpolation oder eine Festkomma-Polynomapproximation benutze, und folge diesen mit einem 32-Bit-Festkomma-Newton-Raphson-Schritt.

Gibt es Möglichkeiten, die Effizienz meines derzeitigen Ansatzes zu verbessern? Portable und semi-portable Wege mit intrinsic für spezifische Plattformen wären von Interesse (insbesondere für x86 und ARM als die derzeit dominierenden CPU-Architekturen). Kompilieren für x86_64 unter Verwendung des Intel-Compilers bei sehr hoher Optimierung (/O3 /QxCORE-AVX2 /Qprec-div-) erfordert die Berechnung der anfänglichen Approximation mehr Anweisungen als die Iteration, die ungefähr 20 Anweisungen benötigt. Unten ist der vollständige Divisionscode als Referenz angegeben, der die Approximation im Kontext zeigt.

uint64_t udiv64 (uint64_t dividend, uint64_t divisor) 
{ 
    uint64_t temp, quot, rem, recip, neg_divisor = 0ULL - divisor; 
    float r, s, t; 

    /* compute initial approximation for reciprocal; must be underestimate! */ 
    t = uint64_to_float_ru (divisor); 
    r = 1.0f/t; 
    s = 0x1.0p64f * nextafterf (r, 0.0f); 
    recip = (uint64_t)s; /* underestimate of 2**64/divisor */ 

    /* perform Halley iteration with cubic convergence to refine reciprocal */ 
    temp = neg_divisor * recip; 
    temp = umul64hi (temp, temp) + temp; 
    recip = umul64hi (recip, temp) + recip; 

    /* compute preliminary quotient and remainder */ 
    quot = umul64hi (dividend, recip); 
    rem = dividend - divisor * quot; 

    /* adjust quotient if too small; quotient off by 2 at most */ 
    if (rem >= divisor) quot += ((rem - divisor) >= divisor) ? 2 : 1; 

    /* handle division by zero */ 
    if (divisor == 0ULL) quot = ~0ULL; 

    return quot; 
} 

umul64hi() würde Karte im Allgemeinen auf ein plattformspezifische intrinsisch, oder ein bisschen von Inline-Assembler-Code. Auf x86_64 verwende ich zur Zeit diese Implementierung:

inline uint64_t umul64hi (uint64_t a, uint64_t b) 
{ 
    uint64_t res; 
    __asm__ (
     "movq %1, %%rax;\n\t" // rax = a 
     "mulq %2;\n\t"   // rdx:rax = a * b 
     "movq %%rdx, %0;\n\t" // res = (a * b)<63:32> 
     : "=rm" (res) 
     : "rm"(a), "rm"(b) 
     : "%rax", "%rdx"); 
    return res; 
} 

Answer 1

Diese Lösung kombiniert zwei Ideen:

• Sie zu Floating-Point, indem Sie einfach neu interpretiert die Bits als Floating-Point und Subtrahieren einer konstanten, solange die umwandeln kann Nummer liegt in einem bestimmten Bereich. Also füge eine Konstante hinzu, reinterpretiere und subtrahiere diese Konstante. Dies ergibt ein abgeschnittenes Ergebnis (das daher immer kleiner oder gleich dem gewünschten Wert ist).
• Sie können den reziproken Wert annähern, indem Sie sowohl den Exponenten als auch die Mantisse negieren. Dies kann erreicht werden, indem die Bits als int interpretiert werden.

Option 1 funktioniert hier nur in einem bestimmten Bereich, daher prüfen wir den Bereich und passen die verwendeten Konstanten an. Dies funktioniert in 64 Bits, da der gewünschte Float nur 23 Bits Genauigkeit hat.

Das Ergebnis in diesem Code wird doppelt sein, aber die Konvertierung in Float ist trivial, und kann auf den Bits oder direkt erfolgen, je nach Hardware.

Danach sollten Sie die Newton-Raphson Iteration (en) tun.

Ein Großteil dieses Codes konvertiert einfach zu magischen Zahlen.

double              
u64tod_inv(uint64_t u64) {         
    __asm__("#annot0");          
    union {              
    double f;             
    struct {             
     unsigned long m:52; // careful here with endianess  
     unsigned long x:11;          
     unsigned long s:1;          
    } u64;             
    uint64_t u64i;           
    } z,              
     magic0 = { .u64 = { 0, (1<<10)-1 + 52, 0 } },   
     magic1 = { .u64 = { 0, (1<<10)-1 + (52+12), 0 } }, 
     magic2 = { .u64 = { 0, 2046, 0 } };     

    __asm__("#annot1");          
    if(u64 < (1UL << 52UL)) {        
    z.u64i = u64 + magic0.u64i;        
    z.f -= magic0.f;          
    } else {             
    z.u64i = (u64 >> 12) + magic1.u64i;      
    z.f -= magic1.f;          
    }               
    __asm__("#annot2");          

    z.u64i = magic2.u64i - z.u64i;        

    return z.f;             
}                

das Kompilieren auf einem Intel-Core-7 gibt eine Reihe von Anweisungen (und einem Zweig), aber natürlich keine Multiplikationen oder überhaupt trennt. Wenn die Umwandlungen zwischen int und double schnell sind, sollte dies ziemlich schnell laufen.

Ich vermute, Schwimmer (mit nur 23 Bit Genauigkeit) wird mehr als 1 oder 2 Newton-Raphson-Iterationen erfordern die Genauigkeit, die Sie wollen zu bekommen, aber ich habe nicht die Mathematik gemacht ...