拆系数FFT(任意模数FFT)
拆系数FFT
对于任意模数 mod
设m=√mod
把多项式A(x)和B(x)的系数都拆成a×m+b的形式,时a,b都小于m
提出,那么一个多项式就可以拆成两个多项式的加法
一个是a∗m的,一个是b的
直接乘法分配律,aa一遍,ab一遍,ba,bb一遍,四遍FFT
乘出来不会超过取模范围
然后合并直接
(a×m+b)(c×m+d)=a×c×m2+(a×c+b×d)m+b×d
这样子要进行 7 次 DFT
如果研究一下 myy 2016 年的集训队论文就会发现有 2 次 或者 1.5 次 DFT 的 FFT 算法
2次的够了吧
myy 巧妙的运用了复数的虚部,优化了算法
具体来说
设
C(x)=A(x)+iB(x)
D(x)=A(x)−iB(x)
假设
c(wkn) 表示将 C(x) DFT 后的点值
d(wkn) 表示将 D(x) DFT 后的点值
w 为 n 次单位复数根
设 conj(x) 表示 x 的共轭复数
那么
c(wk2n)=A(wk2n)+iB(wk2n)=2n−1∑j=0Ajwjk2n+iBjwjk2n
=2n−1∑j=0(Aj+iBj)(cosπjkn+isinπjkn)
d(wk2n)=A(wk2n)−iB(wk2n)=2n−1∑j=0(Aj−iBj)(cosπjkn+isinπjkn)
设 x=πjkn
那么
d(wk2n)=2n−1∑j=0(Ajcosx+Bjsinx)+i(Ajsinx−Bjcosx)
=conj(2n−1∑j=0(Ajcosx+Bjsinx)−i(Ajsinx−Bjcosx))
=conj(2n−1∑j=0(Ajcos(−x)−Bjsin(−x))+i(Ajsin(−x)+Bjcos(−x)))
=conj(2n−1∑j=0(Aj+iBj)(cos(−x)+isin(−x)))
=conj(2n−1∑j=0(Aj+iBj)(cos(2π−x)+isin(2π−x)))
=conj(c(w2n−k2n))
也就是说
只需要一次 DFT 求出 c 就可以求出 d
那么
DFT(A(k))=c(wk)+d(wk)2
DFT(B(k))=c(wk)−d(wk)2i
再用一次 DFT−1 还原出多项式,多项式乘法只要两次 DFT
考虑拆系数 FFT
DFT 两两合并,从 4 次变成 2 次
DFT−1 合并其中一个,从 3 次变成 2 次
一共 4 次
# include <bits/stdc++.h>
using namespace std;
typedef long long ll;
namespace IO {
const int maxn((1 << 21) + 1);
char ibuf[maxn], *iS, *iT, obuf[maxn], *oS = obuf, *oT = obuf + maxn - 1, c, st[65];
int f, tp;
char Getc() {
return (iS == iT ? (iT = (iS = ibuf) + fread(ibuf, 1, maxn, stdin), (iS == iT ? EOF : *iS++)) : *iS++);
}
void Flush() {
fwrite(obuf, 1, oS - obuf, stdout);
oS = obuf;
}
void Putc(char x) {
*oS++ = x;
if (oS == oT) Flush();
}
template <class Int> void In(Int &x) {
for (f = 1, c = Getc(); c < '0' || c > '9'; c = Getc()) f = c == '-' ? -1 : 1;
for (x = 0; c <= '9' && c >= '0'; c = Getc()) x = (x << 3) + (x << 1) + (c ^ 48);
x *= f;
}
template <class Int> void Out(Int x) {
if (!x) Putc('0');
if (x < 0) Putc('-'), x = -x;
while (x) st[++tp] = x % 10 + '0', x /= 10;
while (tp) Putc(st[tp--]);
}
}
using IO :: In;
using IO :: Out;
using IO :: Putc;
using IO :: Flush;
const int maxn(1 << 18);
const double pi(acos(-1));
struct Complex {
double a, b;
inline Complex() {
a = b = 0;
}
inline Complex(double _a, double _b) {
a = _a, b = _b;
}
inline Complex operator +(Complex x) const {
return Complex(a + x.a, b + x.b);
}
inline Complex operator -(Complex x) const {
return Complex(a - x.a, b - x.b);
}
inline Complex operator *(Complex x) const {
return Complex(a * x.a - b * x.b, a * x.b + b * x.a);
}
inline Complex Conj() {
return Complex(a, -b);
}
};
Complex a[maxn], b[maxn], w[maxn], a1[maxn], a2[maxn];
int r[maxn], l, deg, g[maxn], h[maxn], mod;
inline void FFT(Complex *p, int opt) {
register int i, j, k, t;
register Complex wn, x, y;
for (i = 0; i < deg; ++i) if (r[i] < i) swap(p[r[i]], p[i]);
for (i = 1; i < deg; i <<= 1)
for(t = i << 1, j = 0; j < deg; j += t)
for (k = 0; k < i; ++k) {
wn = w[deg / i * k];
if (opt == -1) wn.b *= -1;
x = p[j + k], y = wn * p[i + j + k];
p[j + k] = x + y, p[i + j + k] = x - y;
}
}
inline void Mul(int n, int *p, int *q, int *f) {
register int i, k, v1, v2, v3;
register Complex ca, cb, da1, da2, db1, db2;
for (deg = 1, l = 0; deg < n; deg <<= 1) ++l;
for (i = 0; i < deg; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
for (i = 0; i < deg; ++i) w[i] = Complex(cos(pi * i / deg), sin(pi * i / deg));
for (i = 0; i < n; ++i) a[i] = Complex(p[i] & 32767, p[i] >> 15), b[i] = Complex(q[i] & 32767, q[i] >> 15);
for (FFT(a, 1), FFT(b, 1), i = 0; i < deg; ++i) {
k = (deg - i) & (deg - 1), ca = a[k].Conj(), cb = b[k].Conj();
da1 = (ca + a[i]) * Complex(0.5, 0), da2 = (a[i] - ca) * Complex(0, -0.5);
db1 = (cb + b[i]) * Complex(0.5, 0), db2 = (b[i] - cb) * Complex(0, -0.5);
a1[i] = da1 * db1 + (da1 * db2 + da2 * db1) * Complex(0, 1), a2[i] = da2 * db2;
}
for (FFT(a1, -1), FFT(a2, -1), i = 0; i < deg; ++i) {
v1 = (ll)(a1[i].a / deg + 0.5) % mod, v2 = (ll)(a1[i].b / deg + 0.5) % mod;
v3 = (ll)(a2[i].a / deg + 0.5) % mod, f[i] = (((ll)v3 << 30) + ((ll)v2 << 15) + v1) % mod;
if (f[i] < 0) f[i] += mod;
}
}
int main() {
register int len, i, n, m;
In(n), In(m), In(mod), ++n, ++m;
for (i = 0; i < n; ++i) In(h[i]), h[i] %= mod;
for (i = 0; i < m; ++i) In(g[i]), g[i] %= mod;
for (len = 1, n += m - 1; len < n; len <<= 1);
for (Mul(len, h, g, g), i = 0; i < n; ++i) Out(g[i]), Putc(' ');
return Flush(), 0;
}
分类:
数论数学--FFT\NTT
, A--模板\算法\知识点总结
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· [.NET]调用本地 Deepseek 模型
· 一个费力不讨好的项目,让我损失了近一半的绩效!
· .NET Core 托管堆内存泄露/CPU异常的常见思路
· PostgreSQL 和 SQL Server 在统计信息维护中的关键差异
· C++代码改造为UTF-8编码问题的总结
· 【.NET】调用本地 Deepseek 模型
· CSnakes vs Python.NET:高效嵌入与灵活互通的跨语言方案对比
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· Plotly.NET 一个为 .NET 打造的强大开源交互式图表库