a*b problem 卡常
出了 P8482 Number 这题之后，有人问我 spj 是咋写的。实际上这题的 spj 使用的是暴力卷积，将两个 5e6 级别的数直接相乘。
实际上，该题的数据只有读入数据，也就是答案其实是 spj 现场跑的，所以说 spj 在 1s 内完成了两轮
该博客不会对卡常技巧进行具体描述，对于具体描述可见卡常(小)技巧。
首先，高精度乘法是基于卷积的，因此我们需要有一个快“亿”点的卷积模板，我使用的卷积模板是自己写的 ntt，相较于 P3803 的最优解仍旧慢了 1.3 倍，因此实际上 spj 还有很大的优化空间。
下面将对 spj 使用的 ntt 的优化进行讲解。
注意到 ntt 过程中会大量使用到原根，因此我们预处理原根。
我们知道，ntt 的本质是带入 n 个点将多项式的系数表示法转换为点值表示法，再将点值表示法转换为系数表示法，与带入顺序无关，只要两次转换互相对称即可。因此我们的点值表示无需顺序带入，存在一种非顺序带入的方式会使得运算更为简便，我们可以使用 dif dit 进行表示法的转换。
注意到 ntt 过程中大量使用定模数取模，且所有数据均在一个
对于每次 dif/dit 有大量连续的数进行相同的运算，因此我们可以使用单指令流多数据流指令集进行批量操作。
另外，对于余下的三轮连续数字不满 8 个的循环，AVX2 支持打包及拆分操作，可以将数据重排使得需要运算的数字连续。
于是我们就得到了一个跑的相对较快的 ntt。下面给出 dif 的部分代码，dit 类似。
for(unsigned i=lim>>1,i2=lim;i>=8;i2=i,i>>=1)
for(ty *p=num,*pi;p!=end;p+=i)
{
__m256i num1,num2,num3,num4,num5;
ty *rt=(ty*)root+i;
for(unsigned j=i;j+8<=i2;rt+=8,p+=8,j+=8)//对连续的数进行8个一组的指令集批量运算
{
pi=p+i;
num1=_mm256_loadu_si256((__m256i*)p);
num2=_mm256_loadu_si256((__m256i*)pi);
_mm256_storeu_si256((__m256i*)p,add(num1,num2));
num4=_mm256_sub_epi32(num1,num2);
num5=_mm256_cmpgt_epi32(a0,num4);
num3=_mm256_add_epi32(_mm256_and_si256(num5,mod2),num4);
num5=_mm256_loadu_si256((__m256i*)rt);
num1=num3;num2=num5;
num2=_mm256_mul_epu32(num1,num2);
num1=_mm256_mul_epu32(_mm256_mul_epu32(num2,R),mod1);
num4=_mm256_srli_epi64(_mm256_add_epi64(num1,num2),32);
num1=_mm256_srli_si256(num3,4);num2=_mm256_srli_si256(num5,4);
num2=_mm256_mul_epu32(num1,num2);
num1=_mm256_mul_epu32(_mm256_mul_epu32(num2,R),mod1);
num1=_mm256_and_si256(_mm256_add_epi64(num1,num2),hi32);
_mm256_storeu_si256((__m256i*)pi,_mm256_or_si256(num1,num4));
}
}
for(unsigned i=4,i2=8;i;i2=i,i>>=1)
{
if(i==4)//对零散的进行重排再批量操作
{
ty *p=num;
__m256i _rt=_mm256_set_epi32(0,0,0,0,(ty)root[7],(ty)root[6],(ty)root[5],(ty)root[4]),num1,num2,num3,num4;
for(;p+8<=end;p+=8)
{
num1=_mm256_loadu_si256((__m256i*)p),num2=_mm256_loadu_si256((__m256i*)(p+4));
num3=add(num1,num2);
num4=sub(num1,num2);
num4=mul(num4,_rt);
memcpy(p,&num3,16);
memcpy(p+4,&num4,16);
}
}
//以下省略对 i=2,i=1 时操作的代码有了一份较为快速的卷积模板，我们就可以开始写高精度乘法了。
首先考虑最朴素的高精度乘法，将两个数字直接卷积，然后处理进位。这种方法在这里显然是不可行的，即使不考虑速度，能做 1e7 位的 ntt 模数也难以找到。
因此考虑压位来减少需要卷积的长度，我采用了压 8 位高精，但是压 8 位高精又会带来另一个问题，单模 ntt 的精度不够，所以采用三模 ntt。当然，即使不存在长度问题，三模 ntt 仍快于直接 ntt，三模 ntt 用三倍常数换来了序列长度除八，最终会比直接 ntt 快一倍左右。
不同于拆系数 fft,三模 ntt 几乎不存在（至少我不知道）减少 ntt 次数的办法。因此我们直接进行 3 次卷积。
接下来进行三模合并，发现会进行 n 次的对 _uint128_t 的取模，十分慢，用时甚至超过了一次
下面是我写的 uint256 的板子（由于除法不会写且用的不多直接进行了暴力实现）。
struct __uint256_t{
__uint128_t hi,lo;
#define il __attribute__((always_inline))//强制内联，用编译时间换运行时间
__uint256_t(){}
template<typename ty>
__uint256_t(ty x){hi=0;lo=x;}
__uint256_t(__uint128_t x,__uint128_t y){hi=x;lo=y;}
template<typename ty>
il void operator = (ty x){hi=0;lo=x;}
il void operator = (__uint256_t x){hi=x.hi;lo=x.lo;}
template<typename ty>
il operator ty () {return lo;}
il void operator += (__uint256_t x){lo+=x.lo;hi+=x.hi;if(lo<x.lo)hi++;}
il void operator -= (__uint256_t x){if(lo<x.lo)hi--;lo-=x.lo;hi-=x.hi;}
il __uint256_t operator + (__uint256_t x){__uint256_t y=*this;y+=x;return y;}
il __uint256_t operator - (__uint256_t x){__uint256_t y=*this;y-=x;return y;}
il void operator *= (__uint256_t x)
{
__uint128_t _x[4];
_x[0]=(lo&-1ull)*(x.lo&-1ull);
_x[1]=(lo&-1ull)*(x.lo>>64)+(lo>>64)*(x.lo&-1ull);
_x[2]=(lo>>64)*(x.lo>>64)+(hi&-1ull)*(x.lo&-1ull)+(x.hi&-1ull)*(lo&-1ull);
_x[3]=(hi>>64)*(x.lo&-1ull)+(x.hi>>64)*(lo&-1ull)+(lo>>64)*(x.hi&-1ull)+(x.lo>>64)*(hi&-1ull);
lo=_x[0]+(_x[1]<<64);
hi=(_x[1]>>64)+_x[2]+(_x[3]<<64);
}
il __uint256_t operator * (__uint256_t x){__uint256_t y=*this;y*=x;return y;}
il void operator |= (__uint256_t x){lo|=x.lo;hi|=x.hi;}
il __uint256_t operator | (__uint256_t x){return {hi|x.hi,lo|x.lo};}
il void operator ^= (__uint256_t x){lo^=x.lo;hi^=x.hi;}
il __uint256_t operator ^ (__uint256_t x){return {hi^x.hi,lo^x.lo};}
il void operator &= (__uint256_t x){lo&=x.lo;hi&=x.hi;}
il __uint256_t operator & (__uint256_t x){return {hi&x.hi,lo&x.lo};}
il void operator <<= (unsigned x){hi=(hi<<x)^(lo>>(128-x));lo<<=x;}
il __uint256_t operator << (unsigned x){return {(hi<<x)^(lo>>(128-x)),lo<<x};}
il void operator >>= (unsigned x){lo=(lo>>x)^(hi<<(128-x));hi>>=x;}
il __uint256_t operator >> (unsigned x){return {hi>>x,(lo>>x)^(hi<<(128-x))};}
il __uint256_t operator ~ () {return {~hi,~lo};}
il bool operator < (__uint256_t x){return hi!=x.hi?hi<x.hi:lo<x.lo;}
il bool operator <=(__uint256_t x){return !(hi!=x.hi?hi>x.hi:lo>x.lo);}
il bool operator > (__uint256_t x){return hi!=x.hi?hi>x.hi:lo>x.lo;}
il bool operator >=(__uint256_t x){return !(hi!=x.hi?hi<x.hi:lo<x.lo);}
il operator bool () {return hi||lo;};
il void operator /= (__uint256_t x)
{
__uint256_t ans=0,fl=~__uint256_t(0);
unsigned cnt=256;
while(cnt--)
{
fl>>=1;
if(x&fl) continue;
if(x<<cnt<=*this) *this-=(x<<cnt),ans|=__uint256_t(1)<<cnt;
}
*this=ans;
}
il __uint256_t operator / (__uint256_t x)
{
__uint256_t ans=0,fl=~__uint256_t(0);
unsigned cnt=256;
while(cnt--)
{
fl<<=1;
if(x&fl) continue;
if(x<<cnt<=*this) *this-=(x<<cnt),ans|=__uint256_t(1)<<cnt;
}
return ans;
}
#undef il
};这样我们就解决了取模的问题。
观察到下面将数字转换为字符时会有大量运算，因此我们可以将 __uint128_t 先拆分成 unsigned long long 进行运算，以减小 128 位整数运算的大常数。
下面是多模合并的代码。
for(int i=0;i<lim;i++)
{
__uint256_t modans=((__uint128_t)((int)num1_1[i])*inv1*mod23+
(__uint128_t)((int)num1_2[i])*inv2*mod13+(__uint128_t)((int)num1_3[i])*inv3*mod12);
modans=modans-__uint256_t((modans*brt).hi)*mod123;//取模
while(modans>=mod123) modans-=mod123;//取模
num+=__uint128_t(modans);
unsigned long long fl1=num/100000000,fl=num-fl1*100000000;num=fl1;//将 num 末 8 位取出到 ull 运算而不是直接用 num 运算
for(int k=0;k<sz;k++)
pr1[--len]=fl%10^48,fl/=10;
if(i>(lena+lenb+sz)/sz&&num==0) break;
}于是我们就得到了一个较快的高精度乘法模板，两个
实际上，该代码仍有很大的卡常空间，例如为了保证通用性多项式板子有很多无用的特判，删除对高精乘没用的部分后还可以快不少。另外，我使用的 ntt 模板也并非最优解，还存在比我的快很多的板子。