任意模数多项式乘法MTT(可拆系数FFT、三模数NTT)笔记
任意模数多项式乘法
前言:
在教练讲的时候脑子并不清醒,所以没听懂。后来自己看博客学会了,但目前只学了一种方法:可拆系数FFT。为了方便日后复习,决定先写下这个的笔记,关于三模数NTT下次再补。建议:准备好演算纸和笔,本篇含有大量推算部分。
注:本篇文章是本蒟写的,dalao随便看看就好,不必争论对错(但是欢迎指出文章存在的小错误)。
为什么不直接使用NTT/FFT
此处的模数是任意的,所以我们使用NTT时,有局限性。只有当模数满足下列情况时才可使用NTT。
模数为\(P\),\(P=r\times 2^k + 1\),其中\(k\)足够大,满足\(2^k \geq N\),其中\(N\)为多项式扩充后的项数,多项式乘法都需要将项数扩充到\(2\)的幂。
如果不满足上述条件,就不能直接使用NTT。
那为什么不使用FFT带模运算?
首先不考虑模数的情况下,多项式结果的系数不能超过\(double\)能表示的精度(一般在\(10^{16}\))。超过后,\(double\)所表示的结果将不再精确。
那怎么办?目前可以给出两种解决方法
- 可拆系数FFT
- 三模数NTT
(可是我还没学懂QAQ)
所以,向可拆系数FFT进发!
可拆系数FFT
怎么用可拆系数FFT
首先我们还是先假设一个情景:
求:\(c_1\)与\(c_2\)的卷积
我们可以把一个数拆成\(a\times 2^{15}+b\)的形式(不一定是\(2^{15}\),大概在\(\sqrt{值域}\)内)。
那这俩的积为
好耶!那我们可以直接做4次FFT(\(a_1b_1\)、\(a_1b_2\)、\(a_2b_1\)、\(a_2b_2\))!
然后你发现正逆变换总共做了8次常数爆炸然后炸了
所以,我们需要优化!
优化可拆系数FFT
注意:我们这里的优化会用到复数,你可能会害怕得逃走,但是你无需害怕,因为(我也是这样过来的)我会简单地讲一讲。
小资料(可能不太学术规范):
啥是复数?其实复数,是一种含实部和虚部的一种数。我们知道\(i=\sqrt{-1}\),\(i\)就是虚数。那我们以实部建立\(x\)轴、虚部建立\(y\)轴。
那我们假设在这个平面直角坐标系上有一个点\(A(2,3)\),那这个点的复数表示为\(2+3i\)。那你就基本知道什么是复数了,让我们学一下基本运算吧。(这个自己记一记好了)
我们可以利用一下复数的乘法运算:
现在令\(P_1=a_1+a_2i\)、\(P_2=a_1-a_2i\)、\(Q=b_1+b_2i\)。
计算:
如果我们将上式除以\(2\),那我们可以得到\(a_1b_1,a_1b_2\)(分别通过“实部相加/2”、“虚部相加/2”可得)。
我们把得到的\(a_1b_1\)代入\(P_2Q\)的实部,可得\(a_2b_2\)。
类似地,将\(a_1b_2\)代入\(P_1Q\)的虚部,可得\(a_2b_1\)。
(注:这里的运算请自己用演算纸推一下)
我们就可以带回\(c_1\times c_2\)了。
那我们的任务就完成啦!!
代码实现
注:此处的FFT部分与FFT版题的部分是一模一样的,可参照本人之前所写的FFT笔记。
code:
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define rp(i,o,p) for(ll i=o;i<=p;++i)
#define pr(i,o,p) for(ll i=o;i>=p;--i)
#define cp complex<long double>
const ll MAXN=1e5+5,P=1e9+7;
const ll lim=32000; // lim = sqrt(值域) <- 1e9
const long double PI=acos(-1.0);
cp p1[MAXN<<2],p2[MAXN<<2],q[MAXN<<2];
ll n,m,limit;
cp p[MAXN<<2],omg[MAXN<<2],inv[MAXN<<2];
void init() {
for (ll i = 0; i < limit; ++i) {
omg[i] = cp(cos(2 * PI * i / limit), sin(2 * PI * i / limit));
inv[i] = conj(omg[i]);
}
}
void fft(cp *a, cp *omg) {
for (ll i = 0, j = 0; i < limit; ++i) {
if (i > j)
swap(a[i], a[j]);
for (ll l = limit >> 1; (j ^= l) < l; l >>= 1)
;
}
for (ll l = 2; l <= limit; l <<= 1) {
ll m = l >> 1;
for (cp *p = a; p != a + limit; p += l) {
rp(i, 0, m - 1) {
cp t = omg[limit / l * i] * p[i + m];
p[i + m] = p[i] - t;
p[i] += t;
}
}
}
}
int main()
{
scanf("%lld%lld",&n,&m);
rp(i,0,n)
{
ll x;
scanf("%lld",&x);
p1[i]=cp(x/lim,x%lim);
p2[i]=cp(x/lim,-x%lim);
}
rp(i,0,m)
{
ll x;
scanf("%lld",&x);
q[i]=cp(x/lim,x%lim);
}
limit=1;
while(limit<=n+m) limit<<=1;
init();
fft(p1,omg),fft(p2,omg),fft(q,omg);
rp(i,0,limit-1)
{
long double xx,xy;
xx=p1[i].real()*q[i].real(),yy=p1[i].imag()*q[i].imag();
xy=p1[i].real()*q[i].imag(),yx=p1[i].imag()*q[i].real();
p1[i]=cp(xx/limit-yy/limit,xy/limit+yx/limit);
xx=p2[i].real()*q[i].real(),yy=p2[i].imag()*q[i].imag();
xy=p2[i].real()*q[i].imag(),yx=p2[i].imag()*q[i].real();
p2[i]=cp(xx/limit-yy/limit,xy/limit+yx/limit);
}
/*
上面的循环等价于这两行循环,但是因为c++中complex模板为const类型,
单独的real()或imag()值不可以直接修改,只可以两个同时赋值来修改,
所以上面的循环还用到了复数的除法运算,自己看看吧
rp(i,0,limit-1) q[i].real()/=limit,q[i].imag()/=limit;
rp(i,0,limit-1) p1[i]*=q[i],p2[i]*=q[i];
*/
fft(p1,inv),fft(p2,inv);
rp(i,0,n+m)
{
ll a1b1=(ll)((p1[i].real()+p2[i].real())/2+0.5)%P;
ll a1b2=(ll)((p1[i].imag()+p2[i].imag())/2+0.5)%P;
ll a2b1=(ll)((p1[i].imag()+0.5)-a1b2)%P;
ll a2b2=(ll)((p2[i].real()+0.5)-a1b1)%P;
ll ans=(a1b1*lim*lim+(a1b2+a2b1)*lim+a2b2)%P;
ans=(ans%P+P)%P;
printf("%lld ",ans);
}
return 0;
}
后记
三模数NTT的内容会尽快补上的。