整数sqrt的准确性
我有一个像这样的循环：
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 // eg 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的所有数字，因此它们可以为sqrt获得不同的结果。 但是sqrt的float版本是一个更大的整数。 这就是我要的。 对于64位整数，只要double的sqrt总是向上舍入就可以了。
首先，整数乘法非常便宜。 只要每个循环迭代和一个备用执行槽有多个工作周期，就应该通过在大多数非微小处理器上重新排序来完全隐藏它。
如果你的处理器具有极慢的整数乘法，那么真正聪明的编译器可能会将你的循环转换为：
for (uint64_t i = 0, j = 0; j < cut; j += 2*i+1, i++)
用lea或shift替换乘法，然后加两次。
除了这些注释，让我们看看你提出的问题。 不，你不能只使用i < sqrt(n) 。 反例： n = 0x20000000000000 。 假设遵守IEEE-754，你将cut = 0x5a82799 ， cut*cut为0x1ffffff8eff971 。
但是，基本的浮点错误分析表明，计算sqrt(n) （转换为整数之前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 |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]所需的精确答案，并且必然square_root(double_n)到该值。
[1]我在这里滥用相对误差估计是快速和宽松的，其中ULP的相对大小实际上在一个binade中变化 - 我试图给出一些证据的味道而不会太过于陷入困境详细说明。 这一切都可以完全严格，它只是有点罗嗦Stack Overflow。
如果您能够访问IEEE 754双精度浮点数，那么我的所有答案都是无用的，因为Stephen Canoncertificate了这两点
- 一种避免imul in循环的简单方法
- 计算天花板sqrt的简单方法
否则，如果由于某种原因你有一个非IEEE 754兼容平台，或者只有一个单精度，你可以用一个简单的Newton-Raphson循环得到平方根的整数部分。 例如在Squeak Smalltalk中，我们在Integer中有这个方法：
sqrtFloor "Return the integer part of the square root of self" | guess delta | guess := 1 bitShift: (self highBit + 1) // 2. [ delta := (guess squared - self) // (guess + guess). delta = 0 ] whileFalse: [ guess := guess - delta ]. ^guess - 1
其中//是整数除法商的运算符。
最终后卫guess*guess <= self ifTrue: [^guess]. 如果初始猜测超过精确解，则可以避免这种情况。
使用近似浮点sqrt初始化不是一个选项，因为整数是任意大的并且可能溢出
但是在这里，你可以用浮点sqrt近似来初始猜测，我的赌注是在非常少的循环中找到确切的解。 在C中将是：
uint32_t sqrtFloor(uint64_t n) { int64_t diff; int64_t delta; uint64_t guess=sqrt(n); /* implicit conversions here... */ while( (delta = (diff=guess*guess-n) / (guess+guess)) != 0 ) guess -= delta; return guess-(diff>0); }
这是一些整数乘法和除法，但在主循环之外。
您正在寻找的是一种计算自然数的平方根的有理上界的方法。 你需要继续分数[维基百科] [1]。
对于x> 0，有[！[平方根公式] [2]] [2]。
为了使符号更紧凑，重写上述公式为
通过在每个递归深度处移除尾项（x-1）/ 2来截断连续分数，得到一个近似sqrt（x）的序列，如下所示：
上限出现在具有奇数行数的行上，并且变得更紧。 当上限与其相邻下界之间的距离小于1时，该近似值就是您所需要的。 使用该值作为cut的值，此处cut必须是一个浮点数，才能解决问题。
对于非常大的数，应该使用有理数，因此在整数和浮点数之间的转换期间不会丢失精度。