Как работает это сглаживание с квадратным корнем?

Ответ 1

(*(int*)&f >> 1) сдвигает вправо поразрядное представление f. Это почти делит показатель на два, что приблизительно эквивалентно взятию квадратного корня. ¹

Почему почти? В IEEE-754 фактический показатель равен e-127. ² Чтобы разделить это на два, нам понадобится e/2 - 64, но приведенное выше приближение дает только e/2 - 127. Поэтому нам нужно добавить 63 к полученному экспоненту. Это обеспечивается битами 30-23 этой магической константы (0x1fbb4000).

Я бы предположил, что оставшиеся биты волшебной константы были выбраны для минимизации максимальной ошибки в диапазоне мантиссы или что-то в этом роде. Однако неясно, было ли оно определено аналитически, итеративно или эвристически.

Стоит отметить, что этот подход несколько не переносимый. Это делает (по крайней мере) следующие предположения:

• На платформе используется IEEE-754 с одной точностью для float.
• Консистенция представления float.
• Чтобы вы не пострадали от поведения undefined из-за того, что этот подход нарушает правила строгого сглаживания C/С++ .

Таким образом, этого следует избегать, если вы не уверены, что он дает предсказуемое поведение на вашей платформе (и действительно, что он обеспечивает полезное ускорение и sqrtf!).

<суб > 1. sqrt (a ^ b) = (a ^ b) ^ 0,5 = a ^ (b/2)

<суб > 2. См. https://en.wikipedia.org/wiki/Single-precision_floating-point_format#Exponent_encoding

Ответ 2

См. описание Оливера Чарльворта, почему это почти работает. Я рассматриваю проблему, поднятую в комментариях.

Поскольку несколько человек указали на непереносимость этого, вот несколько способов сделать его более переносимым или, по крайней мере, заставить компилятор сказать вам, не будет ли это работать.

Во-первых, С++ позволяет вам проверять std::numeric_limits<float>::is_iec559 во время компиляции, например, в static_assert. Вы также можете проверить, что sizeof(int) == sizeof(float), что не будет истинным, если int - 64-разрядный, но то, что вы действительно хотите сделать, это использовать uint32_t, который, если он существует, всегда будет ровно 32 бит в ширину, будет иметь четко определенное поведение со сдвигами и переполнением, и вызовет ошибку компиляции, если ваша странная архитектура не имеет такого интегрального типа. В любом случае, вы также должны static_assert(), чтобы типы имели одинаковый размер. Статические утверждения не имеют затрат времени исполнения, и вы всегда должны проверять свои предпосылки таким образом, если это возможно.

К сожалению, проверка того, преобразует ли биты в float в uint32_t и смещение, является big-endian, little-endian или никоим образом не может быть вычислен как выражение постоянной времени компиляции. Здесь я поставил проверку времени выполнения части кода, которая зависит от нее, но вы можете поместить ее в инициализацию и сделать это один раз. На практике, как gcc, так и clang могут оптимизировать этот тест во время компиляции.

Вы не хотите использовать небезопасный указатель, и есть некоторые системы, над которыми я работал в реальном мире, где это может привести к сбою программы с ошибкой шины. Максимально переносимый способ преобразования представлений объектов - memcpy(). В моем примере ниже я пишу каламбуром с помощью union, который работает с любой фактически существующей реализацией. (Юридические юристы возражают против этого, но ни один успешный компилятор никогда не сломает этот многозначный код молча.) Если вы должны сделать преобразование указателя (см. Ниже), есть alignas(). Но как бы вы это ни делали, результат будет определяться реализацией, поэтому мы проверяем результат преобразования и изменения тестового значения.

Во всяком случае, не то, что вы, скорее всего, будете использовать его на современном процессоре, есть русифицированная версия С++ 14, которая проверяет эти непереносимые предположения:

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

Update

Ниже приведено альтернативное определение reinterpret<T,U>(), которое позволяет избежать путаницы. Вы также можете реализовать тип-каламбур в современном C, где его разрешено по стандарту, и вызывать функцию как extern "C". Я думаю, что тип-punning более элегантен, безопасен по типу и соответствует квазифункциональному стилю этой программы, чем memcpy(). Я также не думаю, что вы набираете много, потому что у вас все еще может быть поведение undefined из гипотетического представления ловушки. Кроме того, clang++ 3.9.1 -O -S способен статически анализировать версию типа Punning, оптимизировать переменную is_little_endian до константы 0x1 и исключить тест времени выполнения, но она может оптимизировать эту версию только к одиночной инструкции.

Но что более важно, этот код не гарантированно работает на каждом компиляторе. Например, некоторые старые компьютеры не могут даже адресовать ровно 32 бита памяти. Но в этих случаях он не должен компилироваться и говорить вам почему. Ни один компилятор просто не собирается разрушать огромное количество устаревшего кода. Хотя стандарт технически дает разрешение на это и все еще говорит, что он соответствует С++ 14, это произойдет только в архитектуре, очень отличающейся от ожидаемой. И если наши допущения настолько недействительны, что какой-то компилятор собирается превратить ключевое слово между float и 32-разрядным целым без знака в опасную ошибку, я действительно сомневаюсь, что логика этого кода задержит, если мы просто используем memcpy(). Мы хотим, чтобы этот код не работал во время компиляции, и чтобы сообщить нам, почему.

#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 и др., в Основные принципы С++, рекомендуется вместо reinterpret_cast:

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

Компиляторы, которые я тестировал, также могут оптимизировать это до сложенной константы. Примером Страустрапа является [sic]:

Доступ к результату reinterpret_cast к другому типу из объявленного типа объектов по-прежнему выполняется undefined, но по крайней мере мы можем видеть, что происходит что-то сложное.

Ответ 3

Пусть y = sqrt (x),

из свойств логарифмов следует, что log (y) = 0,5 * log (x) (1)

Интерпретация нормального float как целого дает INT (x) = Ix = L * (log (x) + B - σ) (2)

где L = 2 ^ N, N - число бит значащего, B - смещение экспоненты, а σ - свободный коэффициент для настройки аппроксимации.

Сочетание (1) и (2) дает: Iy = 0,5 * (Ix + (L * (B - σ)))

Что написано в коде как (*(int*)&x >> 1) + 0x1fbb4000;

Найдите σ так, чтобы константа равнялась 0x1fbb4000 и определяла, является ли она оптимальной.

Ответ 4

Добавление тестового жгута для проверки всех float.

Аппроксимация составляет не более 4% для многих float, но очень низкая для под нормальных чисел. YMMV

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

Обратите внимание, что с аргументом +/- 0.0 результат не равен нулю.

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

Тестовый код

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