지진 역 제곱근 : 정확도

Question

분명히 퀘이크 게임을 거슬러 올라가는 역 제곱근을 계산하는 "마법"방법은 많은 출처에서 설명됩니다. 위키 백과는 그것에 좋은 기사가 있습니다 https://cs.uwaterloo.ca/~m32rober/rsqrt.pdf

나는이 논문에서 이러한 결과의 일부를 복제하기 위해 노력하고있어 : 내가 특히 발견 https://en.wikipedia.org/wiki/Fast_inverse_square_root

을 알고리즘의 아주 좋은 쓰기 업 및 분석을 할 다음 ,하지만 정확성에 문제가 있습니다. 다음 C 코딩 알고리즘은, 주어진다

#include <math.h> 
#include <stdio.h> 

float Q_rsqrt(float number) { 
    long i; 
    float x2, y; 
    const float threehalfs = 1.5F; 

    x2 = number * 0.5F; 
    y = number; 
    i = *(long *) &y; 
    i = 0x5f3759df - (i >> 1); 
    y = *(float *) &i; 
    y = y * (threehalfs - (x2 * y * y)); 
    // y = y * (threehalfs - (x2 * y * y)); 
    return y; 
} 

paper 상태 상대 오차가 모두 양의 정상 수레 0.0017522874 많아야된다. (. 코드, 및 1.4 절에서 논의 부록 2 참조)

그러나, 나는 "플러그인"때 숫자 1.4569335e-2F, 내가 얻을 오류가이보다 큰 예측 오차 :

int main() 
{ 

    float f = 1.4569335e-2F; 

    double tolerance = 0.0017522874; 
    double actual = 1.0/sqrt(f); 
    float magic  = Q_rsqrt(f); 
    double err  = fabs (sqrt(f) * (double) magic - 1); 

    printf("Input : %a\n", f); 
    printf("Actual : %a\n", actual); 
    printf("Magic : %a\n", magic); 
    printf("Err  : %a\n", err); 
    printf("Tolerance: %a\n", tolerance); 
    printf("Passes : %d\n", err <= tolerance); 

    return 0; 
} 

출력은 :

Input : 0x1.dd687p-7 
Actual : 0x1.091cc953ea828p+3 
Magic : 0x1.08a5dcp+3 
Err  : 0x1.cb5b716b7b6p-10 
Tolerance: 0x1.cb5a044e0581p-10 
Passes : 0 

그래서, 이러한 특정 입력은 종이로 이루어지는 제 위배 보인다.

이것이 용지 자체에 문제가 있는지 또는 코딩 과정에서 실수를 한 것인지 궁금합니다. 나는 어떤 의견을 주셔서 감사합니다!

Answer 1

이의이 (가) 상대 오차에 바인딩과는 thesis of Matthew Robertson에서보다 약간 더 큰하다고 보여 재 계산하는 코드의 작은 조각을 해보자 :

여기에 내가 0x5f375a86을 사용하여 얻은 결과이다. 실제로, @skamamishossifrage의 대답에서 처음 발견되었고 the thesis of Matthew Robertson에 기록 된 바와 같이이 구현은 Quake III 소스에서 공개 된 것입니다. 특히 Quake III의 상수 값은 Quake III 소스에서 파일 번호 q_math.c의 561 행에 있습니다.

먼저 64 비트 플레이트 폼에서 작동하도록 코드를 수정해야합니다. 수정해야 할 수있는 유일한 유형은 정수 유형입니다. long은 플레이트 형태와 관련이 없습니다. 제 리눅스 컴퓨터에서 sizeof(long)은 8을 반환합니다 ... 49 페이지의 논문에서 업데이트 된 것처럼 uint32_t 유형은 float과 같은 크기의 정수 유형을 보장합니다. 여기

코드, gcc main.c -o main -lm -Wall에 의해 컴파일 ./main에 의해 실행되는 이동 : 바운드를 들어

#include <math.h> 
#include <stdio.h> 
#include <inttypes.h> 

float Q_rsqrt(float number) { 
    uint32_t i; 
    float x2, y; 
    const float threehalfs = 1.5F; 

    x2 = number * 0.5F; 
    y = number; 
    i = *(uint32_t *) &y; 
    i = 0x5f3759df - (i >> 1); // 0x5f3759df 0x5f375a86 
    y = *(float *) &i; 
    y = y * (threehalfs - (x2 * y * y)); 
    // y = y * (threehalfs - (x2 * y * y)); 
    return y; 
} 

int main() 
{ 

    printf("%ld %ld\n",sizeof(long),sizeof(uint32_t)); 

    uint32_t i; 
    float y; 
    double e, max = 0.0; 
    float maxval=0; 
    for(i = 0x0000000; i < 0x6f800000; i++) { 
     y = *(float *) &i; 
     if(y>1e-30){ 
      e = fabs(sqrt((double)y)*(double)Q_rsqrt(y) - 1); 
      if(e > max){ 
       max = e; 
       maxval=y; 
      } 
     } 
    } 
    printf("On value %2.8g == %a\n", maxval, maxval); 
    printf("The bound is %2.12g == %a\n", max, max); 

    return 0; 
} 

을, 나는 0.0017523386721 == 0x1.cb5d752717ep-10를 얻을. 당신이 주목했듯이, 그것은 신문 (0.001752287)에서보고 된 것보다 약간 큽니다. double 대신 float을 사용하여 오류를 평가해도 결과가 크게 바뀌지 않습니다.

Answer 2

잘못된 마법 번호를 사용 중입니다.

0x5f3759df은 원래 Quake III에서 사용 된 값이지만 나중에 0x5f375a86이 더 나은 결과를 제공한다는 것을 알게되었습니다. 인용 한 논문의 40 페이지에있는 그림 6.1을 보면 향상된 상수를 사용하고 있음을 알 수 있습니다.

Input : 0x1.dd687p-7 
Actual : 0x1.091cc953ea828p+3 
Magic : 0x1.08a5fap+3 
Err  : 0x1.cae79153f2cp-10 
Tolerance: 0x1.cb5a044e0581p-10 
Passes : 1