整数sqrt的准确性

Question

我有一个像这样的循环:

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是一个完美的正方形,那么n = y^2我的问题就是 - 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  //solve
dx < 2*sqrt(x) + 1 
// e.g for x = 100 dx < 21
// sqrt(100+20) < sqrt(100) + 1

Float可以存储23位,所以我们让y = 2 ^ 23 + 2 ^ 9.这已经足够了,因为2 ^ 9 <2*sqrt(2 ^ 23)+ 1.很容易将此显示为双精度以及64位整数.因此,尽管只要他们可以存储的sqrt是准确的,他们就不能存储所有数字,那么sqrt(分数)就足够了.现在让我们看看接近INT_MAX和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

由于float不能存储2 ^ 31-2的所有数字而且double可以为sqrt获得不同的结果.但是sqrt的float版本是一个更大的整数.这就是我要的.对于64位整数,只要double的sqrt总是向上舍入就可以了.

Answer 1

首先,整数乘法非常便宜.只要每个循环迭代和一个备用执行槽有多个工作周期,就应该通过在大多数非微小处理器上重新排序来完全隐藏它.

如果你的处理器具有极慢的整数乘法,那么一个真正聪明的编译器可能会将你的循环变换为:

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

用一个lea或一个移位替换乘法,然后加两个.

除了这些注释之外,让我们看看你提出的问题.不,你不能只是使用i < sqrt(n).反例:n = 0x20000000000000.假设遵守IEEE-754,您将拥有cut = 0x5a82799,而且cut*cut是0x1ffffff8eff971.

但是,基本的浮点错误分析表明计算中的错误sqrt(n)(在转换为整数之前)受到ULP的3/4的限制.所以你可以安全地使用:

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:

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

因此,数学上精确的值square_root(double_n)小于ULP的一半,远离[1]期望的精确答案,并且必然舍入到该值.

[1]我在这里滥用相对误差估计是快速和宽松的,其中ULP的相对大小实际上在一个binade中变化 - 我试图给出一些证据的味道而不会太过于陷入困境细节.这一切都可以完全严格,它只是对Stack Overflow有点罗嗦.