上一篇：
（请无视此图）
我最近学了一个常数小还不用背三个模数的做法：拆系数法
（以下默认多项式项数 \(N=10^5\) ，系数不超过 \(M=10^9\) 且为非负整数）
我们放弃「数论变换」「利用原根性质」之类的想法，来点简单粗暴的：用实数 FFT 把原始结果算出来，然后直接取模。
为什么过去我们没有这样做呢？因为卷积结果的系数最大可能达到 \(NM^2=10^{24}\) ，long double 的精度也不够（通常情况下 double 的精度约为 15 位十进制，long double 的精度约为 19 位十进制）。考虑「拆系数」来牺牲时间保证精度。
设相乘的两个多项式为 \(A(x)=\sum a_ix^i\) 和 \(B(x)=\sum b_ix^i\) ，结果为 \(C(x)=\sum c_ix^i\) 。把\(A(x)\) 拆成两个多项式 \(A_0(x)=\sum {a_0}_ix^i\) 和 \(A_1(x)=\sum {a_1}_ix^i\) ，其中 \({a_1}_i=\lfloor\frac{a_i}{S}\rfloor\) ，\({a_0}_i=a_i-S\cdot {a_1}_i\) ，即 \(a_i=S\cdot {a_1}_i+{a_0}_i\) 。对 \(B(x)\) 也作同样的操作。这样，就有
\[\begin{aligned} a_ib_j&=(S\cdot {a_1}_i+{a_0}_i)(S\cdot {b_1}_j+{b_0}_j)\\ &={a_1}_i{b_1}_jS^2+({a_1}_i{b_0}_j+{a_0}_i{b_1}_j)S+{a_0}_i{b_0}_j \end{aligned}\]
于是直接计算 \(C_1=A_1*B_1\) ，\(C_2=A_1*B_0+A_0*B_1\) ，\(C_3=A_0*B_0\) ，然后 \(c_i={c_1}_iS^2+{c_2}_iS+{c_3}_i\) ，算最后一步的时候对「任意模数」取模即可。
这样，大致估计一下 \(C_1\) 的最大系数是 \(N\cdot(\frac{M}{S})^2=\frac{NM^2}{S^2}\) ，\(C_2\) 最大系数是 \(2N\cdot\frac{M}{S}\cdot S=2NM\) ，\(C_3\) 最大系数是 \(NS^2\) 。当 \(S=\sqrt{M}\) 时，以上三项均为 \(NM\) 级别，即 \(10^{14}\) 左右，足以保证精度（跟瓜学的一般偷懒直接 \(S=32768\) ）。
代码：
事实上由于只跟 7 个多项式有关，所以只需要进行 7 次 FFT 。我写的常数特别大了不要跟我学 ……
题目：
#include#include #include #include #include using namespace std;namespace zyt{ template inline bool read(T &x) { char c; bool f = false; x = 0; do c = getchar(); while (c != EOF && c != '-' && !isdigit(c)); if (c == EOF) return false; if (c == '-') f = true, c = getchar(); do x = x * 10 + c - '0', c = getchar(); while (isdigit(c)); if (f) x = -x; return true; } template inline void write(T x) { static char buf[20]; char *pos = buf; if (x < 0) putchar('-'), x = -x; do *pos++ = x % 10 + '0'; while (x /= 10); while (pos > buf) putchar(*--pos); } const int N = 1e5 + 10, S = 1 << 15, p = 1e9 + 7, B = 18; typedef long double ld; typedef long long ll; int power(int a, int b) { int ans = 1; while (b) { if (b & 1) ans = (ll)ans * a % p; a = (ll)a * a % p; b >>= 1; } return ans; } int get_inv(const int a) { return power(a, p - 2); } ll dtol(const ld x) { return ll(fabs(x) + 0.5) * (x < 0 ? -1 : 1); } namespace Polynomial { const int LEN = 1 << B; const ld PI = acos(-1.0L); struct cpx { ld x, y; cpx(const ld _x = 0, const ld _y = 0) : x(_x), y(_y) {} cpx conj() { return cpx(x, -y); } }omega[LEN], winv[LEN]; ll ctol(const cpx &a) { return dtol(a.x); } int rev[LEN]; cpx operator + (const cpx &a, const cpx &b) { return cpx(a.x + b.x, a.y + b.y); } cpx operator - (const cpx &a, const cpx &b) { return cpx(a.x - b.x, a.y - b.y); } cpx operator * (const cpx &a, const cpx &b) { return cpx(a.x * b.x - a.y * b.y, a.y * b.x + a.x * b.y); } void init(const int n, const int lg2) { cpx w = cpx(cos(2.0L * PI / n), sin(2.0L * PI / n)), wi = w.conj(); omega[0] = winv[0] = 1; for (int i = 1; i < n; i++) omega[i] = omega[i - 1] * w, winv[i] = winv[i - 1] * wi; for (int i = 0; i < n; i++) rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) << (lg2 - 1))); } void FFT(cpx *const a, const cpx * const w, const int n) { for (int i = 0; i < n; i++) if (i < rev[i]) swap(a[i], a[rev[i]]); for (int l = 1; l < n; l <<= 1) for (int i = 0; i < n; i += (l << 1)) for (int k = 0; k < l; k++) { cpx x = a[i + k], y = a[i + l + k] * w[n / (l << 1) * k]; a[i + k] = x + y, a[i + l + k] = x - y; } } void mul(const cpx *const a, const cpx *const b, cpx *const c, const int n) { static cpx x[LEN], y[LEN]; int m = 1, lg2 = 0; while (m < (n + n - 1)) m <<= 1, ++lg2; init(m, lg2); memcpy(x, a, sizeof(cpx[n])); memcpy(y, b, sizeof(cpx[n])); for (int i = n; i < m; i++) x[i] = y[i] = 0; FFT(x, omega, m), FFT(y, omega, m); for (int i = 0; i < m; i++) x[i] = x[i] * y[i]; FFT(x, winv, m); for (int i = 0; i < n; i++) c[i] = cpx(x[i].x / m, 0.0); } void MTT(const int *const a, const int *const b, int *const ans, const int n) { const int S = 1 << 15; static cpx a0[LEN], a1[LEN], b0[LEN], b1[LEN], c1[LEN], c2[LEN], c3[LEN], c4[LEN]; for (int i = 0; i < n; i++) { a0[i] = cpx(a[i] % S, 0), a1[i] = cpx(a[i] / S, 0); b0[i] = cpx(b[i] % S, 0), b1[i] = cpx(b[i] / S, 0); } mul(a0, b0, c1, n), mul(a0, b1, c2, n), mul(a1, b0, c3, n), mul(a1, b1, c4, n); for (int i = 0; i < n; i++) { int x1 = ctol(c1[i]) % p, x2 = ctol(c2[i]) % p, x3 = ctol(c3[i]) % p, x4 = ctol(c4[i]) % p; ans[i] = (x1 + ll(x2 + x3) * S % p + ll(x4) * S % p * S % p) % p; } } void _inv(const int *const a, int *b, const int n) { if (n == 1) return void(b[0] = get_inv(a[0])); static int tmp[LEN]; _inv(a, b, (n + 1) >> 1); memset(b + ((n + 1) >> 1), 0, sizeof(int[n - ((n + 1) >> 1)])); MTT(a, b, tmp, n), MTT(tmp, b, tmp, n); for (int i = 0; i < n; i++) b[i] = (2LL * b[i] % p - tmp[i] + p) % p; } void inv(const int *const a, int *b, const int n) { static int tmp[LEN]; memcpy(tmp, a, sizeof(int[n])); _inv(tmp, b, n); } } int A[N << 1]; int work() { using namespace Polynomial; int n; read(n); for (int i = 0; i < n; i++) read(A[i]); inv(A, A, n); for (int i = 0; i < n; i++) write(A[i]), putchar(' '); return 0; }}int main(){#ifdef BlueSpirit freopen("4239.in", "r", stdin);#endif return zyt::work();}