64位整数乘法

前置知识：

• 在一些平台上，long long 可以无损转 long double（但是无论如何都不可以无损转 double）
• 对于无符号整数来说，如果\(x<y\)，那么\(x-y\)会得到一个极大的数，相当于绕了一圈

然后可以得到正确代码如下：

#include<bits/stdc++.h>
#define ll long long
#define ull unsigned long long
using namespace std;
ll mul(ull a, ull b, ull p)
{
    a %= p;
    b %= p;

    ull c = (long double)a * b / p;
    ull res = a * b - c * p;

    return (ll)((res + p) % p);
}
int main()
{
    ll a,b,p;
    scanf("%lld%lld%lld",&a,&b,&p);
    printf("%lld",mul(a,b,p));
    return 0;
}

这是因为：
整个计算过程可以写成：

\[\widehat m=\operatorname{RN}(a\cdot b) \]

\[\widehat x=\operatorname{RN}\left(\frac{\widehat m}{p}\right) \]

\[c=\operatorname{trunc}(\widehat x) \]

其中 RN 表示 IEEE 754 默认的“舍入到最近值，距离相同时取偶数”。

算法并不保证 \(c\) 一定是正确商，只保证它与正确商最多相差 1，然后再通过加减一次 \(p\) 修正。

1. a、b 和 p 的转换是否有误差

在常见的 x86 80 位扩展精度 long double 中，有 64 个二进制有效位。

而一个有符号 64 位整数的绝对值最多约为：

\[2^{63} \]

64 个二进制有效位足以精确表示所有 long long 整数。因此在这种平台上：

(long double)a

是精确转换，没有误差。

表达式：

(long double)a * b / p

执行时，b 和 p 也会先隐式转换为 long double，因此实际运算是：

(long double)a * (long double)b / (long double)p

所以转换阶段没有丢失整数位。

2. a*b 最多有 126 位，为什么只保留 64 位仍然可以

假设：

\[0\le a,b<p<2^{63} \]

那么精确乘积 \(ab\) 最多需要约 126 个二进制位。

例如精确乘积规格化后可能是：

\[ab=1.b_1b_2\cdots b_{63}b_{64}b_{65}\cdots\times 2^E \]

80 位扩展精度只能保留 64 个有效位：

\[1.b_1b_2\cdots b_{63}\times 2^E \]

后面的位确实会被丢弃，但不是简单截断。处理器会查看 guard、round、sticky 等尾部信息，把结果舍入到最近的可表示值。

设舍入后的乘积为：

\[\widehat m=ab$1+\delta_1$ \]

对于 64 位有效精度，在“舍入到最近值”模式下：

\[|\delta_1|\le 2^{-64} \]

也就是说，乘法虽然丢掉很多低位，但它的相对误差非常小。

3. 除以 \(p\) 后，乘法误差会变成多大

真正关心的是商：

\[x=\frac{ab}{p} \]

因为 \(a,b<p\)，所以：

\[ab<p^2 \]

从而：

\[x=\frac{ab}{p}<p<2^{63} \]

舍入后的乘积除以 \(p\)：

\[\frac{\widehat m}{p}=\frac{ab(1+\delta_1)}{p} x(1+\delta_1) \]

乘法舍入造成的商的绝对误差是：

\[\left|\frac{\widehat m}{p}-x\right|= x|\delta_1| \]

因为：

\[x<2^{63},\qquad |\delta_1|\le 2^{-64} \]

所以：

\[x|\delta_1|<2^{63}\cdot 2^{-64}=\frac12 \]

这就是最重要的一步：

a*b 虽然可能丢失很多低位，但这些误差除以 \(p\) 后，对最终商造成的误差小于 \(0.5\)。

这里并不需要恢复精确的 \(ab\)，只需要商的误差足够小。

4. 除法本身还会再舍入一次

处理器还需要计算：

\[\widehat x= \operatorname{RN}\left(\frac{\widehat m}{p}\right) \]

除法也会产生一次相对误差。设：

\[\widehat x=\frac{\widehat m}{p}(1+\delta_2), \qquad |\delta_2|\le2^{-64} \]

代入 \(\widehat m=ab(1+\delta_1)\)：

\[\widehat x=\frac{ab}{p}(1+\delta_1)(1+\delta_2) \]

也就是：

\[\widehat x=x(1+\delta_1)(1+\delta_2) \]

因此：

\[|\widehat x-x| \le x\left(2\cdot 2^{-64}+2^{-128}\right) \]

因为 \(x<p\le2^{63}-1\)，所以：

\[|\widehat x-x| < (2^{63}-1) \left(2^{-63}+2^{-128}\right) <1 \]

因此最终浮点商与精确商满足：

\[\boxed{|\widehat x-x|<1} \]

这才是算法正确性的核心。

5. 最后转换成 long long 才是截断

代码：

long long c = (long double)a * b / p;

右侧先完成浮点乘除，得到 \(\widehat x\)，然后转换为 long long。

因为这里是非负数，所以向零截断等价于下取整：

\[c=\lfloor\widehat x\rfloor \]

令精确商为：

\[x=\frac{ab}{p}=q+\frac rp \]

其中：

\[q=\left\lfloor\frac{ab}{p}\right\rfloor, \qquad 0\le r<p \]

由于：

\[|\widehat x-x|<1 \]

所以截断后的 \(c\) 最多只可能是：

\[c\in{q-1,q,q+1} \]

注意，这里并没有保证：

\[c=q \]

算法只保证最多偏差 1。

因此数学意义的 \(ab-cp\) 的范围在 \([-p,2p)\)，所以最后对于 res 来说需要加上 p 之后再取模