卡常(小)技巧

bh1234666

2022-08-03 17:30:13

Tech. & Eng.

由于 cpu 架构的更新及计算机运算速度的加快，本文部分内容不适用于所有评测机或具有时效性。

[WC2017]挑战 是一道十分经典的卡常题，从中我们可以学习到很多卡常知识，且基本不涉及高深的硬件级别的卡常，因此本文从这题开始。

为避免浪费评测资源，本文部分数据来自本地运行。

Subtask 1：

其用到主要卡常思路是更换算法，由于篇幅及作者能力限制，不会涉及具体算法实现。

  1. 实际上，更换算法很多时候并不能称为“卡常”，例如使用 O(n) 的基数排序替换 O(n\log n) 的快速排序，使用 O(n\sqrt {\log n}) 的基数堆替换普通堆等。其本质是使用了时间复杂度更小的算法。

  2. 还有一种更换算法不改变时间复杂度，而是换用常数更小的算法。例如使用树状数组更换线段树，使用红黑树代替 AVL 树等。

  3. 我们甚至可以使用时间复杂度更大的小常数算法替换时间复杂度小的大常数算法。例如使用 O(n\sqrt n) 的分块替换 O(n\log n) 的平衡树。

上述三种方式下面各举出一个例子。

  1. 即本题第一问，在 2e8 的数据范围下基数排序用时四秒多，而快速排序用时十秒多。

  2. 1e6 数据范围的单点加法区间求和，树状数组一百多毫秒，而线段树约五百毫秒。

  3. 【模板】普通平衡树,1e5 的平衡树操作，分块约十毫秒，pbds的平衡树约一百五十毫秒。

可以看出，更换算法带来的速度提升通常是很大的，但是更换算法要求选手有较深厚的知识储备。

实际上这个点还需要用到高速缓存，后面会提到。

Subtask 2：

其用到主要卡常思路是二进制压位。电脑一次能处理的位不止一个，因此很多对于单个位操作的操作都可以使用二进制压位。

二进制压位可以节省空间及提升运算速度。由于电脑最小存储单位是 byte 而非 bit，因此将多个位存储在一个字节可以节省 8 倍的空间。由于使用内存的减小，同时也增加了内存调度的速度。对于位运算相当于可以一次做 8 位，也提升了速度。

常见的例子如 bitset。

由于该 sub 不使用压位耗时过长，这里举出【模板】线性筛素数的例子，对于同样的埃式筛筛取 1e8 的素数，使用 bitset 耗时约七百毫秒, 使用 bool 数组耗时约两秒。

Subtask 3：

该点主要的思路是在复杂度不变的情况下优化算法。使用指针实现整体移动或高速的空间操作替换较慢的通过赋值位移的操作。

对于 n=2666666 的括号序列求解，使用指针实现整体移动用时约两秒，使用赋值操作移动用时三秒多。

由于此题特殊原因（构造数据可以使得优化无效），上述例子优化并不明显。

下面给出更明显的例子，对于【模板】快速排序，使用插入排序实现。

很多人认为使用 vector 实现插入排序十分迅速，其本质其实是 vector 的位移使用了 memmove，实际上 string.h 下的 mem 操作大多十分迅速。

此题使用 memmove 实现 O(n^2) 的插排耗时约三百毫秒，使用赋值移动耗时约七秒（不开 O2，因为 O2 会识别赋值移动并优化成 memmove）。

至此，我们顺利的完成了[WC2017]挑战，除 T1 外并不需要使用到硬件级别优化。

更多常用的卡常技巧

从 【模板】快速阶乘算法 看更高级的卡常技巧

该部分可能涉及部分硬件知识，同时部分技巧优化效果与 cpu 架构及工艺有关。

该题正解为 O(\sqrt n\log n) 的多项式解法，但是我们可以对 O(n) 做法进行卡常。为便于比较速度，我们取模数为 998244353（并非在程序内直接定义，因为固定模数时编译器会优化取模），并计算出 998244352 的阶乘。

直接枚举 1\sim n 相乘。

scanf("%d%d",&n,&mod);
start=clock();
for(int i=1;i<=n;i++)
    ans=1ll*ans*i%mod;
finish=clock();
printf("%d %.3f",ans,(double)(finish-start)/CLOCKS_PER_SEC);

耗时约十二秒。

实际上部分架构的 cpu 支持零开销硬件循环，但是仅仅会节省循环变量调用时间，由于上述代码主要耗时不在循环变量的加减，因此对耗时几乎没影响。

上面提到固定模数时编译器会优化取模，但是由于这里模数是读入的，编译器无法优化，因此我们手动进行优化。

我们采用 barrett 约减（固定模数时编译器也会采用相同的方法优化），其原理是通过成本较低的乘法运算和位运算替换取模（整除）运算。

预处理 brt=\lfloor 2^{r}/len\rfloor，可以通过 brt\times x /2^r 得到 \lfloor x/len \rfloor 的近似值，当 r 足够大时，得到的近似值与精确值至多差 1。

我们可以通过 x-brt\times x /2^r 来得到 x\bmod len 的近似值，若得到的值大于 len 则减 len 即可。

即使用低成本的加减乘及位运算代替高成本的取模运算。

scanf("%d%d",&n,&mod);
start=clock();
__uint128_t brt=((__uint128_t)1<<64)/mod;
for(int i=1;i<=n;i++)
{
    ans*=i;
    ans=ans-mod*(brt*ans>>64);
    while(ans>=mod) ans-=mod;
}
finish=clock();
printf("%d %.3f",ans,(double)(finish-start)/CLOCKS_PER_SEC);

耗时约三秒。

将上述代码 while 改成 if 不会影响结果，但是会使得耗时增加到三点四秒左右。具体原理不明，但作者怀疑是分支预测的影响。

程序中的条件分支是根据程序指令在流水线处理后结果再执行的，所以当CPU等待指令结果时，流水线前级处于空闲状况，因此部分 cpu 会加入分支预测提前执行部分指令来避免流水线的等待以提高 cpu 的运行效率。

但是 barrett 约减中由于精度较高， ans\ge mod 几乎不会成立，此时分支预测提前执行的指令几乎全部作废，因此不会提高效率，反而因为需要将已经装入流水线执行的指令和结果全部清除，然后再装入正确指令重新处理导致效率更低。

有多个分支时可以使用 __builtin_expect(x,0/1) 来告诉编译器 x 更倾向于 0（假）还是 1（真），以控制分支预测走向，不过这里条件判断语句只有一个分支因此这个方法没用。

上述代码使用 while 时不会触发分支预测，因此效率反而高于 if。

循环展开是一种牺牲程序的尺寸来加快程序的执行速度的优化方法，常用来降低循环开销，为具有多个功能单元的处理器提供指令级并行。也有利于指令流水线的调度。

展开可以积极调度（或管道化）循环以掩盖一些延迟。如果有足够的空闲寄存器使变量保持活动状态，因为通过展开相关性链展露了关键路径，这将非常有用。

要特别注意，循环展开时要尽量避免上下句使用相同变量，否则会导致 cpu 并行不充分。

scanf("%d%d",&n,&mod);
start=clock();
const __uint128_t brt=((__uint128_t)1<<64)/mod;
for(i=1;i+8<=n;i+=8)
{
    ans1*=i;ans2*=(i+1);ans3*=(i+2);ans4*=(i+3);ans5*=(i+4);ans6*=(i+5);ans7*=(i+6);ans8*=(i+7);
    ans1=ans1-mod*(brt*ans1>>64);
    ans2=ans2-mod*(brt*ans2>>64);
    ans3=ans3-mod*(brt*ans3>>64);
    ans4=ans4-mod*(brt*ans4>>64);
    ans5=ans5-mod*(brt*ans5>>64);
    ans6=ans6-mod*(brt*ans6>>64);
    ans7=ans7-mod*(brt*ans7>>64);
    ans8=ans8-mod*(brt*ans8>>64);
    if(ans1>=mod) ans1-=mod;
    if(ans2>=mod) ans2-=mod;
    if(ans3>=mod) ans3-=mod;
    if(ans4>=mod) ans4-=mod;
    if(ans5>=mod) ans5-=mod;
    if(ans6>=mod) ans6-=mod;
    if(ans7>=mod) ans7-=mod;
    if(ans8>=mod) ans8-=mod;
}
ans1=ans1*ans2%mod;
ans3=ans3*ans4%mod;
ans5=ans5*ans6%mod;
ans7=ans7*ans8%mod;
ans1=ans1*ans3%mod;
ans5=ans5*ans7%mod;
ans=ans1*ans5%mod;
for(;i<=n;i++)
    ans=ans*i%mod;
finish=clock();
printf("%d %.3f",ans,(double)(finish-start)/CLOCKS_PER_SEC);

耗时约 0.95s。

注意此时我把上面的 while 更换成了 if，因为 while 会影响 cpu 并发。使用 while 耗时约一秒。

循环展开本质是触发 cpu 并行，那是否存在一种办法直接让 cpu 并行呢，答案是肯定的。

AVX2 是一种实现 SIMD 操作的指令集。SIMD的全称为：Single Instruction stream Multiple Data streams，对应的中文名为 单指令流多数据流。SIMD 为并行计算中的一种。

AVX 全称为：Advanced Vector Extensions（又名，Sandy Bridge New Extensions）,是Intel和AMD微服务器x86指令集的extension扩展，AVX2 为其升级版本，支持的位数扩展到了 256 位。

其核心代码较长，这里就不放了，可以参考我在这题提交的题解 link。

耗时约 0.85s。

我们通过一步步加强优化的力度，将耗时从 12s 一步步压到了 0.85s，并且使用 O(n) 的暴力通过了这道正解 O(\sqrt n\log n) 的题。相信从中读者也领会到了卡常的魅力及它的强大。

然后出题人就把时限从原来的 4s 改成了现在的 2s。

当然指令集也是可以循环展开的，但是一般不用，否则码量过大及可读性过差，同时使用指令集和循环展耗时可以到 0.2s 左右于是又通过了此题。

简单考虑阶乘的一些简单的性质，假如我们要求 n 的阶乘，且我们已经通过指令集及循环展开得到了 n/2 的阶乘。

我们发现，由于实际上我们得到的是 8 个形如若干 8n+k 的数分别相乘的结果，那么我们可以很轻易的从中分离出奇偶项。

考虑将 n/2 的阶乘乘 2^{n/2} ，我们就能得到 n 以内所有偶数的乘积，且我们可以从 n/2 的阶乘中分离出所有奇数的乘积，那么我们只需要再求出 n/2 到 n 之间所有奇数的乘积即可。

于是我们又节省了 1/4 的常数，且求出 n/2 的阶乘的过程也可以通过上述做法优化，最终可以优化掉一半的常数。

实际上为了方便，我将 n 先缩小到了 2 的整数次幂，最后再将未乘上去的部分暴力乘上去。可以发现进行第 t 轮优化时缩小的常数为 1/2^{t+1} ，因此轮数过多以后对常数几乎不再产生影响，反而会因为 n 缩的过小导致之后暴力的成本增加，因此我只进行了 8 轮优化。

部分核心代码：

#define mv 8//进行8轮 
int n=N-(N&(256*(1<<mv)-1));//先将n缩小到2的整数次幂 
...//初始化蒙哥马利约减 
unsigned as=mon_in(1),as2=mon_in(1),*fl;//as记录偶数的乘积，as2记录奇数的乘积
{
__m256i ans[8]={...};//将答案初始化为1 
ad=_mm256_set1_epi32(mon_in(64));//将每次步长初始化为64(指令集8位乘循环展开8位) 
__m256i ml[8]={...};//将乘数初始化为1-64 
for(unsigned i=1;i+63<=n>>mv;i+=64)
{
    ans[0]=mul(ans[0],ml[0]);
    ...
    ans[7]=mul(ans[7],ml[7]); 
    ml[0]=add(ml[0],ad);
    ...
    ml[7]=add(ml[7],ad);
}//得到前n>>mv位的乘积
fl=(unsigned*)(ans+0);
as=mul(as,mul(mul(fl[1],fl[3]),mul(fl[5],fl[7])));//分离偶数 
as2=mul(as2,mul(mul(fl[0],fl[2]),mul(fl[4],fl[6])));//分离奇数 
...
fl=(unsigned*)(ans+7);
as=mul(as,mul(mul(fl[1],fl[3]),mul(fl[5],fl[7])));
as2=mul(as2,mul(mul(fl[0],fl[2]),mul(fl[4],fl[6])));
}
for(int j=mv-1;j>=0;j--)
{
as=mul(as,as2);as=mul(as,mon_in(qpow(2,n>>(j+1),p)));//as偶数项,as2前半奇数项 
__m256i ans[8]={...};//初始化ans为1 
ad=_mm256_set1_epi32(mon_in(128));//初始化步长为128(从这里开始仅需处理奇数位了) 
const unsigned add_=n>>(j+1);//初始化起始位置为n>>(j+1)，即上一次结束的位置 
__m256i ml[8]={...};//初始化乘数为从add_开始的64个奇数 
for(unsigned i=add_;i+127<=n>>j;i+=128)
{
    ans[0]=mul(ans[0],ml[0]);
    ...
    ans[7]=mul(ans[7],ml[7]);
    ml[0]=add(ml[0],ad);
    ...
    ml[7]=add(ml[7],ad);
}//得到当前处理段后半奇数的乘积 
fl=(unsigned*)(ans+0);
as2=mul(as2,mul(mul(mul(fl[0],fl[1]),mul(fl[2],fl[3])),mul(mul(fl[4],fl[5]),mul(fl[6],fl[7]))));
...
fl=(unsigned*)(ans+7);//合并当前处理段后半奇数的乘积 
as2=mul(as2,mul(mul(mul(fl[0],fl[1]),mul(fl[2],fl[3])),mul(mul(fl[4],fl[5]),mul(fl[6],fl[7]))));
}
as=mul(as,as2);//将奇数偶数部分相乘得到答案 
as=mon_out(as);//离开数域 
for(int i=n+1;i<=N;i++)//暴力将最后一段乘上去 
    as=1ull*as*i%p;
printf("%u\n",as);//输出答案

于是我们又卡掉了一半的常数，耗时达到了约 0.1s，在原题中总时间 5.11s，已经超过了大部分 O(\sqrt n\log n) 的正解。

从最开始的直接循环枚举取模的 12s，到最后的指令集加循环展开加小技巧的 0.1s，我们通过各种技巧卡掉了 120 倍的常数同时码长翻了 200 倍，让一个 O(n) 的暴力吊打了常数较大的 O(\sqrt n\log n) 的正解。

其他高级卡常技巧

其全称为高速缓冲存储器，高速缓冲存储器是存在于主存与CPU之间的一级存储器，由静态存储芯片(SRAM)组成，容量比较小但速度比主存高得多，接近于CPU的速度。

部分题目我们可以通过增加部分操作运行时间压缩数组，使得数组大小得以进入高速缓存。高速缓存的速度比内存快上百倍，因此频繁调用内存时可以尝试减少调用的内存空间，增加调用次数，使得访问的数据尽量位于高速缓存中。

使用实例例如基数排序。经典的基数排序分松式基排和紧式基排，松式基排进行了四次循环（假设数据范围为 32 位无符号整型），但是计数数组大小为 256，紧式基排进行了两次循环，但是计数数组为 65536。

试验结果表明松式基排虽然循环次数较多，但是速度远快于紧式基排（对于部分cpu）。

上文提到的 [WC2017]挑战 subtask1 就需要使用这个技巧，紧式基排无法通过。

对于部分极为连续的数据，我们甚至可以将本身就为 0 的数据（全局数组）再 memset 一遍，memset 本身耗时极短可以忽略不计，但是会将数据放入高速缓存，提升后续访问速度。

这里主要列举出 x86 架构下常见的指令集，感兴趣的读者可以自行查找资料(微软上有一些资料)。

SIMD 系列:MMX SSE SSE2 SSE3 AVX AVX2，上述指令集指令均类似，多用于整型及浮点型数据的批量运算。

通用指令:数据传输指令，二进制算术运算指令，十进制算数运算指令，逻辑指令，移位和旋转指令，位和字节指令，控制传输指令，串指令，I/O指令，进入和离开指令，标志控制指令，段寄存器指令等。

流扩展指令：水平加减指令，批量求绝对值，批量乘、加等

cpu(x86 系列 64 位处理器)更多相关信息可以参考《64位微处理器系统编程和应用编程》（感觉这本书写的挺好的至少它是中文的）。

对于内存调度方面实际上有极大的卡常空间，但是本人对于内存调度了解甚少，这里不做赘述，感兴趣的读者可以自行查询资料。

OI 内禁止使用，不做赘述。

写在最后

卡常是一门极其深奥且有魅力的学问，深入了解卡常可以让我们更加了解计算机，在写代码时自觉或不自觉的采用更为快速的写法。

同时卡常也可以让我们完成一些远超自己实力的题目，例如上文提到的快速阶乘，我们避开了困难的多项式，使用暴力通过了此题。比赛时遇到不会的题也可以尝试对暴力卡常得到意想不到的分数。

最后，卡常是永无止境的。举一个例子，之前作者给一份代码卡常，发现看似完全不会影响代码运行时间的操作会对代码产生极大的影响。例如对 namespace 内的变量加上 static 关键字，交换 if else 的顺序等。

一份代码看似已经卡常彻底，但是可能在某个细枝末节的地方进行略微修改，就会导致代码效率有质的飞跃。

本人只是个高二的 OIER，对计算机并不是非常了解，上面只是讲了我知道的一些卡常技巧，卡常无穷无尽，还有更多技巧等着大家探索。