double을 사용하는 정수 sqrt의 정확도

Question

uint64_t의 정수 부분을 계산하고 싶습니다. 32 비트 uint32_t의 경우 먼저 double, sqrt으로 캐스팅 한 다음 uint32_t으로 다시 캐스팅하는 것이 좋습니다.

double은 정확히 2^53까지만 숫자를 수용 할 수 있다고 가정하면 uint64_t에서도 작동합니까? 즉, 항상 다음과 같은 올바른 대답을 줄 것입니다 : 심지어

#include <math.h> 
uint64_t x = ...; 
uint64_t result = (uint64_t)sqrt((double)x); 

또는 :

#include <math.h> 
uint64_t x = ...; 
uint32_t result = (uint32_t)sqrt((double)x); 

Answer 1

경험적으로, 대답은 더 없습니다. 입력 4503599761588224에 대한 결과는 67108864 대신 67108865로 잘못 계산됩니다.

다음 코드는이 사례를 식별합니다. 물론 break;을 제거하여 다른 경우를 관찰 할 수 있습니다.

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

int main(void) { 
    for (uint32_t y = 1; y != 0; y++) { 
     // *Just* smaller than a perfect square 
     uint64_t x = ((uint64_t)y * (uint64_t)y) - 1; 

     // We expect the floor of the result  
     uint32_t expected = y - 1; 

     uint32_t result = (uint32_t)sqrt((double)x); 

     if (result != expected) { 
      printf("Incorrect: x = %llu, result = %u\n", x, result); 
      break; 
     } 
    } 
    return 0; 
} 

값 4503599761588224의 특별한 점은 무엇입니까? 음, 정확하게 (2 +1) -1, AKA (2)입니다. 정확하게 double으로 표시 될 수 있으므로 long ->double 변환으로 인한 오류가 아닙니다.

대신 오류는 sqrt 구현 내부에 있습니다. 여기서 델타 (대 완벽한 사각형)는 제곱근을 약 2 ^-27으로 줄입니다.이 값은 2 번인 result 자체입니다. 이것은 배정도가 처리 할 수있는 한계에 달려 있으므로 당연히 오프 - 바이 - 원 오류를 볼 것으로 예상됩니다.

---

_{1. Live demo 2}

. 당신이 당신의 수학 라이브러리의`sqrt` 좋은 및 C 알고있는 경우에도`uint32_t`에 대한 :)