我的乘法取模方式果然有问题

· · Tech. & Eng.

众所周知，模数为 long long 的乘法取模可以使用 __int128 ，也可以使用一些手法操作一下（即浮点数运算得到下取整的结果并使用自然溢出），后者的具体实现如下面这份代码所示。

inline long long mul(long long a,long long b)
{
    long long ret=a*b-((long long)(double(a)/Mod*b+0.5))*Mod;
    return ret+(ret<0)*Mod;
}

特别的，我们默认模数不是特别大（如 2^{50} 左右），否则需要使用 long double 提高精度。

考虑下面这份代码，它比较了上述两种乘法取模的实现方式的运行效率。

#include<iostream>
#include<random>
using namespace std;
mt19937_64 Rand(123);
const int N=3e4,V=50;
const long long mod=(1ll<<V)+123456789;
long long Mod=mod,a[N+10],b[N+10];
inline long long mul(long long a,long long b)
{
    long long ret=a*b-((long long)(double(a)/Mod*b+0.5))*Mod;
    return ret+(ret<0)*Mod;
}
inline long long Mul(long long a,long long b)
{
    return __int128(a)*b%mod;
}
int main()
{
    for(int i=1;i<=N;i++)
    {
        a[i]=Rand()&(1ll<<V)-1;
        b[i]=Rand()&(1ll<<V)-1;
    }
    long long ans=0;
    clock_t Start=clock();
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++)
            ans^=mul(a[i],b[j]);
    cout<<clock()-Start<<endl;
    Start=clock();
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++)
            ans^=Mul(a[i],b[j]);
    cout<<clock()-Start<<endl;
    cout<<ans<<endl;
    return 0;
}

在 AtCoder 的 Custom Test 上测试，语言选的是 C++ 23 (Clang 16.0.6) 。

运行后输出如下结果。

362830
3579920
0

（注：该评测机的 CLOCKS_PER_SEC 是 10^6 ，即总运行时间不到 4 秒。）

不难发现， __int128 的运行时间大约是 double 的十倍甚至九倍，我对此还是比较惊讶的。

顺带一提，就算是使用 long double 也是比 __int128 快的。

使用上述评测环境，输出结果如下，大概是 __int128 运行时间的三分之二。

2534783
3675875
0

我们回到第一份代码的输出结果，仔细观察后发现， double 竟然只跑了 0.36 秒左右。

这是什么概念？这又意味着什么呢？

我们在没有使用循环展开等奇技淫巧且只开启了 O2 优化的情况下，计算 (3\times 10^4)^2=0.9\times 10^9 次模数为 long long 的乘法取模只使用了 0.36 秒！换句话说，在一秒的时间内，我们的程序可以进行约 2.5\times 10^9 次模数为 long long 乘法取模！

你可能对乘法取模的速度没什么概念，那么我们来看看我们熟知的模数在 int 级别的乘法取模。

#include<iostream>
#include<random>
using namespace std;
mt19937 Rand(123);
#define rand Rand
const int N=3e4,V=29,mod=1e9+7;
int Mod=mod,a[N+10],b[N+10];
inline int mul(int a,int b)
{
    int ret=a*b-(int(double(a)/Mod*b+0.5))*Mod;
    return ret+(ret<0)*Mod;
}
inline int Mul(int a,int b) 
{
    return (long long)a*b%mod;
}
int main()
{
    for(int i=1;i<=N;i++)
    {
        a[i]=Rand()&(1<<V)-1;
        b[i]=Rand()&(1<<V)-1;
    }
    int ans=0;
    clock_t Start=clock();
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++)
            ans^=mul(a[i],b[j]);
    cout<<clock()-Start<<endl;
    Start=clock();
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++)
            ans^=Mul(a[i],b[j]);
    cout<<clock()-Start<<endl;
    cout<<ans<<endl;
    return 0;
}

使用同样的评测环境，运行后的输出结果如下。

257907
963987
0

？？？

我不懂，但我大受震撼。

下面是笔者自己的一点猜测。

我怀疑编译器（ Clang ）进行了一些操作使得大量对 mul() 函数的调用得以并行计算。

可能与循环展开本质类似。

（注：疑似该程序使用 gcc 编译器时的运行效率会大大降低。）

（更新于 2025/06/21 。）

昨天凌晨破防的时候我发现，由于 a[i] 的取值与 j 无关，故编译器可能会优化掉一些重复的运算（如 double(a)/Mod 等），对此我们可以考虑改为调用 mul(a[i]^b[j],b[j]) 。若编译环境保持不变，上述三份代码运行后的输出结果依次如下。

554456
3351854
0
3485189
3347934
0
526950
1211967
0

不难发现，更改后运行效率确实有显著降低，看上去也更符合我对乘法取模的想象了。

（虽然感觉还是有点快。）

特别的， long double 实现的乘法取模不再能跑过 __int128 。

（注：现在疑似该程序使用 gcc 编译器时的运行效率与使用 Clang 时无显著区别。）

详细揭秘：关于我在无意间被浮点数的指令集变成废柴这件事。

wow 好厉害。