尽快比较形式（a + sqrt（b））的两个值？

Question

作为我编写的程序的一部分，我需要比较形式为a + sqrt(b)where a和bunsigned integer的两个值。由于这是紧密循环的一部分，因此我希望此比较尽可能快地运行。（如果重要的话，我正在x86-64机器上运行代码，并且无符号整数不大于10 ^ 6。此外，我知道这样的事实a1<a2。）

作为独立功能，这就是我要优化的功能。我的数字是足够小的整数，可以double（或什至float）精确地表示它们，但是sqrt结果中的舍入误差不能改变结果。

// known pre-condition: a1 < a2  in case that helps
bool is_smaller(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
    return a1+sqrt(b1) < a2+sqrt(b2);  // computed mathematically exactly
}

测试用例：is_smaller(900000, 1000000, 900001, 998002)应返回true，但如@wim注释所示，sqrtf()将其返回false。因此将(int)sqrt()截断为整数。

a1+sqrt(b1) = 90100和a2+sqrt(b2) = 901000.00050050037512481206。最接近的浮点数就是90100。

---

由于sqrt()即使在现代的x86-64上作为sqrtsd指令完全内联时，该函数通常也非常昂贵，所以我试图避免sqrt()尽可能地调用。

通过平方运算删除sqrt还可以通过使所有计算都精确来避免舍入错误的任何危险。

如果相反，功能是这样的...

bool is_smaller(unsigned a1, unsigned b1, unsigned x) {
    return a1+sqrt(b1) < x;
}

...那么我可以做 return x-a1>=0 && static_cast<uint64_t>(x-a1)*(x-a1)>b1;

但是现在既然有两个sqrt(...)术语，我就不能进行相同的代数运算。

我可以使用以下公式对值进行两次平方：

      a1 + sqrt(b1) = a2 + sqrt(b2)
<==>  a1 - a2 = sqrt(b2) - sqrt(b1)
<==>  (a1 - a2) * (a1 - a2) = b1 + b2 - 2 * sqrt(b1) * sqrt(b2)
<==>  (a1 - a2) * (a1 - a2) = b1 + b2 - 2 * sqrt(b1 * b2)
<==>  (a1 - a2) * (a1 - a2) - (b1 + b2) = - 2 * sqrt(b1 * b2)
<==>  ((b1 + b2) - (a1 - a2) * (a1 - a2)) / 2 = sqrt(b1 * b2)
<==>  ((b1 + b2) - (a1 - a2) * (a1 - a2)) * ((b1 + b2) - (a1 - a2) * (a1 - a2)) / 4 = b1 * b2

无符号除以4很便宜，因为它只是一个移位，但是由于我将数字平方两次，所以我将需要使用128位整数，并且需要引入一些>=0检查（因为我正在比较不平等而不是相等）。

感觉可能有一种更快的方法，通过对这个问题应用更好的代数。有没有办法更快地做到这一点？

Answer 1

Here's a version without sqrt, though I'm not sure whether it is faster than a version which has only one sqrt (it may depend on the distribution of values).

Here's the math (how to remove both sqrts):

ad = a2-a1
bd = b2-b1

a1+sqrt(b1) < a2+sqrt(b2)              // subtract a1
   sqrt(b1) < ad+sqrt(b2)              // square it
        b1  < ad^2+2*ad*sqrt(b2)+b2    // arrange
   ad^2+bd  > -2*ad*sqrt(b2)

Here, the right side is always negative. If the left side is positive, then we have to return true.

If the left side is negative, then we can square the inequality:

ad^4+bd^2+2*bd*ad^2 < 4*ad^2*b2

这里要注意的关键是，如果a2>=a1+1000，则is_smaller始终返回true（因为的最大值sqrt(b1)是1000）。如果a2<=a1+1000，ad则为小数，因此ad^4将始终适合64位（无需128位算术运算）。这是代码：

bool is_smaller(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
    int ad = a2 - a1;
    if (ad>1000) {
        return true;
    }

    int bd = b2 - b1;
    if (ad*ad+bd>0) {
        return true;
    }

    int ad2 = ad*ad;

    return (long long int)ad2*ad2 + (long long int)bd*bd + 2ll*bd*ad2 < 4ll*ad2*b2;
}

编辑：正如彼得·科德斯（Peter Cordes）所注意到的，第一个if不是必需的，因为第二个可以处理它，因此代码变得越来越小和更快：

bool is_smaller(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
    int ad = a2 - a1;
    int bd = b2 - b1;
    if ((long long int)ad*ad+bd>0) {
        return true;
    }

    int ad2 = ad*ad;
    return (long long int)ad2*ad2 + (long long int)bd*bd + 2ll*bd*ad2 < 4ll*ad2*b2;
}