整数sqrt的准确度 [英] accuracy of sqrt of integers

问题描述

  for（uint64_t i = 0; i * i< n; i ++）{ 
  

这需要在每次迭代中进行乘法运算。如果我可以在循环之前计算sqrt，那么我可以避免这种情况。

$ $ $ $ $ $ $ $> unsigned cut = sqrt（n）
for（uint64_t i = 0; i< cut; i ++）{

sqrt函数向上舍入到下一个整数，但是如果它舍入，它不是好的。

我的问题是：sqrt函数是否足够准确，可以用于所有情况？

编辑：让我列出一些情况。如果n是一个完美的正方形，那么我的问题就是 cut = sqrt（n）> = y 对于所有的n？如果cut = y-1那么就有问题了。例如。如果n = 120，cut = 10，那么没关系，但是如果n = 121（11 ^ 2），那么cut仍然是10，那么它就不起作用了。



是浮点数的小数部分只有23位和双52，所以他们不能存储一些32位或64位整数的所有数字。不过，我不认为这是一个问题。假设我们想要某个数字y的sqrt，但是我们不能存储y的所有数字。如果我们让y的分数可以存储为x，我们可以写y = x + dx，那么我们要确保我们选择的任何dx都不会将我们移动到下一个整数。



  sqrt（x + dx）< sqrt（x）+ 1 //解
 dx< 2 * sqrt（x）+ 1 
 //例如对于x = 100dx< 21 
 // sqrt（100 + 20）< sqrt（100）+ 1 
  

浮点数可以存储23位，所以我们让y = 2 ^ 23 + 2 ^ 9。这是绰绰有余的，因为2 ^ 9 < 2 * sqrt（2 ^ 23）+ 1.使用64位整数也很容易显示双精度值。所以尽管只要他们可以存储的sqrt是准确的，那么他们就不能存储所有的数字，那么sqrt（分数）就足够了。现在让我们来看看接近INT_MAX和sqrt的整数：
$ b

  unsigned xi = -1-1; 
 printf（％u％u \ n，xi，（unsigned）（float）xi）; // 4294967294 4294967295 
 printf（％u％u \ n，（unsigned）sqrt（xi），（unsigned）sqrtf（xi））; // 65535 65536 
  

由于float不能存储2 ^ 31-2的所有数字，并且double他们可以得到不同的结果sqrt。但是sqrt的float版本是一个更大的整数。这就是我要的。对于64位整数，只要double的sqrt总是舍入就可以了。


解决方案

首先，整数乘法是真的挺便宜的。只要你每个循环迭代有多个工作周期和一个备用执行槽，它应该完全被大多数非微小处理器上的重新排序隐藏。





如果你确实有一个处理器的速度慢得多，那么一个真正聪明的编译器可能会把你的循环转换为：
$ b $ pre code> for（uint64_t i = 0，j = 0; j

<用$ lea 替换乘数，或者用shift和两个加法替换乘数。

$ hr
$ b $撇开这些说明，让我们看看你说的问题。不，你不能只用 i< SQRT（n）的。反例： n = 0x20000000000000 。假设遵守IEEE-754，你将有 cut = 0x5a82799 ，并且 cut * cut 是 0x1ffffff8eff971 。

然而，基本的浮点错误分析表明计算 sqrt（n） （在转换为整数之前）以ULP的3/4为界。所以你可以安全地使用：

$ $ $ $ $ c $ uint32_t cut = sqrt（n）+ 1;


，您最多可以执行一次额外的循环迭代，这可能是可以接受的。如果你想完全精确的话，可以使用：

  uint32_t cut = sqrt（n）; 
 cut + =（uint64_t）cut * cut< N; 
  

---

编辑： z boson 澄清，为了他的目的，这只有当 n 是一个确切的正方形（否则，得到 cut的值即太小是可以接受的）。在这种情况下，不需要调整，就可以安全地使用：

  uint32_t cut = sqrt（n）; 
  

为什么这是真的？实际上，这很简单。将 n 转换为 double 引入扰动：

  double_n = n *（1 + e）
  

满足 | E | < 2 ^ -53 。这个值的数学平方根可以扩展如下：

  square_root（double_n）= square_root（n）* square_root（1 + e）
  

现在，既然 n 是假定是一个最多64位的完美平方， square_root（n）是一个精确的整数，最多有32位，是我们希望计算的数学精确值。要分析 square_root（1 + e）项，请使用关于 1 的泰勒级数：

$ p code $ square_root（1 + e）= 1 + e / 2 + O（e ^ 2）
= 1 + d with | d | <〜2 ^ -54

因此，数学确切值 square_root double_n）小于ULP远离[1]所需的确切答案的一半，并且必然舍入到该值。在我滥用相对误差估计的情况下，我正在快速和松散，ULP的相对大小实际上在不同情况下有所不同 - 我试图给出一些证据的味道，而不会陷入细节。这一切都可以做得非常严格，它只是堆栈溢出有点罗嗦。

I have a loop like this:

for(uint64_t i=0; i*i<n; i++) {

This requires doing a multiplication every iteration. If I could calculate the sqrt before the loop then I could avoid this.

unsigned cut = sqrt(n)
for(uint64_t i=0; i<cut; i++) {

In my case it's okay if the sqrt function rounds up to the next integer but it's not okay if it rounds down.

My question is: is the sqrt function accurate enough to do this for all cases?

Edit: Let me list some cases. If n is a perfect square so that n = y^2 my question would be - is cut=sqrt(n)>=y for all n? If cut=y-1 then there is a problem. E.g. if n = 120 and cut = 10 it's okay but if n=121 (11^2) and cut is still 10 then it won't work.

My first concern was the fractional part of float only has 23 bits and double 52 so they can't store all the digits of some 32-bit or 64-bit integers. However, I don't think this is a problem. Let's assume we want the sqrt of some number y but we can't store all the digits of y. If we let the fraction of y we can store be x we can write y = x + dx then we want to make sure that whatever dx we choose does not move us to the next integer.

sqrt(x+dx) < sqrt(x) + 1  //solve
dx < 2*sqrt(x) + 1 
// e.g for x = 100 dx < 21
// sqrt(100+20) < sqrt(100) + 1

Float can store 23 bits so we let y = 2^23 + 2^9. This is more than sufficient since 2^9 < 2*sqrt(2^23) + 1. It's easy to show this for double as well with 64-bit integers. So although they can't store all the digits as long as the sqrt of what they can store is accurate then the sqrt(fraction) should be sufficient. Now let's look at what happens for integers close to INT_MAX and the sqrt:

unsigned xi = -1-1;
printf("%u %u\n", xi, (unsigned)(float)xi);  //4294967294 4294967295
printf("%u %u\n", (unsigned)sqrt(xi), (unsigned)sqrtf(xi));  //65535 65536

Since float can't store all the digits of 2^31-2 and double can they get different results for the sqrt. But the float version of the sqrt is one integer larger. This is what I want. For 64-bit integers as long as the sqrt of the double always rounds up it's okay.

解决方案

First, integer multiplication is really quite cheap. So long as you have more than a few cycles of work per loop iteration and one spare execute slot, it should be entirely hidden by reorder on most non-tiny processors.

If you did have a processor with dramatically slow integer multiply, a truly clever compiler might transform your loop to:

for (uint64_t i = 0, j = 0; j < cut; j += 2*i+1, i++)

replacing the multiply with an lea or a shift and two adds.

---

Those notes aside, let’s look at your question as stated. No, you can’t just use i < sqrt(n). Counter-example: n = 0x20000000000000. Assuming adherence to IEEE-754, you will have cut = 0x5a82799, and cut*cut is 0x1ffffff8eff971.

However, a basic floating-point error analysis shows that the error in computing sqrt(n) (before conversion to integer) is bounded by 3/4 of an ULP. So you can safely use:

uint32_t cut = sqrt(n) + 1;

and you’ll perform at most one extra loop iteration, which is probably acceptable. If you want to be totally precise, instead use:

uint32_t cut = sqrt(n);
cut += (uint64_t)cut*cut < n;

---

Edit: z boson clarifies that for his purposes, this only matters when n is an exact square (otherwise, getting a value of cut that is "too small by one" is acceptable). In that case, there is no need for the adjustment and on can safely just use:

uint32_t cut = sqrt(n);

Why is this true? It’s pretty simple to see, actually. Converting n to double introduces a perturbation:

double_n = n*(1 + e)

which satisfies |e| < 2^-53. The mathematical square root of this value can be expanded as follows:

square_root(double_n) = square_root(n)*square_root(1+e)

Now, since n is assumed to be a perfect square with at most 64 bits, square_root(n) is an exact integer with at most 32 bits, and is the mathematically precise value that we hope to compute. To analyze the square_root(1+e) term, use a taylor series about 1:

square_root(1+e) = 1 + e/2 + O(e^2)
                 = 1 + d with |d| <~ 2^-54

Thus, the mathematically exact value square_root(double_n) is less than half an ULP away from[1] the desired exact answer, and necessarily rounds to that value.

[1] I’m being fast and loose here in my abuse of relative error estimates, where the relative size of an ULP actually varies across a binade — I’m trying to give a bit of the flavor of the proof without getting too bogged down in details. This can all be made perfectly rigorous, it just gets to be a bit wordy for Stack Overflow.

这篇关于整数sqrt的准确度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！