整数の精度sqrt doubleを使用

Question

uint64_tの整数部分を計算したいと思います。 32ビットのuint32_tの場合は、最初にdouble、sqrtにキャストしてからuint32_tにキャストすることをお勧めします。

は、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 + 2 ）です。これは正確にdoubleで表すことができるので、エラーはlong - >doubleの変換によるものではありません。

代わりに、エラーはsqrt実装の内部にあります。ここでのデルタ（完全な正方形に対する）は、平方根を約2×^-27だけ減少させ、それは約2 倍であり、それ自体resultより小さい。これは、倍精度で扱うことができる限界にあります。この時点では、当然、オフ・バイ・ワンのエラーが発生することが予想されます。

---

_{1. Live demo。あなたは数学ライブラリの `sqrt`が良いですし、あなたのCを知っている場合でも、` uint32_t`ため:)}