Wie funktioniert diese Float-Quadratwurzel-Näherung?

Ich fand eine ziemlich seltsame, aber funktionierende Quadratwurzel-Näherung für floats; Ich verstehe es wirklich nicht. Kann mir jemand erklären warum dieser Code funktioniert?

float sqrt(float f)
{
    const int result = 0x1fbb4000 + (*(int*)&f >> 1);
    return *(float*)&result;   
}

Ich habe es ein bisschen getestet und es gibt Werte von aus std::sqrt() um etwa 1 bis 3 %. Ich kenne die Quake III schnelle inverse Quadratwurzel und ich denke, es ist hier etwas Ähnliches (ohne die Newton-Iteration), aber ich würde mich sehr über eine Erklärung freuen wie es funktioniert.

(Anmerkung: Ich habe es sowohl mit c als auch mit c++ gekennzeichnet, da es sich um gültigen (siehe Kommentare) C- und C++-Code handelt.)

(*(int*)&f >> 1) verschiebt die bitweise Darstellung von nach rechts f. Dies fast dividiert den Exponenten durch zwei, was ungefähr dem Ziehen der Quadratwurzel entspricht.¹

Warum fast? In IEEE-754 ist der eigentliche Exponent e-127.² Um dies durch zwei zu teilen, bräuchten wir e/2 – 64aber die obige Annäherung gibt uns nur e/2 – 127. Also müssen wir 63 zum resultierenden Exponenten addieren. Dies wird durch die Bits 30-23 dieser magischen Konstante (0x1fbb4000).

Ich würde mir vorstellen, dass die verbleibenden Bits der magischen Konstante ausgewählt wurden, um den maximalen Fehler über den Mantissenbereich oder so ähnlich zu minimieren. Es ist jedoch unklar, ob es analytisch, iterativ oder heuristisch bestimmt wurde.

Es ist erwähnenswert, dass dieser Ansatz etwas nicht portabel ist. Sie macht (mindestens) folgende Annahmen:

• Die Plattform verwendet IEEE-754 mit einfacher Genauigkeit für float.
• Die Endianität von float Darstellung.
• Dass Sie von undefiniertem Verhalten nicht betroffen sind, da dieser Ansatz gegen die strikten Aliasing-Regeln von C/C++ verstößt.

Daher sollte es vermieden werden, es sei denn, Sie sind sich sicher, dass es ein vorhersehbares Verhalten auf Ihrer Plattform liefert (und tatsächlich eine nützliche Beschleunigung gegenüber sqrtf!).

_{1. sqrt(a^b) = (a^b)^0,5 = a^(b/2)}

_{2. Siehe zB https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Exponent_encoding}

Siehe Oliver Charlesworths Erklärung, warum das so ist fast funktioniert. Ich spreche ein Problem an, das in den Kommentaren angesprochen wurde.

Da mehrere Leute darauf hingewiesen haben, dass dies nicht portierbar ist, finden Sie hier einige Möglichkeiten, wie Sie es portierbarer machen oder zumindest den Compiler dazu bringen können, Ihnen mitzuteilen, ob es nicht funktioniert.

Erstens erlaubt Ihnen C++, dies zu überprüfen std::numeric_limits<float>::is_iec559 zur Kompilierzeit, wie in a static_assert. Das kannst du auch überprüfen sizeof(int) == sizeof(float)was nicht wahr ist, wenn int ist 64-Bit, aber was Sie wirklich tun möchten, ist zu verwenden uint32_t, die, wenn sie existiert, immer genau 32 Bit breit ist, ein wohldefiniertes Verhalten bei Verschiebungen und Überläufen hat und einen Kompilierungsfehler verursacht, wenn Ihre seltsame Architektur keinen solchen integralen Typ hat. So oder so, sollten Sie auch static_assert() dass die Typen die gleiche Größe haben. Statische Zusicherungen haben keine Laufzeitkosten und Sie sollten Ihre Vorbedingungen nach Möglichkeit immer auf diese Weise überprüfen.

Leider ist der Test, ob die Konvertierung der Bits in a float zu einem uint32_t und Shifting Big-Endian, Little-Endian oder keines von beiden ist, kann nicht als konstanter Ausdruck zur Kompilierzeit berechnet werden. Hier habe ich die Laufzeitprüfung in den Teil des Codes eingefügt, der davon abhängt, aber Sie möchten sie vielleicht in die Initialisierung einfügen und nur einmal ausführen. In der Praxis können sowohl gcc als auch clang diesen Test zur Kompilierzeit optimieren.

Sie möchten die unsichere Zeigerumwandlung nicht verwenden, und es gibt einige Systeme, an denen ich in der realen Welt gearbeitet habe, bei denen dies das Programm mit einem Busfehler zum Absturz bringen könnte. Die maximal tragbare Methode zum Konvertieren von Objektdarstellungen ist with memcpy(). In meinem Beispiel unten schreibe ich ein Wortspiel mit a union, die auf jeder tatsächlich existierenden Implementierung funktioniert. (Sprachanwälte lehnen dies ab, aber kein erfolgreicher Compiler wird jemals so viel Legacy-Code brechen schweigend.) Wenn Sie eine Zeigerkonvertierung durchführen müssen (siehe unten), gibt es eine alignas(). Aber wie auch immer Sie es tun, das Ergebnis wird implementierungsabhängig sein, weshalb wir das Ergebnis der Konvertierung und Verschiebung eines Testwerts überprüfen.

Wie auch immer, nicht dass Sie es wahrscheinlich auf einer modernen CPU verwenden werden, hier ist eine verspielte C++ 14-Version, die diese nicht-portablen Annahmen überprüft:

#include <cassert>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <vector>

using std::cout;
using std::endl;
using std::size_t;
using std::sqrt;
using std::uint32_t;

template <typename T, typename U>
  inline T reinterpret(const U x)
/* Reinterprets the bits of x as a T.  Cannot be constexpr
 * in C++14 because it reads an inactive union member.
 */
{
  static_assert( sizeof(T)==sizeof(U), "" );
  union tu_pun {
    U u = U();
    T t;
  };
  
  const tu_pun pun{x};
  return pun.t;
}

constexpr float source = -0.1F;
constexpr uint32_t target = 0x5ee66666UL;

const uint32_t after_rshift = reinterpret<uint32_t,float>(source) >> 1U;
const bool is_little_endian = after_rshift == target;

float est_sqrt(const float x)
/* A fast approximation of sqrt(x) that works less well for subnormal numbers.
 */
{
  static_assert( std::numeric_limits<float>::is_iec559, "" );
  assert(is_little_endian); // Could provide alternative big-endian code.
  
 /* The algorithm relies on the bit representation of normal IEEE floats, so
  * a subnormal number as input might be considered a domain error as well?
  */
  if ( std::isless(x, 0.0F) || !std::isfinite(x) )
    return std::numeric_limits<float>::signaling_NaN();
  
  constexpr uint32_t magic_number = 0x1fbb4000UL;
  const uint32_t raw_bits = reinterpret<uint32_t,float>(x);
  const uint32_t rejiggered_bits = (raw_bits >> 1U) + magic_number;
  return reinterpret<float,uint32_t>(rejiggered_bits);
}

int main(void)
{  
  static const std::vector<float> test_values{
    4.0F, 0.01F, 0.0F, 5e20F, 5e-20F, 1.262738e-38F };
  
  for ( const float& x : test_values ) {
    const double gold_standard = sqrt((double)x);
    const double estimate = est_sqrt(x);
    const double error = estimate - gold_standard;
    
    cout << "The error for (" << estimate << " - " << gold_standard << ") is "
         << error;

    if ( gold_standard != 0.0 && std::isfinite(gold_standard) ) {
      const double error_pct = error/gold_standard * 100.0;
      cout << " (" << error_pct << "%).";
    } else
      cout << '.';

    cout << endl;
  }

  return EXIT_SUCCESS;
}

Aktualisieren

Hier ist eine alternative Definition von reinterpret<T,U>() das vermeidet Wortspielereien. Sie könnten das Typ-Wortspiel auch in modernem C implementieren, wo es standardmäßig erlaubt ist, und die Funktion als aufrufen extern "C". Ich denke, Typ-Wortspiel ist eleganter, typsicherer und im Einklang mit dem quasi-funktionalen Stil dieses Programms als memcpy(). Ich glaube auch nicht, dass Sie viel gewinnen, weil Sie immer noch ein undefiniertes Verhalten von einer hypothetischen Fallendarstellung haben könnten. Außerdem ist clang++ 3.9.1 -O -S in der Lage, die Type-Puning-Version statisch zu analysieren und die Variable zu optimieren is_little_endian zur Konstante 0x1und den Laufzeittest eliminieren, aber er kann diese Version nur bis zu einem Stub mit einer einzigen Anweisung optimieren.

Aber was noch wichtiger ist, es ist nicht garantiert, dass dieser Code auf jedem Compiler portabel funktioniert. Beispielsweise können einige alte Computer nicht einmal genau 32 Bit Speicher adressieren. Aber in diesen Fällen sollte es nicht kompilieren und Ihnen sagen, warum. Kein Compiler wird einfach ohne Grund plötzlich eine riesige Menge an Legacy-Code brechen. Obwohl der Standard dies technisch erlaubt und immer noch sagt, dass er C++14 entspricht, wird dies nur auf einer Architektur geschehen, die sich stark von der erwarteten unterscheidet. Und wenn unsere Annahmen so ungültig sind, dass ein Compiler ein Wortspiel zwischen a umwandelt float und eine 32-Bit-Ganzzahl ohne Vorzeichen in einen gefährlichen Fehler, bezweifle ich wirklich, dass die Logik hinter diesem Code bestehen bleibt, wenn wir nur verwenden memcpy() stattdessen. Wir möchten, dass dieser Code zur Kompilierzeit fehlschlägt und uns mitteilt, warum.

#include <cassert>
#include <cstdint>
#include <cstring>

using std::memcpy;
using std::uint32_t;

template <typename T, typename U> inline T reinterpret(const U &x)
/* Reinterprets the bits of x as a T.  Cannot be constexpr
 * in C++14 because it modifies a variable.
 */
{
  static_assert( sizeof(T)==sizeof(U), "" );
  T temp;
  
  memcpy( &temp, &x, sizeof(T) );
  return temp;
}

constexpr float source = -0.1F;
constexpr uint32_t target = 0x5ee66666UL;

const uint32_t after_rshift = reinterpret<uint32_t,float>(source) >> 1U;
extern const bool is_little_endian = after_rshift == target;

Stroustrup et al., in der C++ Core-Richtlinienempfehle a reinterpret_cast stattdessen:

#include <cassert>

template <typename T, typename U> inline T reinterpret(const U x)
/* Reinterprets the bits of x as a T.  Cannot be constexpr
 * in C++14 because it uses reinterpret_cast.
 */
{
  static_assert( sizeof(T)==sizeof(U), "" );
  const U temp alignas(T) alignas(U) = x;
  return *reinterpret_cast<const T*>(&temp);
}

Die von mir getesteten Compiler können dies auch zu einer gefalteten Konstante optimieren. Stroustrups Argumentation ist [sic]:

Zugriff auf das Ergebnis einer reinterpret_cast auf einen anderen Typ als den deklarierten Typ des Objekts zu ändern, ist immer noch ein undefiniertes Verhalten, aber zumindest können wir sehen, dass etwas Verzwicktes vor sich geht.

Aktualisieren

Aus den Kommentaren: C++20 führt ein std::bit_castdie eine Objektdarstellung in einen anderen Typ mit konvertiert nicht spezifiziertnicht nicht definiert, Verhalten. Dies garantiert nicht, dass Ihre Implementierung dasselbe Format verwendet float und int das dieser Code erwartet, aber es gibt dem Compiler keinen Freibrief, Ihr Programm willkürlich zu unterbrechen, weil es technisch undefiniertes Verhalten in einer Zeile davon gibt. Es kann Ihnen auch eine geben constexpr Wandlung.

Sei y = sqrt(x),

aus den Eigenschaften von Logarithmen folgt log(y) = 0,5 * log(x) (1)

Interpretieren eines Normalen float als ganze Zahl ergibt INT(x) = Ix = L * (log(x) + B – σ) (2)

wobei L = 2 ^ N, N die Anzahl der Bits des Signifikanten, B die Exponentenabweichung und σ ein freier Faktor zum Abstimmen der Annäherung ist.

Die Kombination von (1) und (2) ergibt: Iy = 0,5 * (Ix + (L * (B – σ)))

Was im Code so geschrieben ist (*(int*)&x >> 1) + 0x1fbb4000;

Finden Sie das σ so, dass die Konstante gleich 0x1fbb4000 ist, und bestimmen Sie, ob es optimal ist.

Hinzufügen eines Wiki-Testrahmens, um alles zu testen float.

Die Annäherung liegt für viele innerhalb von 4% floataber sehr schlecht für subnormale Zahlen. YMMV

Worst:1.401298e-45 211749.20%
Average:0.63%
Worst:1.262738e-38 3.52%
Average:0.02%

Beachten Sie, dass bei einem Argument von +/-0,0 das Ergebnis nicht Null ist.

printf("% e % e\n", sqrtf(+0.0), sqrt_apx(0.0));  //  0.000000e+00  7.930346e-20
printf("% e % e\n", sqrtf(-0.0), sqrt_apx(-0.0)); // -0.000000e+00 -2.698557e+19

Code testen

#include <float.h>
#include <limits.h>
#include <math.h>
#include <stddef.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>

float sqrt_apx(float f) {
  const int result = 0x1fbb4000 + (*(int*) &f >> 1);
  return *(float*) &result;
}

double error_value = 0.0;
double error_worst = 0.0;
double error_sum = 0.0;
unsigned long error_count = 0;

void sqrt_test(float f) {
  if (f == 0) return;
  volatile float y0 = sqrtf(f);
  volatile float y1 = sqrt_apx(f);
  double error = (1.0 * y1 - y0) / y0;
  error = fabs(error);
  if (error > error_worst) {
    error_worst = error;
    error_value = f;
  }
  error_sum += error;
  error_count++;
}

void sqrt_tests(float f0, float f1) {
  error_value = error_worst = error_sum = 0.0;
  error_count = 0;
  for (;;) {
    sqrt_test(f0);
    if (f0 == f1) break;
    f0 = nextafterf(f0, f1);
  }
  printf("Worst:%e %.2f%%\n", error_value, error_worst*100.0);
  printf("Average:%.2f%%\n", error_sum / error_count);
  fflush(stdout);
}

int main() {
  sqrt_tests(FLT_TRUE_MIN, FLT_MIN);
  sqrt_tests(FLT_MIN, FLT_MAX);
  return 0;
}