关于O(1)快速乘和关于其特判的原因

· · 个人记录

O(1)快速乘

在模超过 int 范围的大质数的意义下, 乘法运算可能会超过 long long 的范围. 在不能使用 __int128 的情况下, 我们可以使用 \mathcal O(\log n) 的慢速乘来避免爆 long long .

long long mul(long long a, long long b) {
    long long rt = 0;
    while (b) {
        if (b) (rt += a) %= mod;
        (a += a) %= mod, b >>= 1;
    }
    return rt;
}

然而这个让复杂度带上一个 \log 的方法并不令人十分满意. 人类运用智慧, 利用模数在 long long 范围内的性质, 发明了一个复杂度 \mathcal O(1) 的方法(下称 \mathcal O(1) 快速乘).

long long mul(long long a, long long b) {
    long long rt = a * b - ((long long)((long double)a / mod * b + 0.5)) * mod;
    return rt < 0 ? rt + mod : rt;
}

上述代码的原理在于 a\bmod b=a-\lfloor \frac{a}{b}\rfloor \times b .

由于代码中 (long long)((long double)a / mod * b + 0.5)) 肯定在 long long 范围内, 并且乘上 mod 后与 a * b 的差不超过 mod, 故正确性是有保证的.

特判

注意到代码中有两个奇怪的地方:

\

可以想到返回值的特判的原因是 "五入" 导致减数比向下取整要大了 mod , 那么我们将代码改为如下 Code 1 的向下取整而不特判的形式是否是可以的呢?

// Code 1
long long mul(long long a, long long b) {
    long long rt = a * b - ((long long)((long double)a / mod * b)) * mod;
    return rt;
}

本人经过测试, 得到的答案是否定的, 在网上查找资料也未能得出结论, 网上的某些错误结论甚至使我在错误的道路越走越远.

所幸最后我还是找到了问题的结症所在, 我将向下取整的代码调整一下后的得到了如下正确的代码:

long long mul(long long a, long long b) {
    long long rt = a * b - ((long long)((long double)a / mod * b)) * mod;
    if (rt >= mod) rt -= mod;
    if (rt < 0) rt += mod;
    return rt;
}

这说明我们直接得出的答案与正确值之差的绝对值是不大于 mod 的, 并且是它的一个整数倍, 这令我感到鼓舞, 冷静分析之后, 发现了 Code 1 错误的原因——精度误差.

考虑当 (long double)a / mod * b 的真实值是一个十分接近但小于某个整数 x 的值, 在计算过程中的精度损失会使其超过 x 并且使向下取整后的结果大 1 , 同理其十分接近但大于某个整数会使结果小 1 , 故我们得到的答案与正确值之差的绝对值是 mod 或 0.

至此, 特判的原因已经完全清楚了.

回过头来看最初的的写法, 我们意识到 + 0.5 的真正原因不是 '四舍五入', 而是使计算出来的值不会比原来大, 这样就只需要特判一种情况, 从而达到减小常数的目的.

总而言之, 挺妙的*<|:-).