이 플로트 제곱근 근사치는 어떻게 작동합니까?

이 플로트 제곱근 근사치는 어떻게 작동합니까?

하고 있는 요.floats; 정말 이해가 안 돼요.누가 이 코드가 왜 작동하는지 설명해 줄 수 있나요?

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

테스트를 좀 해봤더니 1~3% 정도 값이 출력됩니다.나는 퀘이크 III의 빠른 역제곱근에 대해 알고 있고, 나는 이것이 여기서 비슷한 것이라고 추측한다(뉴턴 반복 없이). 하지만 나는 그것이 어떻게 작동하는지 설명해주면 정말 좋겠다.

(주의: c와 c++ 둘 다 유효하기 때문에 태그를 붙였습니다.(댓글 참조) C와 C++ 코드)

(*(int*)&f >> 1)의 f이것은 지수를 2로 거의 나누는데,¹ 이는 제곱근을 구하는 것과 거의 같습니다.

왜 거의?IEEE-754에서는 실제 지수는 e - ²127입니다.이를 2로 나누려면 e/2 - 64가 필요한데 위의 근사치에서는 e/2 - 127만 나옵니다.따라서 63을 더해야 합니다.이는 해당 매직 상수의 비트 30-23에 의해 발생합니다(0x1fbb4000를 참조해 주세요.

마법 상수의 나머지 비트는 가수 범위에 걸쳐 최대 오차를 최소화하기 위해 선택되었을 것입니다.그러나, 그것이 분석적으로 결정되었는지, 반복적으로 결정되었는지, 또는 휴리스틱하게 결정되었는지 불분명하다.

---

이 접근방식은 다소 휴대성이 없다는 것을 지적할 필요가 있습니다.적어도 다음과 같은 가정을 합니다.

• 에서는, 「IEEE-754」에 를 사용합니다.float.
• ★★★★★★★★의 엔디안성float★★★★★★ 。
• 이 접근방식은 C/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 등 참조}

Oliver Charlesworth의 설명을 참조하십시오.댓글로 제기된 문제를 다루고 있습니다.

여러 사람이 이 기능을 휴대할 수 없다고 지적했기 때문에 이 기능을 보다 휴대하기 쉽게 하거나 최소한 컴파일러가 동작하지 않을 경우 알려 주도록 할 수 있는 몇 가지 방법이 있습니다.

하면 C++를 체크할 수 .std::numeric_limits<float>::is_iec559를 컴파일 시static_assert 것도 할 수 sizeof(int) == sizeof(float) 경우는 되지 않습니다.int이지만, 하고 싶은 은 64비트입니다.uint32_t이는 항상 정확히 32비트 너비이며, 시프트와 오버플로우가 있는 명확한 동작을 하며, 이상한 아키텍처에 이러한 통합 유형이 없는 경우 컴파일 오류가 발생합니다.쪽이든, 은 '', '아예', '아예', '아예'도 해야 합니다.static_assert()두 종류가 같은 크기라는 거죠.정적 어설션에는 런타임 비용이 들지 않으므로 가능하면 항상 이 방법으로 사전 조건을 확인해야 합니다.

하고 있는지 입니다.float a까지uint32_t이동은 빅 엔디안, 리틀 엔디안 또는 둘 다 컴파일 시간 상수식으로 계산할 수 없습니다.여기서는 런타임 체크를 코드 부분에 넣지만 초기화에 넣어 한 번 실행하는 것이 좋을지도 모릅니다.실제로 gcc와 clang 모두 컴파일 시 이 테스트를 최적화할 수 있습니다.

안전하지 않은 포인터 캐스트를 사용하고 싶지 않을 것입니다.또한 실제 환경에서 작업한 시스템 중에는 버스 오류로 인해 프로그램이 크래시 될 수 있는 시스템도 있습니다.하기 위한 은 " " " 입니다.memcpy()다음 예에서는 type-pun과union(언어 변호사들은 반대하지만 성공한 컴파일러는 그 많은 레거시 코드를 소리 없이 깨지는 일은 없을 것입니다.)포인터 변환을 실시할 필요가 있는 경우(아래 참조),alignas()그러나 어떤 방법으로 실행해도 결과는 구현 정의되기 때문에 테스트 값의 변환 및 이동 결과를 확인합니다.

어쨌든 최신 CPU에서 사용할 가능성이 높은 것은 아니지만, 다음은 이러한 포터블 이외의 전제 조건을 체크하는 Gused Up C++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;
}

갱신하다

또 .reinterpret<T,U>()을 사용하다 표준으로를 type-pun으로 도 있습니다.extern "C" than of, 이 . . . . . ., . . . . . . . . . . . . . . . . . . . . . . . . . . . . .보다 이 프로그램의 에 더하며 일관성이 memcpy()그리고 아직 가상의 함정 표현에서 행동을 정의하지 못했을 수도 있기 때문에 얻는 것이 많지 않다고 생각합니다.9.1 는, 타입 해, 「clang++ 3.9.1 -O -S」를 할 수 .is_little_endian0x1실행 시 테스트를 배제합니다.단, 이 버전을 최적화할 수 있는 것은 단일 명령 스터브뿐입니다.

그러나 더 중요한 것은 이 코드가 모든 컴파일러에서 휴대할 수 있도록 보장되지 않는다는 것입니다.예를 들어, 일부 오래된 컴퓨터는 정확히 32비트의 메모리를 주소 지정할 수도 없습니다.그러나 이 경우 컴파일에 실패하여 이유를 알 수 없습니다.어떤 컴파일러도 이유 없이 엄청난 양의 레거시 코드를 갑자기 망가뜨릴 수는 없습니다.이 규격은 기술적으로 이를 허용하는 한편 C++14에 준거하고 있다고 말하고 있지만, 이는 우리가 예상하는 것과는 매우 다른 아키텍처에서만 발생합니다.또한 어떤 을 type-pun으로 하려고 할 이 유효하지 float 정수를 버그에 이할 수 입니다.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 등에서는 C++ 핵심 가이드라인에서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);
}

테스트한 컴파일러도 이 값을 접힌 상수로 최적화할 수 있습니다.Strustrup의 논리는 다음과 같습니다.

「 」의 결과에의 reinterpret_cast오브젝트와는 다른 타입으로 선언된 타입은 아직 정의되지 않은 동작이지만 적어도 뭔가 교묘한 일이 일어나고 있는 것을 알 수 있습니다.

갱신하다

코멘트로부터:C++20은 오브젝트 표현을 미지정 동작(정의되지 않은 동작)으로 다른 타입으로 변환합니다.따라서 구현에서 동일한 형식을 사용하는 것은 아닙니다.float ★★★★★★★★★★★★★★★★★」int이 코드는 예상하지만, 컴파일러가 임의로 프로그램을 중단시킬 수 있는 것은 아닙니다.그 한 줄에 기술적으로 정의되지 않은 동작이 있기 때문입니다., ,, 능, 능, 이, 이, constexpr★★★★★★ 。

y = sqrt(x),

log(y) = 0.5 * log(x) (1)인 로그의 속성에서 비롯된다.

의 「」의 float는 INT=* ( + - (2)의 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이 되도록 so를 찾아 최적 여부를 판단합니다.

Wiki를 Wiki 추가float.

4%입니다.float단, 정상 이하의 수치에는 매우 좋지 않습니다.YMMV

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

인수가 +/-0.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;
}

언급URL : https://stackoverflow.com/questions/43120045/how-does-this-float-square-root-approximation-work