使用 sqrt 和 floor 时的近似错误

Question

我必须枚举一个方程的解，我知道y < x *( sqrt(n) - 1 )，其中 x、y 和 n 是整数。

我天真的方法是寻找小于或等于floor( x * ( sqrt( (float)n ) - 1 ) ) 的y。

---

• 我应该担心近似误差吗？

• 例如，如果我的表达式有点大于整数m，我应该担心最后得到m-1吗？

• 如何检测此类错误？

Answer 1

您绝对应该担心近似误差，但担心程度取决于您关注的 x 和 n 值的范围。

在 IEEE 4 字节浮点表示中的计算会出现大约 2^23 到 2^24 的一部分的错误；对于 8 字节表示（即double），它将大约是 2^52 到 2^53 的一部分。您可能会期望您需要使用 doubles 而不是 floats 来获得 32 位整数 x 和 n 的准确结果，并且对于 64 位整数，即使是 double 也不够。

例如，考虑代码：

template <typename F,typename V>
F approxub(V x,V n) {
    return std::floor(x*std::sqrt(F(n))-x);
}

uint64_t n=1000000002000000000ull; // (10^9 + 1)^2 - 1
uint64_t x=3;
uint64_t y=approxub<double>(x,n);

这给出了 y=3000000000 的值，但正确的值是 2999999999。

当 x 很大而 n 很小时，情况更糟：在 IEEE doubles 中不能精确表示大的 64 位整数：

uint64_t n=9;
uint64_t x=5000000000000001111; // 5e18 + 1111
uint64_t y=approxlb<double>(x,n);

y 的正确值（将 n 是一个完美正方形的问题放在一边——在这种情况下，真正的上限将小一）是 2 x em> = 10000000000000002222，即 1e19 + 2222。但是，计算出的 y 是 10000000000000004096。

避免浮点近似

假设您有一个函数isqrt，它精确计算整数平方根的整数部分。那你可以说

y = isqrt(x*x*n) - x

如果产品 x*x*n 适合您的整数类型，您将有一个精确的上限（如果 n 是一个完美的正方形，则比上限多一个。）编写isqrt 函数的一种方法；这是一个基于material at code codex的示例实现：

template <typename V>
V isqrt(V v) {
    if (v<0) return 0;

    typedef typename std::make_unsigned<V>::type U;
    U u=v,r=0;

    constexpr int ubits=std::numeric_limits<U>::digits;
    U place=U(1)<<(2*((ubits-1)/2));

    while (place>u) place/=4;
    while (place) {
        if (u>=r+place) {
            u-=r+place;
            r+=2*place;
        }
        r/=2;
        place/=4;
    }
    return (V)r;
}

如果 x 太大了怎么办？例如，如果我们最大的整数类型有 64 位，并且 x 大于 2^32。最直接的解决方案是进行二分搜索，以 x r - x 和 x r 为边界，其中 r = [√ n] 是整数平方根。