Comment fonctionne cette approximation de la racine carrée flottante ?

Question

J'ai trouvé une approximation de racine carrée plutôt étrange mais qui fonctionne pour float s ; je ne comprends vraiment pas. Quelqu'un peut-il m'expliquer pourquoi ce code fonctionne ?

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

Je l'ai testé un peu et il produit des valeurs à partir de std::sqrt() d'environ 1 à 3 %. . Je connais les Quake III Racine carrée inverse rapide et je suppose que c'est quelque chose de similaire ici (sans l'itération de Newton) mais j'apprécierais vraiment une explication de comment cela fonctionne .

(nota : Je l'ai étiqueté à la fois c y c++ puisqu'il s'agit d'un code C et C++ à la fois valide (voir les commentaires).

Answer 1

(*(int*)&f >> 1) décale vers la droite la représentation binaire de f . Ce site presque divise l'exposant par deux, ce qui équivaut approximativement à prendre la racine carrée. 1

Pourquoi presque ? Dans IEEE-754, l'exposant réel est e - 127 . 2 Pour diviser ce chiffre par deux, il faudrait e/2 - 64 mais l'approximation ci-dessus ne nous donne que e/2 - 127 . Nous devons donc ajouter 63 à l'exposant résultant. Cet ajout est apporté par les bits 30 à 23 de cette constante magique ( 0x1fbb4000 ).

J'imagine que les bits restants de la constante magique ont été choisis pour minimiser l'erreur maximale sur la plage de la mantisse, ou quelque chose comme ça. Cependant, il n'est pas clair si cela a été déterminé analytiquement, itérativement ou heuristiquement.

---

Il convient de souligner que cette approche est quelque peu non portable. Elle fait (au moins) les hypothèses suivantes :

• La plate-forme utilise la précision simple IEEE-754 pour le calcul de l'impôt sur le revenu. float .
• L'endiannité de float représentation.
• que vous ne serez pas affecté par un comportement non défini, car cette approche viole les règles du C/C++. règles strictes d'anticrénelage .

Ainsi, il devrait être évité à moins que vous ne soyez certain qu'il donne un comportement prévisible sur votre plate-forme (et en fait, qu'il fournit un gain de vitesse utile par rapport à l'autre. sqrtf !).

---

_{1. <em>sqrt(a^b) = (a^b)^0.5 = a^(b/2)</em>}

_{2. Voir par exemple <a href="https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Exponent_encoding" rel="noreferrer">https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Exponent_encoding</a>}

Answer 2

Voir l'explication d'Oliver Charlesworth sur la raison pour laquelle ceci presque fonctionne. Je réponds à une question soulevée dans les commentaires.

Comme plusieurs personnes ont souligné la non-portabilité de ce système, voici quelques moyens de le rendre plus portable, ou au moins de faire en sorte que le compilateur vous indique si cela ne fonctionne pas.

Tout d'abord, le C++ vous permet de vérifier std::numeric_limits<float>::is_iec559 au moment de la compilation, comme dans un static_assert . Vous pouvez également vérifier que sizeof(int) == sizeof(float) ce qui ne sera pas vrai si int est de 64 bits, mais ce que vous voulez vraiment faire c'est utiliser uint32_t qui, s'il existe, sera toujours d'une largeur de 32 bits exactement, aura un comportement bien défini avec les décalages et les débordements, et provoquera une erreur de compilation si votre architecture bizarre ne possède pas un tel type d'intégrale. Dans tous les cas, vous devriez également static_assert() que les types ont la même taille. Les assertions statiques n'ont aucun coût d'exécution et vous devriez toujours vérifier vos préconditions de cette façon si possible.

Malheureusement, le test pour savoir si la conversion des bits d'une float à un uint32_t et le décalage est big-endian, little-endian ou ni l'un ni l'autre ne peut pas être calculé comme une expression constante du temps de compilation. Ici, j'ai mis la vérification au moment de l'exécution dans la partie du code qui en dépend, mais vous pourriez vouloir la mettre dans l'initialisation et la faire une fois. En pratique, gcc et clang peuvent tous deux optimiser ce test au moment de la compilation.

Vous ne voulez pas utiliser le "unsafe pointer cast", et il y a certains systèmes sur lesquels j'ai travaillé dans le monde réel où cela pourrait faire planter le programme avec une erreur de bus. La façon la plus portable de convertir les représentations d'objets est avec memcpy() . Dans mon exemple ci-dessous, je tape-pun avec un union qui fonctionne sur toute implémentation existante. (Les avocats spécialisés dans les langues s'y opposent, mais aucun compilateur performant ne cassera jamais autant de code existant. en silence .) Si vous devez faire une conversion de pointeur (voir ci-dessous), il y a alignas() . Mais quelle que soit la façon dont vous le faites, le résultat sera défini par l'implémentation, c'est pourquoi nous vérifions le résultat de la conversion et du décalage d'une valeur de test.

Quoi qu'il en soit, bien que vous ne soyez pas susceptible de l'utiliser sur un processeur moderne, voici une version C++14 améliorée qui vérifie ces hypothèses non portables :

#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;
}

Mise à jour

Voici une autre définition de reinterpret<T,U>() qui évite le détournement de type. Vous pourriez également implémenter le détournement de type en C moderne, où il est autorisé par la norme, et appeler la fonction comme suit extern "C" . Je pense que la ruse de type est plus élégante, plus sûre et plus cohérente avec le style quasi-fonctionnel de ce programme que la ruse de type. memcpy() . Je ne pense pas non plus que vous gagniez beaucoup, car vous pourriez toujours avoir un comportement indéfini à partir d'une hypothétique représentation piège. Par ailleurs, clang++ 3.9.1 -O -S est capable d'analyser statiquement la version de type-punning, d'optimiser la variable is_little_endian à la constante 0x1 et éliminer le test d'exécution, mais il ne peut optimiser cette version que jusqu'à un stub à instruction unique.

Mais plus important encore, ce code n'est pas garanti pour fonctionner de manière portable sur tous les compilateurs. Par exemple, certains vieux ordinateurs ne peuvent même pas adresser exactement 32 bits de mémoire. Mais dans ces cas-là, il devrait échouer à la compilation et vous dire pourquoi. Aucun compilateur ne va soudainement casser une énorme quantité de code hérité sans raison. Bien que la norme donne techniquement la permission de le faire et de continuer à dire qu'il est conforme à C++14, cela ne se produira que sur une architecture très différente de celle à laquelle nous nous attendons. Et si nos hypothèses sont si invalides qu'un compilateur va transformer un type-pun entre un float et un entier non signé de 32 bits en un bug dangereux, je doute vraiment que la logique derrière ce code tienne si nous utilisons simplement memcpy() à la place. Nous voulons que ce code échoue à la compilation et qu'il nous dise pourquoi.

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

Toutefois, Stroustrup et al. Directives de base du C++ recommande un reinterpret_cast à la place :

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

Les compilateurs que j'ai testés peuvent également optimiser cette fonction en la réduisant à une constante repliée. Le raisonnement de Stroustrup est [sic] :

Accéder au résultat d'un reinterpret_cast à un type différent du type déclaré de l'objet est toujours un comportement non défini, mais au moins nous pouvons voir que quelque chose de délicat se passe.

Answer 3

Soit y = sqrt(x),

il découle des propriétés des logarithmes que log(y) = 0,5 * log(x) (1)

Interprétation d'une normale float en tant que nombre entier donne INT(x) = Ix = L * (log(x) + B - ) (2)

où L = 2^N, N le nombre de bits du significande, B est le biais de l'exposant, et est un facteur libre pour ajuster l'approximation.

La combinaison de (1) et (2) donne : Iy = 0,5 * (Ix + (L * (B - )))

Ce qui est écrit dans le code comme (*(int*)&x >> 1) + 0x1fbb4000;

Trouvez la solution pour que la constante soit égale à 0x1fbb4000 et déterminez si elle est optimale.

Answer 4

Ajout d'un harnais de test wiki pour tester toutes les float .

L'approximation est à 4% près pour de nombreux float mais très mauvais pour les nombres inférieurs à la normale. YMMV

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

Notez qu'avec un argument de +/-0.0, le résultat n'est pas nul.

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 de test

#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;
}