设计外贸英文网站,简述网站开发的流程,男女做那个视频的网站,下载百度导航app更好的阅读体验
任意模数多项式乘法 前言#xff1a; 在教练讲的时候脑子并不清醒#xff0c;所以没听懂。后来自己看博客学会了#xff0c;但目前只学了一种方法#xff1a;可拆系数FFT。为了方便日后复习#xff0c;决定先写下这个的笔记#xff0c;关于三模数NTT下次…更好的阅读体验
任意模数多项式乘法 前言 在教练讲的时候脑子并不清醒所以没听懂。后来自己看博客学会了但目前只学了一种方法可拆系数FFT。为了方便日后复习决定先写下这个的笔记关于三模数NTT下次再补。 建议准备好演算纸和笔本篇含有大量推算部分。 为什么不直接使用NTT/FFT
此处的模数是任意的所以我们使用NTT时有局限性。只有当模数满足下列情况时才可使用NTT。 模数为 P P P P r × 2 k 1 Pr\times 2^k 1 Pr×2k1其中 k k k足够大满足 2 k ≥ N 2^k \geq N 2k≥N其中 N N N为多项式扩充后的项数多项式乘法都需要将项数扩充到 2 2 2的幂。 如果不满足上述条件就不能直接使用NTT。
那为什么不使用FFT带模运算 首先不考虑模数的情况下多项式结果的系数不能超过 d o u b l e double double能表示的精度一般在 1 0 16 10^{16} 1016。超过后 d o u b l e double double所表示的结果将不再精确。
那怎么办目前可以给出两种解决方法
可拆系数FFT三模数NTTs可是我还没学懂QAQ/s
所以向可拆系数FFT进发
可拆系数FFT
怎么用可拆系数FFT
首先我们还是先假设一个情景 求 c 1 c_1 c1与 c 2 c_2 c2的卷积 我们可以把一个数拆成 a × 2 15 b a\times 2^{15}b a×215b的形式不一定是 2 15 2^{15} 215大概在 值域 \sqrt{值域} 值域 内。 c 1 a 1 × 2 15 a 2 c 2 b 1 × 2 15 b 2 c_1a_1\times 2^{15}a_2 \\ c_2b_1\times 2^{15}b_2 c1a1×215a2c2b1×215b2
那这俩的积为 c 1 × c 2 ( a 1 × 2 15 a 2 ) × ( b 1 × 2 15 b 2 ) a 1 b 1 × 2 30 ( a 1 b 2 a 2 b 1 ) × 2 15 a 2 b 2 \begin{aligned} c_1\times c_2(a_1\times 2^{15}a_2)\times(b_1\times 2^{15}b_2)\\a_1b_1\times 2^{30}(a_1b_2a_2b_1)\times 2^{15}a_2b_2 \end{aligned} c1×c2(a1×215a2)×(b1×215b2)a1b1×230(a1b2a2b1)×215a2b2
好耶那我们可以直接做4次FFT a 1 b 1 a_1b_1 a1b1、 a 1 b 2 a_1b_2 a1b2、 a 2 b 1 a_2b_1 a2b1、 a 2 b 2 a_2b_2 a2b2 s然后你发现正逆变换总共做了8次常数爆炸然后炸了/s
所以我们需要优化
优化可拆系数FFT
注意我们这里的优化会用到复数你可能会害怕得逃走但是你无需害怕因为s我也是这样过来的/s我会简单地讲一讲。 小资料可能不太学术规范 啥是复数 其实复数是一种含实部和虚部的一种数。我们知道 i − 1 i\sqrt{-1} i−1 i i i就是虚数。那我们以实部建立 x x x轴、虚部建立 y y y轴。 那我们假设在这个平面直角坐标系上有一个点 A ( 2 , 3 ) A(2,3) A(2,3)那这个点的复数表示为 2 3 i 23i 23i。 那你就基本知道什么是复数了让我们学一下基本运算吧。这个自己记一记好了 我们可以利用一下复数的乘法运算 ( a b i ) ( c d i ) ( a c − b d ) ( a d b c ) i ( a − b i ) ( c d i ) ( a c b d ) ( a d − b c ) i (abi)(cdi)(ac-bd)(adbc)i\\ (a-bi)(cdi)(acbd)(ad-bc)i (abi)(cdi)(ac−bd)(adbc)i(a−bi)(cdi)(acbd)(ad−bc)i
现在令 P 1 a 1 a 2 i P_1a_1a_2i P1a1a2i、 P 2 a 1 − a 2 i P_2a_1-a_2i P2a1−a2i、 Q b 1 b 2 i Qb_1b_2i Qb1b2i。 计算 P 1 Q P 2 Q ( a 1 b 1 − a 2 b 2 ) ( a 1 b 2 a 2 b 1 ) i ( a 1 b 1 a 2 b 2 ) ( a 1 b 2 − a 2 b 1 ) i 2 ( a 1 b 1 a 1 b 2 i ) \begin{aligned} P_1QP_2Q(a_1b_1-a_2b_2)(a_1b_2a_2b_1)i(a_1b_1a_2b_2)(a_1b_2-a_2b_1)i\\ 2(a_1b_1a_1b_2i) \end{aligned} P1QP2Q(a1b1−a2b2)(a1b2a2b1)i(a1b1a2b2)(a1b2−a2b1)i2(a1b1a1b2i)
如果我们将上式除以 2 2 2那我们可以得到 a 1 b 1 , a 1 b 2 a_1b_1,a_1b_2 a1b1,a1b2分别通过“实部相加/2”、“虚部相加/2”可得。
我们把得到的 a 1 b 1 a_1b_1 a1b1代入 P 2 Q P_2Q P2Q的实部可得 a 2 b 2 a_2b_2 a2b2。 类似地将 a 1 b 2 a_1b_2 a1b2代入 P 1 Q P_1Q P1Q的虚部可得 a 2 b 1 a_2b_1 a2b1。
注这里的运算请自己用演算纸推一下
我们就可以带回 c 1 × c 2 c_1\times c_2 c1×c2了。
那我们的任务就完成啦
代码实现
注此处的FFT部分与FFT版题的部分是一模一样的可参照本人之前所写的FFT笔记。
板子P4245 【模板】任意模数多项式乘法
code:
#includebits/stdc.h
using namespace std;#define ll long long
#define rp(i,o,p) for(ll io;ip;i)
#define pr(i,o,p) for(ll io;ip;--i)
#define cp complexlong doubleconst ll MAXN1e55,P1e97;
const ll lim32000; // lim sqrt(值域) - 1e9
const long double PIacos(-1.0);cp p1[MAXN2],p2[MAXN2],q[MAXN2];
ll n,m,limit;
cp p[MAXN2],omg[MAXN2],inv[MAXN2];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);}limit1;while(limitnm) limit1;init();fft(p1,omg),fft(p2,omg),fft(q,omg);rp(i,0,limit-1){long double xx,xy;xxp1[i].real()*q[i].real(),yyp1[i].imag()*q[i].imag();xyp1[i].real()*q[i].imag(),yxp1[i].imag()*q[i].real();p1[i]cp(xx/limit-yy/limit,xy/limityx/limit);xxp2[i].real()*q[i].real(),yyp2[i].imag()*q[i].imag();xyp2[i].real()*q[i].imag(),yxp2[i].imag()*q[i].real();p2[i]cp(xx/limit-yy/limit,xy/limityx/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,nm){ll a1b1(ll)((p1[i].real()p2[i].real())/20.5)%P;ll a1b2(ll)((p1[i].imag()p2[i].imag())/20.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(a1b2a2b1)*lima2b2)%P;ans(ans%PP)%P;printf(%lld ,ans);}return 0;
}后记
三模数NTT的内容会尽快补上的。
