任意模数NTT

三模数NTT

不会

本来想写三模数NTT，发现常数太TM大了，还不知道怎么优化，于是就写拆系数FFT吧

首先一个简单的方法就是直接FFT开long double 去跑，发现精度掉得十分严重（值域大于long long），于是考虑如何解决

通过学习bitset的经验我们知道要压位

考虑拆系数

假设多项式$A\ast B,模数是P$

$确定一个常数C一般为2^{15}=32,768或\sqrt{P}$

$考虑把A_i拆成a_i\ast C+b_i$
$把B_i拆成c_i\ast C+d_i$
即
$A_i=a_i\ast C+b_i$
$B_i=c_i\ast C+d_i$
考虑两个数$X,Y$相乘

$设X=a\ast C+b$
$Y=c\ast C+d$
$X\ast Y=(a\ast C+b)\ast (c\ast C+d)$
$=ac\ast C^2+(ad+bc)C+bd$
发现C先不用管，要求的是$ac,ad+bc,bd$
然后把$a,b,c,d$分别求值(DFT)，乘起来后再对$ac,ad+bc,bd$做插值(IDFT)就行了
一共是4次DFT+3次IDFT=7次FFT
常数巨大

发现一开始的虚数部并没有用到，参考FFT3次变2次的优化，试试这个能不能行

$设X=a+bi$
$Y=c+di$

$X\ast Y=(a+bi)\ast (c+di)=ac-bd+(bc+ad)i$

然鹅和我们要求的$ac,bd,ad+bc$
现在已经就出来了$ad+bc$

考虑如何求剩下两个

$设X^{\prime }=a-bi$
$Y=c+di$
$X^{\prime }\ast Y=(a-bi)\ast (c+di)=ac+bd+(ad-bc)i$

$然后就知道了实数部$ac+bd$的值$
$结合上面ac-bd的值就可以求出来ac和bd了$
$这样貌似只用对X,Y,X^{\prime }$做3次求值（DFT）
$再对X\ast Y,X^{\prime }\ast Y$做2次插值就可以了

也就是说一共5次

据说myy有4次的%%%%% orz

code:

#include<bits/stdc++.h>
#define N 4000005
#define ll long long
#define double long double
using namespace std;
const double pi = acos(-1);
const double eps = 0.49;
const int C = 32768;
struct cp {
	double a, b;
}X[N], Y[N], Xd[N];
cp operator +(cp x, cp y) {return cp{x.a + y.a, x.b + y.b};}
cp operator -(cp x, cp y) {return cp{x.a - y.a, x.b - y.b};}
cp operator *(cp x, cp y) {return cp{x.a * y.a - x.b * y.b, x.b * y.a + x.a * y.b};}
int rev[N], n, m, P;
void fft(cp *a, int len, int o) {
	for(int i = 1; i <= len; i ++) rev[i] = (rev[i >> 1] >> 1) | ((i&1) * len >> 1);
	for(int i = 1; i <= len; i ++) if(i < rev[i]) swap(a[i], a[rev[i]]);

	for(int i = 2; i <= len; i <<= 1) {
		cp wn = cp{cos(2 * pi / i), o * sin(2 * pi / i)};
		for(int j = 0, p = i / 2; j + i - 1 <= len; j += i) {
			cp w0 = cp{1, 0};
			for(int k = j; k < j + p; k ++, w0 = w0 * wn) {
				cp x = a[k], y =  w0 * a[k + p];
				a[k] = x + y;
				a[k + p] = x - y;
			}
		}
	}
	if(o == -1) for(int i = 0; i <= len; i ++) a[i].a /= len, a[i].b /= len;
}
int main() {
	scanf("%d%d%d", &n, &m, &P);
	for(int i = 0; i <= n; i ++) {
		int x;
		scanf("%d", &x);
		X[i].a = x / C;
		X[i].b = x % C;
		Xd[i].a = x / C;
		Xd[i].b = -(x % C);
	}
	for(int i = 0; i <= m; i ++) {
		int x;
		scanf("%d", &x);
		Y[i].a = x / C;
		Y[i].b = x % C;
	}
	int len = 1;
	for(; len <= n + m;) len <<= 1;
	fft(X, len, 1), fft(Xd, len, 1), fft(Y, len, 1);
	for(int i = 0; i <= len; i ++) X[i] = X[i] * Y[i], Xd[i] = Xd[i] * Y[i];
	fft(X, len, -1), fft(Xd, len, -1);
	for(int i = 0; i <= n + m; i ++) {
		ll ac = (X[i].a + Xd[i].a) / 2 + eps;
		ll bd = Xd[i].a - ac + eps;
		ll bcad = X[i].b + eps;
		printf("%lld ", (ac * C % P * C % P + bcad * C % P + bd) % P);
	}
	return 0;
}

板子一遍过，感动.jpg

任意模数多项式求逆

code:

#include<bits/stdc++.h>
#define N 4000005
#define ll long long
#define double long double
using namespace std;
const double pi = acos(-1);
const double eps = 0.49;
const int C = 32768;
struct cp {
	double a, b;
}X[N], Y[N], Xd[N];
cp operator +(cp x, cp y) {return cp{x.a + y.a, x.b + y.b};}
cp operator -(cp x, cp y) {return cp{x.a - y.a, x.b - y.b};}
cp operator *(cp x, cp y) {return cp{x.a * y.a - x.b * y.b, x.b * y.a + x.a * y.b};}
int rev[N], P;
void fft(cp *a, int len, int o) {
	for(int i = 1; i <= len; i ++) rev[i] = (rev[i >> 1] >> 1) | ((i&1) * len >> 1);
	for(int i = 1; i <= len; i ++) if(i < rev[i]) swap(a[i], a[rev[i]]);

	for(int i = 2; i <= len; i <<= 1) {
		cp wn = cp{cos(2 * pi / i), o * sin(2 * pi / i)};
		for(int j = 0, p = i / 2; j + i - 1 <= len; j += i) {
			cp w0 = cp{1, 0};
			for(int k = j; k < j + p; k ++, w0 = w0 * wn) {
				cp x = a[k], y =  w0 * a[k + p];
				a[k] = x + y;
				a[k + p] = x - y;
			}
		}
	}
	if(o == -1) for(int i = 0; i <= len; i ++) a[i].a /= len, a[i].b /= len;
}
void mul(ll *a, ll *b, int n, int m) {
	for(int i = 0; i <= n; i ++) {
		X[i].a = a[i] / C;
		X[i].b = a[i] % C;
		Xd[i].a = a[i] / C;
		Xd[i].b = -(a[i] % C);
	}
	for(int i = 0; i <= m; i ++) {
		Y[i].a = b[i] / C;
		Y[i].b = b[i] % C;
	}
	int len = 1;
	for(; len <= n + m;) len <<= 1;
	fft(X, len, 1), fft(Xd, len, 1), fft(Y, len, 1);
	for(int i = 0; i <= len; i ++) X[i] = X[i] * Y[i], Xd[i] = Xd[i] * Y[i];
	fft(X, len, -1), fft(Xd, len, -1);
	for(int i = 0; i <= n + m; i ++) {
		ll ac = (X[i].a + Xd[i].a) / 2 + eps;
		ll bd = Xd[i].a - ac + eps;
		ll bcad = X[i].b + eps;
		a[i] = (ac % P * C % P * C % P + bcad % P * C % P + bd % P) % P;
	}
	for(int i = 0; i <= len; i ++) X[i].a = X[i].b = Y[i].a = Y[i].b = Xd[i].a = Xd[i].b = 0;
}
int qpow(ll x, int y) {
	ll ret = 1;
	for(; y; y >>= 1, x = x * x % P) if(y & 1) ret = ret * x % P;
	return ret;
}
ll c[N], d[N], e[N], n;
void inv(ll *a, ll *b, int sz){
	if(sz == 0) {b[0] = qpow(a[0], P - 2); return;}
	inv(a, b, sz / 2);
	int len = 1;
	for(; len <= sz + sz; len <<= 1);
	for(int i = 0; i <= sz; i ++) c[i] = a[i];
	for(int i = sz + 1; i <= len; i ++) c[i] = 0;
	
	//for(int i = 0; i <= len; i ++) b[i] = (b[i] * 2 % mod - b[i] * b[i] % mod * c[i] % mod + mod) % mod;//Ö±½ÓËã
	
	for(int i = 0; i <= sz; i ++) e[i] = b[i] * 2 % P;
	for(int i = sz + 1; i <= len; i ++) e[i] = 0;
	
	for(int i = 0; i <= sz; i ++) d[i] = b[i]; 
	for(int i = sz + 1; i <= len; i ++) d[i] = 0;
	
	mul(b, d, sz, sz); mul(b, c, sz, sz);
	for(int i = 0; i <= sz; i ++) b[i] = (e[i] - b[i] + P) % P;
	
	for(int i = sz + 1; i <= len; i ++) b[i] = 0;
	
}
ll a[N], b[N];
int main() {
	scanf("%d", &n);
	n --;
	P = 1000000007;
	for(int i = 0; i <= n; i ++) scanf("%lld", &a[i]);
	inv(a, b, n);
	for(int i = 0; i <= n; i ++) printf("%lld ", b[i]);
	return 0;
}