隐藏页面特效

2333

多项式全集之二 任长任模的FFT:

三模NTT实现任模FFT:

1.为什么要用MTT:当p不是NTT模数或者多项式长度大于模数限制时,就要使用MTT。

2.MTT的使用原理:我们对初始多项式取模,那么如果在不取模卷积情况下,答案x不会超过N×p2。我们取三个NTT模数p1,p2,p3,分别做多项式乘法,得到x分别mod p1,p2,p3的答案,通过CRT合并可以得到x mod p1p2p3的答案,如果x<p1p2p3那么就可以得到准确的答案,再对p取模即可。

3.CRT合并的小优化:

step 0:初始式子

{xc1(mod p1)xc2(mod p2)xc3(mod p3)

step 1:把一式二式合并(LL范围内)。

{xa(mod p1p2)xc3(mod p3)

step 2:再次合并(不需要long double 快速乘)。

4.常用NTT模数:

以下模数的共同g=3189

p=r×2k+1 k g
104857601 22 3
167772161 25 3
469762049 26 3
950009857 21 7
998244353 23 3
1004535809 21 3
2013265921 27 31
2281701377 27 3
3221225473 30 5
拆系数FFT(CFFT)实现任模FFT:

1.实现原理:运用实数FFT不取模做乘法,然后取模回归到整数。但是由于误差较大(值域是1023),我们令t=m把系数ai=kit+bi,对ki,ti交叉做四遍卷积,求出答案按系数贡献取模加入。

2.可按合并DFT的方法优化DFT次数。

bluestein算法实现任长FFT:

m不是2的幂次的时候,我们从式子入手:
Xi=aiwmi22,Yi=wmi22

Ak=j=0m1ajwmjk=j=0m1ajwmj2+k2(kj)22=wmk22j=0m1ajwmj22wm(kj)22=wmk22j=0m1XjYkj

喜闻乐见的模板:

三模NTT模板(注意:不可以MTT回来,因为系数会取模)

namespace MTT{ typedef long long LL; int n, m; LL p, mod; const LL p1 = 998244353; const LL p2 = 1004535809; const LL p3 = 104857601; const int g = 3189; LL a[300005], b[300005], c[300005], cpa[300005], cpb[300005]; LL c3[300005], c1[300005], c2[300005]; LL qpow(LL a, LL b, LL mod) { LL ans = 1; while(b) { if(b & 1) ans = ans * a % mod; a = a * a % mod; b >>= 1; } return ans; } const LL inv12 = qpow(p1, p2 - 2, p2); const LL inv123 = qpow(p1 * p2 % p3, p3 - 2, p3); struct p_l_e{ int wz[300005]; void MTT(LL *a, int N, int op) { for(int i = 0; i < N; i++) if(i < wz[i]) swap(a[i], a[wz[i]]); for(int le = 2; le <= N; le <<= 1) { int mid = le >> 1; LL wn = qpow(g, (mod - 1) / le, mod); if(op == -1) wn = qpow(wn, mod - 2, mod); for(int i = 0; i < N ;i += le) { LL w = 1, x, y; for(int j = 0; j < mid; j++) { x = a[i + j]; y = a[i + j + mid] * w % mod; a[i + j] = (x + y) % mod; a[i + j + mid] = (x - y + mod) % mod; w = w * wn % mod; } } } } void mult(LL *a, LL *b, LL *c, int M) { int N = 1, len = 0; while(N < M) N <<= 1, len++; for(int i = 0; i < N; i++) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1)); MTT(a, N, 1); MTT(b, N, 1); for(int i = 0; i < N; i++) c[i] = a[i] * b[i] % mod; MTT(c, N, -1); LL t = qpow(N, mod - 2, mod); for(int i = 0; i < N; i++) c[i] = c[i] * t % mod; } }PLE; LL CRT(LL c1, LL c2, LL c3) { LL x = (c1 + p1 * ((c2 - c1 + p2) % p2 * inv12 % p2)); LL y = (x % p + p1 * p2 % p * ((c3 - x % p3 + p3) % p3 * inv123 % p3) % p) % p; return y; } void merge(LL *c1, LL *c2, LL *c3, LL *c, int N) { for(int i = 0; i < N; i++) c[i] = CRT(c1[i], c2[i], c3[i]); return; } void main() { scanf("%d%d%lld", &n, &m, &p); n++; m++; for(int i = 0; i < n; i++) scanf("%lld", &a[i]); for(int i = 0; i < m; i++) scanf("%lld", &b[i]); mod = p1; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c1, n + m - 1); mod = p2; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c2, n + m - 1); mod = p3; memcpy(cpa, a, sizeof(a)); memcpy(cpb, b, sizeof(b)); PLE.mult(cpa, cpb, c3, n + m - 1); merge(c1, c2, c3, c, n + m - 1); for(int i = 0; i < n + m - 1; i++) printf("%lld ", (c[i] % p + p) % p); return; } }

拆系数FFT模板(注意:相同系数的两项可以合并一起IDFT。采用共轭优化法,只进行四次DFT)

namespace CFFT{ typedef long long LL; int n, m, p ,sqrp; int a[300005], b[300005]; const long double pi = acos(-1); struct cp{ long double x, y; cp() {x = y = 0;} cp(long double X,long double Y) {x = X; y = Y; } cp conj() {return (cp) {x, -y};} }ka[300005], kb[300005], ta[300005], tb[300005], kk[300005], kt[300005], tt[300005], c[300005], I(0, 1), d[300005]; cp operator+ (const cp &a, const cp &b) {return (cp){a.x + b.x, a.y + b.y}; } cp operator- (const cp &a, const cp &b) {return (cp){a.x - b.x, a.y - b.y}; } cp operator* (const cp &a, const cp &b) {return (cp){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x};} cp operator* (const cp &a, long double b) {return (cp){a.x * b, a.y * b};} cp operator/ (const cp &a, long double b) {return (cp){a.x / b, a.y / b};} struct p_l_e{ int wz[300005]; void FFT(cp *a, int N, int op){ for(int i = 0; i < N; i++) if(i < wz[i]) swap(a[i], a[wz[i]]); for(int le = 2; le <= N; le <<= 1){ int mid = le >> 1; cp x, y, w, wn = (cp){cos(op * 2 * pi / le), sin(op * 2 * pi / le)}; for(int i = 0; i < N; i += le){ w = (cp){1, 0}; for(int j = 0; j < mid; j++){ x = a[i + j]; y = a[i + j + mid] * w; a[i + j] = x + y; a[i + j + mid] = x - y; w = w * wn; } } } } void D_FFT(cp *a, cp *b, int N, int op){ for(int i = 0; i < N; i++) d[i] = a[i] + I * b[i]; FFT(d, N, op); d[N] = d[0]; if(op == 1){ for(int i = 0; i < N; i++){ a[i] = (d[i] + d[N - i].conj()) / 2; b[i] = I * (-1) * (d[i] - d[N - i].conj()) / 2; } } else { for(int i = 0; i < N; i++){ a[i] = cp(d[i].x, 0); b[i] = cp(d[i].y, 0); } } d[N] = cp(0, 0); } void mult(int *a, int *b, int M){ int N = 1, len = 0; while(N < M) N <<= 1, len++; for(int i = 0; i < N; i++) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1)); for(int i = 0; i < N; i++){ ka[i].x = a[i] >> 15; kb[i].x = b[i] >> 15; ta[i].x = a[i] & 32767; tb[i].x = b[i] & 32767; } D_FFT(ta, ka, N, 1); D_FFT(tb, kb, N, 1); for(int i = 0; i < N; i++){ kk[i] = ka[i] * kb[i]; kt[i] = ka[i] * tb[i] + ta[i] * kb[i]; tt[i] = ta[i] * tb[i]; } D_FFT(tt, kk, N, -1); FFT(kt, N, -1); for(int i = 0; i < N; i++){ tt[i] = tt[i] / N; kt[i] = kt[i] / N; kk[i] = kk[i] / N; } } }PLE; void main() { scanf("%d%d%d", &n, &m, &p); n++; m++; for(int i = 0; i < n; i++) scanf("%d", &a[i]),a[i] = a[i] % p; for(int i = 0; i < m; i++) scanf("%d", &b[i]),b[i] = b[i] % p; PLE.mult(a, b, n + m - 1); for(int i = 0; i < n + m - 1; i++) printf("%lld ",(((((LL)round(kk[i].x)) % p) << 30) + ((((LL)round(kt[i].x)) % p) << 15) + ((LL)round(tt[i].x)) % p) % p); } }

blue_stein模板:

struct polynie { CP getw(int m, int k, int op) { return CP(cos(2 * pi * k / m), op * sin(2 * pi * k / m)); } int wz[MAXN]; CP A[MAXN], B[MAXN], C[MAXN]; void FFT(CP *a, int N, int op) { rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]); for(int l = 2; l <= N; l <<= 1) { int mid = l >> 1; CP x, y, w, wn = CP(cos(pi / mid), sin(op * pi / mid)); for(int i = 0; i < N; i += l) { w = CP(1, 0); rop(j, 0, mid) { x = a[i + j]; y = w * a[i + j + mid]; a[i + j] = x + y; a[i + j + mid] = x - y; w = w * wn; } } } } void mult(CP *a, CP *b, CP *c, int M) { int N = 1, len = 0; while(N < M) N <<= 1, len++; rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1)); FFT(a, N, 1); FFT(b, N, 1); rop(i, 0, N) c[i] = a[i] * b[i]; FFT(c, N, -1); rop(i, 0, N) c[i].x = c[i].x / N, c[i].y = c[i].y / N; } void blue_stein(CP *a, int M, int op) { int M2 = M << 1; memset(A, 0, sizeof(A)); memset(B, 0, sizeof(B)); memset(C, 0, sizeof(C)); rop(i, 0, M) A[i] = a[i] * getw(M2, 1ll * i * i % M2, op); rop(i, 0, M2) B[i] = getw(M2, 1ll * (i - M) * (i - M) % M2, -op); mult(A, B, C, M2 + M - 1); rop(i, 0, M) a[i] = C[i + M] * getw(M2, 1ll * i * i % M2, op); if(op == -1) rop(i, 0, M) a[i].x = a[i].x / M, a[i].y = a[i].y / M; } }PLE;

多项式全集之三 多项式求逆与除法:

多项式求逆:

1.问题描述:

已知F(x),且F(x)G(x)1(mod xn),求G(x)

2.推导过程:

\begin{align}
B(x)&\equiv F(x){-1}&(mod~x)\
由于
G(x)&\equiv F(x)^{-1}& (mod~x^n)\
所以
G(x)&\equiv F(x)^{-1}& (mod~x^{\lceil \frac n 2 \rceil})\
G(x)& \equiv B(x)&(mod~x^{\lceil \frac n 2 \rceil})\
G(x)-B(x)& \equiv 0&(mod~x^{\lceil \frac n 2 \rceil})\
\end{align
}
两边平方,得:

由于[G(x)B(x)]2的第k<n项为

i=0k[gixibixi][gkixkibkixki]

i,ki一定有一项<n2,所以
\begin{align}
[G(x)-B(x)]^2&\equiv 0 &(mod~x^n)\
G2(x)+B2(x)-2G(x)B(x)&\equiv 0&(mod~x^n)
\end{align
}
两边同乘A(x),得:
\begin{align}
A(x)G2(x)+A(x)B2(x)-2A(x)G(x)B(x)&\equiv 0& (mod~x^n)\
G(x)+A(x)B^2(x)-2B(x)&\equiv 0& (mod~x^n)\
G(x)&\equiv 2B(x)-A(x)B2(x)&(mod~xn)\
G(x)&\equiv B(x)[2-A(x)B(x)]&(mod~x^n)\
\end{align
}

多项式除法:

1.问题描述:

已知一个n次多项式A(x),一个m次多项式B(x),且A(x)=B(x)C(x)+D(x),求nm次多项式C(x)<m次多项式D(x)

2.推导过程:

A(x)=B(x)C(x)+D(x)得:
\begin{align}
A(\frac 1 x)&=B(\frac 1 x)C(\frac 1 x)+D(\frac 1 x)\
x^nA(\frac 1 x)&=x^nB(\frac 1 x)C(\frac 1 x)+x^nD(\frac 1 x)\
x^nA(\frac 1 x)&=x^mB(\frac 1 x)x^{n-m}C(\frac 1 x)+x{m-1}xD(\frac 1 x)\
A_r(x)&=B_r(x)C_r(x)+x^{n-m+1}D(x)\
A_r(x)&=B_r(x)C_r(x)&(mod ;x^{n-m+1})\
B_r(x)&=A_r^{-1}(x)C_r(x)&(mod ;x^{n-m+1})\
\end{align
}
求逆可得Br(x),再反转得B(x),然后乘C(x)去减A(x)D(x).

喜闻乐见的模板:
namespace INV{ typedef long long LL; int n, a[300005], b[300005]; const int mod = 998244353; const int g = 3189; int qpow(int a, int b){ int ans = 1; while(b){ if(b & 1) ans = 1ll * ans * a % mod; a = 1ll * a * a % mod; b >>= 1; } return ans; } struct p_l_e{ int wz[300005], i_c[300005]; void NTT(int *a, int N, int op){ for(int i = 0; i < N; i++) if(i < wz[i]) swap(a[i], a[wz[i]]); for(int le = 2; le <= N; le <<= 1){ int mid = le >> 1, wn = qpow(g, (mod - 1) / le); if(op == -1) wn = qpow(wn, mod - 2); for(int i = 0; i < N; i += le){ LL w = 1; int x, y; for(int j = 0; j < mid; j++){ x = a[i + j]; y = w * a[i + j + mid] % mod; a[i + j] = (x + y) % mod; a[i + j + mid] = (x - y + mod) % mod; w = w * wn % mod; } } } } int init(int M){ int N = 1, len = 0; while(N < M) N <<= 1, len++; for(int i = 0; i < N; i++) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1)); return N; } void INV(int *a, int *b, int deg){ if(deg == 1){b[0] = qpow(a[0], mod - 2); return;} INV(a, b, (deg + 1) >> 1); int N = init(deg + deg - 1); for(int i = 0; i < deg; i++) i_c[i] = a[i]; for(int i = deg; i < N; i++) i_c[i] = 0; NTT(b, N, 1);NTT(i_c, N, 1); for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * (2 - 1ll * b[i] * i_c[i] % mod + mod) % mod; NTT(b, N, -1); int t = qpow(N, mod - 2); for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod; for(int i = deg; i < N; i++) b[i] = 0; } }PLE; void main(){ scanf("%d", &n); for(int i = 0; i < n; i++) scanf("%d", &a[i]); PLE.INV(a, b, n); for(int i = 0; i < n; i++) printf("%d ",b[i]); } }

多项式全集之四 多项式ln与exp:

多项式ln:

1.做法:

G(x)=lnF(x)

两边求导得

G(x)=F(x)F(x)

积分回去即可。

2.应用:

eF(x)=k0Fk(x)k!=G(x)

F(x)=lnG(x)

这个的组合意义是:无序组合。

F(x)fi表示一些东西,那么这些东西有序组合的方案数为

F0(x)+F1(x)+F2(x)+=11F(x)

而无序组成的方案数为:

F0(x)0!+F1(x)1!+F2(x)2!+=eF(x)

如果无序组合方案数好求,那么求ln就能得到F(x)

例题

多项式exp
喜闻乐见的代码:

多项式ln:

namespace PLE_ln{ struct polyme { int li[SZ], wz[SZ]; void NTT(int *a, int N, int op) { rop(i, 0, N) if(i < wz[i]) swap(a[i], a[wz[i]]); for(int l = 2; l <= N; l <<= 1) { int mid = l >> 1; int x, y, w, wn = qpow(g, (mod - 1) / l); if(op) wn = qpow(wn, mod - 2); for(int i = 0; i < N; i += l) { w = 1; for(int j = 0; j < mid; ++j) { x = a[i + j]; y = 1ll * w * a[i + j + mid] % mod; a[i + j] = (x + y) % mod; a[i + j + mid] = (x - y + mod) % mod; w = 1ll * w * wn % mod; } } } } void qd(int *a, int *b, int n) { rop(i, 0, n) b[i] = 1ll * a[i + 1] * (i + 1) % mod; } void jf(int *a, int *b, int n) { rop(i, 1, n) b[i] = 1ll * a[i - 1] * qpow(i, mod - 2) % mod; } void mult(int *a, int *b, int *c, int M) { int N = 1, len = 0; while(N < M) N <<= 1, len ++; rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1)); NTT(a, N, 0); NTT(b, N, 0); rop(i, 0, N) c[i] = 1ll * a[i] * b[i] % mod; NTT(c, N, 1); int t = qpow(N, mod - 2); rop(i, 0, N) c[i] = 1ll * c[i] * t % mod; } void inv(int *a, int *b, int deg) { if(deg == 1) {b[0] = qpow(a[0], mod - 2) % mod; return;} inv(a, b, (deg + 1) >> 1); rop(i, 0, deg) li[i] = a[i]; int N = 1, len = 0; while(N < deg + deg - 1) N <<= 1, len ++; rop(i, 0, N) wz[i] = (wz[i >> 1] >> 1) | ((i & 1) << (len - 1)); rop(i, deg, N) li[i] = 0; NTT(li, N, 0); NTT(b, N, 0); rop(i, 0, N) b[i] = 1ll * b[i] * (2 - 1ll * li[i] * b[i] % mod + mod) % mod; NTT(b, N, 1); int t = qpow(N, mod - 2); for(int i = 0; i < N; i++) b[i] = 1ll * b[i] * t % mod; rop(i, deg, N) b[i] = 0; } }PLE; int a[SZ], da[SZ], ia[SZ], dla[SZ], la[SZ], n; void main() { scanf("%d", &n); rop(i, 0, n) scanf("%d", &a[i]); PLE.qd(a, da, n); PLE.inv(a, ia, n); PLE.mult(ia, da, dla, n + n - 1); PLE.jf(dla, la, n); rop(i, 0, n) printf("%d ", la[i]); } }

多项式全集之五 多项式快速幂与开根:

多项式快速幂方法一:

多项式全集之六 多项式快速插值与多点求值:

多项式全集之七 拉格朗日插值:


__EOF__

本文作者Smeow
本文链接https://www.cnblogs.com/Smeow/p/10729659.html
关于博主:评论和私信会在第一时间回复。或者直接私信我。
版权声明:本博客所有文章除特别声明外,均采用 BY-NC-SA 许可协议。转载请注明出处!
声援博主:如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。您的鼓励是博主的最大动力!
posted @   Smeow  阅读(382)  评论(0编辑  收藏  举报
编辑推荐:
· Linux glibc自带哈希表的用例及性能测试
· 深入理解 Mybatis 分库分表执行原理
· 如何打造一个高并发系统?
· .NET Core GC压缩(compact_phase)底层原理浅谈
· 现代计算机视觉入门之:什么是图片特征编码
阅读排行:
· 手把手教你在本地部署DeepSeek R1,搭建web-ui ,建议收藏!
· Spring AI + Ollama 实现 deepseek-r1 的API服务和调用
· 数据库服务器 SQL Server 版本升级公告
· 程序员常用高效实用工具推荐,办公效率提升利器!
· C#/.NET/.NET Core技术前沿周刊 | 第 23 期(2025年1.20-1.26)
点击右上角即可分享
微信分享提示