多项式半家桶

多项式全家桶是什么?

1e6多点求值?

5e5快速插值?

很不幸,下面的板子里一个都没有qwq

多项式乘法NTT

FFT是什么?我还真不太会qwq

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int P = 998244353;
const int G = 3;
const int G1 = 332748118;
const int MAXN = 3e6 + 7;
int n, m;
ll A[MAXN], B[MAXN];
int L, len, R[MAXN];

inline ll add(ll x, ll y) {
	x += y;
	if(x >= P)x -= P;
	return x;
}

inline ll ksm(ll x, ll y) {
	ll ans = 1;
	while(y) {
		if(y & 1)ans = ans * x % P;
		x = x * x % P;
		y >>= 1;
	}
	return ans;
}

inline void init() {
	len = 1;
	while(len <= n + m)len <<= 1, L++;
	for(int i = 0; i < len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
	return ;
}

inline void NTT(ll *F, int len, int typ) {
	for(int i = 0; i < len; ++i) {
		if(i < R[i])swap(F[i], F[R[i]]);
	}
	for(int mid = 1; mid < len; mid <<= 1) {
		ll wn = ksm(typ == 1 ? G : G1, (P - 1) / (mid << 1));
		for(int j = 0; j < len; j += (mid << 1)) {
			ll w = 1;
			for(int k = 0; k < mid; ++k, w = w * wn % P) {
				ll x = F[j + k], y = F[j + k + mid] * w % P;
				F[j + k] = add(x, y);
				F[j + k + mid] = add(x, P - y);
			}
		}
	}
	if(typ == -1) {
		ll inv = ksm(len, P - 2);
		for(int i = 0; i < len; ++i) {
			F[i] = F[i] * inv % P;
		}
	}
}

int main() {
	scanf("%d%d", &n, &m);
	for(int i = 0; i <= n; ++i) {
		scanf("%lld", &A[i]);
	}
	for(int i = 0; i <= m; ++i) {
		scanf("%lld", &B[i]);
	}
	init();
	NTT(A, len, 1);
	NTT(B, len, 1);
	for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
	NTT(A, len, -1);
	for(int i = 0; i <= n + m; ++i) {
		printf("%lld ", A[i]);
	}
	return 0;
}

代码分析:

首先是init,注意到我们按顺序合并的序列实际上是二进制翻转后的下标序列,所以对于这个R,我们有其在2l2^l意义下的翻转二进制位,递推方式也很简单,先把之前的22l2到2^{l}平移下来,然后再把最后一位提到2L2^L

其次是NTT过程,其实框架按照单位根反演的式子背诵就好了,因为那个式子连除n个乘上单位根逆元啥的都有....

然后是三层循环,这个很好看,就是我们从长度为1到长度为2n的log层挨个做,然后每层有$2^n/2^g$段,每一段都要和下一段的同相对位置做乘法,而且每一层的某段都要乘上一个wn(段数-1),wn即单位根,相当于线段树的每一层长度之和,就是O(nlogn)O(nlogn)

注意循环卷积要开至少两倍空间

【模板】分治 FFT

生成函数做法:

序列f的OGF为F(x)其中f的第1项是1,g的生成函数为G(x),其中g的第0项是0

F(x)=F(x)G(x)+1F(x)=F(x)G(x)+1

于是可以解得F(x)=11G(x)F(x)=\frac{1}{1-G(x)}多项式求逆即可

分治FFT做法:

CDQ:Yeah!

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int P = 998244353;
const int G1 = 3;
const int G2 = (P + 1) / 3;
const int MAXN = 4e5 + 7;

int n, L, len;
int F[MAXN], G[MAXN], R[MAXN];
ll A[MAXN], B[MAXN];

inline void init(int x) {
	len = 1;
	L = 0;
	while(len <= 2 * x)len <<= 1, ++L;
	for(int i = 1; i <= len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
	return ;
}

inline ll ksm(ll x, ll y) {
	ll ans = 1;
	while(y) {
		if(y & 1)ans = ans * x % P;
		x = x * x % P;
		y >>= 1;
	}
	return ans;
}

inline void add(int &x, ll y) {
	x += y;
	if(x >= P)x -= P;
}

inline void NTT(ll *F, int len, int typ) {
	for(int i = 1; i < len; ++i)if(i < R[i])swap(F[i], F[R[i]]);
	for(int mid = 1; mid < len; mid <<= 1) {
		ll wn = ksm(typ == 1 ? G1 : G2, (P - 1) / (mid << 1));
		for(int j = 0; j < len; j += (mid << 1)) {
			ll w = 1;
			for(int k = 0; k < mid; ++k, w = w * wn % P) {
				int x = F[j + k], y = F[j + k + mid] * w % P;
				F[j + k] = (x + y) % P;
				F[j + k + mid] = (x + P - y) % P;
			}
		}
	}
	if(typ == -1) {
		ll inv = ksm(len, P - 2);
		for(int i = 0; i < len; ++i)F[i] = F[i] * inv % P;
	}
	return ;
}

inline void out1() {
	puts("A->:");
	for(int i = 0; i <= len; ++i) {
		printf("%lld ", A[i]);
	}
	puts("");
	return ;
}
inline void out2() {
	puts("B->:");
	for(int i = 0; i <= len; ++i) {
		printf("%lld ", B[i]);
	}
	puts("");
	return ;
}

inline void cdq(int l, int r) {
	if(l == r)return ;
	int mid = (l + r) >> 1;
	cdq(l, mid);
	for(int i = 0; i < r - l + 1; ++i)A[i] = G[i]; //0~r-l+1,注意我们是多项式乘法,所以这个必须是前...项
	for(int i = l; i <= mid; ++i)B[i - l] = F[i];
	for(int i = mid + 1; i <= r; ++i)B[i - l] = 0;
	init(r - l + 1);//传入项数
	NTT(A, len, 1);
	NTT(B, len, 1);
	for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
	NTT(A, len, -1);
	for(int i = mid + 1; i <= r; ++i)add(F[i], A[i - l]);
	for(int i = 0; i < len; ++i)A[i] = B[i] = 0;
	cdq(mid + 1, r);
	return ;
}

inline void zhqwq() {
	for(int i = 0; i < 10; ++i)A[i] = 1;
	for(int j = 0; j < 10; ++j)B[j] = 1;//1+x+x^2.....,等价于求前缀和
	init(10);
	NTT(A, len, 1);
	NTT(B, len, 1);
	for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
	NTT(A, len, -1);
	return ;
}

int main() {
	scanf("%d", &n);
	for(int i = 1; i < n; ++i) {
		scanf("%d", &G[i]);
	}
	init(n);
	F[0] = 1;
	cdq(0, n - 1);
	for(int i = 0; i < n; ++i) {
		printf("%d ", F[i]);
	}
	return 0;
}

代码分析:调试NTT可以使用zhqwq函数,因为只要生成函数学的好,乘上啥得出啥很容易知道,所以NTT就不会错qwq

分治过程我们要注意是拿出左半部分的F,右半部分都是0,而且g数组要用从0~r-l项,因为我们是用前半部分进行卷积!

然后记得清空要把len范围都清空,因为我们是循环卷积,后面不清空会错的

【模板】多项式乘法逆

证明:

如果我们mod的不是x,就把这个mod项减少一下

假设已知H(x)F(x)1(modxn2)H(x)F(x)\equiv 1 (\mod x^{\lceil\frac{n}{2}\rceil})

显然的F(x)G(x)1(modxn2)F(x)G(x)\equiv 1(\mod x^{\lceil \frac{n}{2}\rceil})

做差有F(x)(G(x)H(x))0(modxn2)F(x)(G(x)-H(x))\equiv 0 (\mod x^{\lceil \frac{n}{2}\rceil})

于是我们除掉这个F,然后同时平方有

(G(x)H(x))2sth(modxn)(G(x)-H(x))^2\equiv sth(\mod x^n)

然后我们发现,无论如何左边都会有一项的系数小于xn2x^{\lceil \frac{n}{2}\rceil}所以就会变没了qwq

所以说这个就做完啦sth=0!

那么G(x)2+H(x)22G(x)H(x)(modxn)G(x)^2+H(x)^2-2G(x)H(x)\equiv (\mod x^n)

同时乘这个FF,因为FG=1F*G=1

G(x)=2H(x)H(x)2F(x)(modxn)G(x)=2*H(x)-H(x)^2*F(x)(\mod x^n)

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN = 4e5 + 7;
const int P = 998244353;
const int G1 = 3;
const int G2 = (P + 1) / 3;
int N, len, L, R[MAXN];
ll G[MAXN], A[MAXN], F[MAXN], B[MAXN];

inline void init(int x) {
	len = 1;
	L = 0;
	while(len <= 2 * x)len <<= 1, ++L;
	for(int i = 1; i < len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
	return ;
}

inline ll ksm(ll x, ll y) {
	ll ans = 1;
	while(y) {
		if(y & 1)ans = ans * x % P;
		x = x * x % P;
		y >>= 1;
	}
	return ans;
}

inline void NTT(ll *F, int len, int typ) {
	for(int i = 1; i < len; ++i)if(i < R[i])swap(F[i], F[R[i]]);
	for(int mid = 1; mid < len; mid <<= 1) {
		ll wn = ksm(typ == 1 ? G1 : G2, (P - 1) / (mid << 1));
		for(int j = 0; j < len; j += (mid << 1)) {
			ll w = 1;
			for(int k = 0; k < mid; ++k, w = w * wn % P) {
				int x = F[j + k], y = F[j + k + mid] * w % P;
				F[j + k] = (x + y) % P;
				F[j + k + mid] = (x + P - y) % P;
			}
		}
	}
	if(typ == -1) {
		ll inv = ksm(len, P - 2);
		for(int i = 0; i < len; ++i)F[i] = F[i] * inv % P;
	}
	return ;
}

inline void out1() {
	printf("A-:");
	for(int i = 0; i < len; ++i)printf("%d ", A[i]);
	puts("\n");
}

inline void getinv(ll *F, int n) {
	if(n == 1) {
		G[0] = ksm(F[0], P - 2);
		return ;
	}//只有一项qwq
	int m = (n + 1) / 2;
	getinv(F, m);
	init(n);
	for(int i = 0; i < n; ++i)A[i] = G[i], B[i] = F[i];
	for(int i = n; i < len; ++i)A[i] = B[i] = 0;
	NTT(A, len, 1);
	NTT(B, len, 1);
	for(int i = 0; i < len; ++i)A[i] = ( 2 * A[i] % P - A[i] * B[i] % P * A[i] % P + P) % P;
	NTT(A, len, -1);
	for(int i = 0; i < n; ++i)G[i] = A[i];//还原G
	return ;
}

inline void zhqwq() {
	for(int i = 0; i < 10; ++i)A[i] = 1;
	for(int j = 0; j < 10; ++j)B[j] = 1;//1+x+x^2.....,等价于求前缀和
	init(10);
	NTT(A, len, 1);
	NTT(B, len, 1);
	for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
	NTT(A, len, -1);
	return ;
}

int main() {
	scanf("%d", &N);
	for(int i = 0; i < N; ++i)scanf("%lld", &F[i]);
	getinv(F, N);//递归这个mod x^n
	for(int i = 0; i < N; ++i) {
		printf("%lld ", G[i]);
	}
	return 0;
}

代码分析:写的时候NTT的w*wn没有取模...服了

然后小心这个getinv我们递归下去是取一半上取整递归哈qwq以及最后到底层是n=1,我们直接求逆即可

【模板】多项式对数函数(多项式 ln)

说到ln,很快啊!,就有chz上课不会初中数学的丢人历史

ln(x)=ln(1x)ln(x)=-ln(\frac{1}{x})

证明:

求导,链式法则(lnF(x))=F(x)F(x)(\ln F(x))'=\frac{F(x)'}{F(x)}

所以我们直接求出F'和F逆,然后乘一下,然后再积分积回去就好了

这意味我们没有这个常数项qwq,而你会发现这个常数项是1?

e0=1e^0=1取ln后的多项式常数项是0....

所以其他有理数不行也很好证明了,因为e是超越数,所以ece^c除了c=0是1外,不可能是其他有理数....超越数不能成为任何整数多项式的根..?

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int MAXN = 4e5 + 7;
const int P = 998244353;
const int G1 = 3;
const int G2 = (P + 1) / 3;
int N, len, L, R[MAXN];
ll F[MAXN], A[MAXN], B[MAXN], C[MAXN], D[MAXN];
ll fac[MAXN], ifac[MAXN], inv[MAXN];

inline ll ksm(ll x, ll y) {
	ll ans = 1;
	while(y) {
		if(y & 1)ans = ans * x % P;
		x = x * x % P;
		y >>= 1;
	}
	return ans;
}

inline void init() {
	fac[0] = 1;
	for(int i = 1; i <= N; ++i)fac[i] = fac[i - 1] * i % P;
	ifac[N] = ksm(fac[N], P - 2);
	for(int i = N - 1; i >= 1; --i)ifac[i] = ifac[i + 1] * (i + 1) % P;
	ifac[0] = 1;
	for(int i = 1; i <= N; ++i)inv[i] = fac[i - 1] * ifac[i] % P;
	return ;
}

inline void qiudao(ll *F, ll *B, int N) {
	for(int i = 1; i < N; ++i) {
		B[i - 1] = F[i] * i % P;
	}
	B[N - 1] = 0;
	return ;
}

inline void jifen(ll *F, ll *B, int N) {
	for(int i = N - 1; i >= 0; --i) {
		B[i + 1] = F[i] * inv[i + 1] % P;
	}
	B[0] = 0;
	return ;
}

inline void init(int x) {
	len = 1;
	L = 0;
	while(len <= 2 * x)len <<= 1, ++L;
	for(int i = 1; i < len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
}

inline void NTT(ll *F, int len, int typ) {
	for(int i = 1; i < len; ++i)if(i < R[i])swap(F[i], F[R[i]]);
	for(int mid = 1; mid < len; mid <<= 1) {
		ll wn = ksm(typ == 1 ? G1 : G2, (P - 1) / (mid << 1));
		for(int j = 0; j < len; j += (mid << 1)) {
			ll w = 1;
			for(int k = 0; k < mid; ++k, w = w * wn % P) {
				int x = F[j + k], y = F[j + k + mid] * w % P;
				F[j + k] = (x + y) % P;
				F[j + k + mid] = (x + P - y) % P;
			}
		}
	}
	if(typ == -1) {
		ll inv = ksm(len, P - 2);
		for(int i = 0; i < len; ++i)F[i] = F[i] * inv % P;
	}
	return ;
}

inline void getinv(ll *F, ll *G, int n) {
	if(n == 1) {
		G[0] = ksm(F[0], P - 2);
		return ;
	}
	getinv(F, G, (n + 1) / 2);
	init(n);
	for(int i = 0; i < n; ++i)D[i] = F[i], C[i] = G[i];
	for(int i = n; i < len; ++i)D[i] = 0, C[i] = 0;
	NTT(D, len, 1);
	NTT(C, len, 1);
	for(int i = 0; i < len; ++i)C[i] = (2 * C[i] % P - C[i] * C[i] % P * D[i] % P + P) % P;
	NTT(C, len, -1);
	for(int i = 0; i < n; ++i)G[i] = C[i];
	return;
}

inline void getln(ll *F, int n) {
	qiudao(F, A, n);
	getinv(F, B, n);
	init(n);
	NTT(A, len, 1);
	NTT(B, len, 1);
	for(int i = 0; i < len; ++i)A[i] = A[i] * B[i] % P;
	NTT(A, len, -1);
	jifen(A, F, n);
	return ;
}

int main() {
	scanf("%d", &N);
	for(int i = 0; i < N; ++i)scanf("%lld", &F[i]);
	init();
	getln(F, N);
	for(int i = 0; i < N; ++i)printf("%lld ", F[i]);
	puts("");
	return 0;
}

代码分析:没什么好说的,有一个很坑我的点,就是变量重名!

比如我们再ln中递归这个inv用A数组,然后inv的辅助数组又用了A....

其实解决方法就是多开两个数组,但是调试的时候可能很痛苦!

注意求导和积分顺序?

【模板】多项式指数函数(多项式 exp)

证明:

先同时取对数吧!

lnB(x)A(x)0(modxn)\ln B(x)-A(x)\equiv 0 (mod x^n)

这个式子直接上牛顿迭代其实就做完了,具体的

牛顿迭代的式子是:

xn+1=xnf(xn)f(xn)x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}

几何意义相当于我们先求出xnx_n号点的切线方程,有

y=f(xn)+f(xn)(xxn)y=f(x_n)+f'(x_n)(x-x_n)

然后xn+1x_{n+1}就是这个式子的根,带入y=0解得上式qwq

这个方法的精度增长的很快,基本上是10的次方增长的,对于迭代多项式,也是每次多两倍的精度

F(B(x))=lnB(x)A(x)F(B(x))=\ln B(x) -A(x)

B1(x)=B(x)F(B(x))F(B(x))B_1(x)=B(x)-\frac{F(B(x))}{F'(B(x))}

注意!我们这个F函数参数是一个多项式B,等价于A是一个常数,没有用!

因为A(x)A(x)的取值不随B(x)B(x)发生变化啊

下面我们有

B1(x)=B(x)F(B(x))F(B(x))=B(x)lnB(x)A(x)1B(x)=B(x)(1lnB(x)+A(x))B_1(x)=B(x)-\frac{F(B(x))}{F'(B(x))}=B(x)-\frac{\ln B(x)-A(x)}{\frac{1}{B(x)}}=B(x)(1-\ln B(x)+A(x))

于是我们对着最后一个式子做即可,时间复杂度就是O(nlogn)O(nlogn)的,当然,和lct一样的"nlogn""nlogn"

注意,题目保证了a0=0a_0=0如果没有这个条件是不能exp的,e^0才是1,其他都是超越数,常数项要是有理数

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int P = 998244353;
const int G1 = 3;
const int G2 = (P + 1) / 3;
const int MAXN = 4e5 + 7;
int N, len, L, R[MAXN];
ll F[MAXN], G[MAXN], A[MAXN], B[MAXN], C[MAXN], D[MAXN], E[MAXN], H[MAXN];

inline ll ksm(ll x, ll y) {
	ll ans = 1;
	while(y) {
		if(y & 1)ans = ans * x % P;
		x = x * x % P;
		y >>= 1;
	}
	return ans;
}

ll ifac[MAXN], fac[MAXN], inv[MAXN];
inline void init2() {
	fac[0] = 1;
	for(int i = 1; i <= N; ++i)fac[i] = fac[i - 1] * i % P;
	ifac[N] = ksm(fac[N], P - 2);
	for(int i = N - 1; i >= 1; --i)ifac[i] = ifac[i + 1] * (i + 1) % P;
	ifac[0] = 1;
	for(int i = 1; i <= N; ++i)inv[i] = ifac[i] * fac[i - 1] % P;
	return ;
}


inline void init(int x) {
	L = 0;
	len = 1;
	while(len < 2 * x)len <<= 1, ++L;
	for(int i = 0; i < len; ++i)R[i] = (R[i >> 1] >> 1) | ((i & 1) << (L - 1));
	return ;
}

inline void NTT(ll *F, int len, int typ) {
	for(int i = 0; i < len; ++i)if(i < R[i])swap(F[i], F[R[i]]);
	for(int mid = 1; mid < len; mid <<= 1) {
		ll wn = ksm(typ == 1 ? G1 : G2, (P - 1) / (mid << 1));
		for(int j = 0; j < len; j += (mid << 1)) {
			ll w = 1;
			for(int k = 0; k < mid; ++k, w = w * wn % P) {
				ll x = F[j + k], y = F[j + k + mid] * w % P;
				F[j + k] = (x + y) % P;
				F[j + k + mid] = (x + P - y) % P;
			}
		}
	}
	if(typ == -1) {
		ll inv = ksm(len, P - 2);
		for(int i = 0; i < len; ++i)F[i] = F[i] * inv % P;
	}
	return ;
}
//(G(x)-H(x))^2=0
//G(x)^2-2GH+H^2=0
//G(x)=2H-H^2F
inline void getinv(ll *F, ll *G, int n) {
	if(n == 1) {
		G[0] = ksm(F[0], P - 2);
		return ;
	}
	getinv(F, G, (n + 1) / 2);
	init(n);
	for(int i = 0; i < n; ++i)A[i] = G[i], B[i] = F[i];
	for(int i = n; i < len; ++i)A[i] = B[i] = 0;
	NTT(A, len, 1);
	NTT(B, len, 1);
	for(int i = 0; i < len; ++i)A[i] = (2 * A[i] % P - A[i] * B[i] % P * A[i] % P + P) % P;
	NTT(A, len, -1);
	for(int i = 0; i < n; ++i)G[i] = A[i];
	return ;
}
inline void jifen(ll *F, ll *B, int n) {
	for(int i = n - 1; i >= 0; --i)B[i + 1] = F[i] * inv[i + 1] % P;
	B[0] = 0;
	return ;
}
inline void qiudao(ll *F, ll *B, int n) {
	for(int i = 1; i < n; ++i) B[i - 1] = F[i] * i % P;
	B[n - 1] = 0;
	return ;
}
//(ln F(x))' = F'(x)/F(x)
inline void getln(ll *F, ll *G, int n) {
	for(int i = 0; i <= n * 2; ++i)D[i] = C[i] = 0;//这句话是关键,如果我们那里init里写的是<=,这个就是错的
	qiudao(F, C, n);
	getinv(F, D, n);
	init(n);
	NTT(C, len, 1);
	NTT(D, len, 1);
	for(int i = 0; i < len; ++i)C[i] = C[i] * D[i] % P;
	NTT(C, len, -1);
	jifen(C, G, n);
	return ;
}
//G_1(x)=G_0(x)*(1-ln (G_0(x))+F(x))) (mod x^n)
inline void getexp(ll *F, ll *G, int n) {
	if(n == 1) {
		G[0] = 1;
		return ;
	}
	getexp(F, G, (n + 1) / 2);
	getln(G, E, n);
	init(n);
	E[0] = (1 - E[0] + F[0] + P) % P;
	for(int i = 1; i < n; ++i)E[i] = (-E[i] + F[i] + P) % P;
	for(int i = 0; i < n; ++i)H[i] = G[i];
	for(int i = n; i < len; ++i)H[i] = E[i] = 0;
	NTT(H, len, 1);
	NTT(E, len, 1);
	for(int i = 0; i < len; ++i)E[i] = E[i] * H[i] % P;
	NTT(E, len, -1);
	for(int i = 0; i < n; ++i)G[i] = E[i];
	for(int i = n; i < len; ++i)E[i] = G[i] = 0;
	return ;
}

int main() {
	scanf("%d", &N);
	for(int i = 0; i < N; ++i)scanf("%lld", &F[i]);
	init2();
	getexp(F, G, N); //坏了!qwq
	for(int i = 0; i < N; ++i) {
		printf("%lld ", G[i]);
	}
	return 0;
}

代码分析:md一开始wqy他式子是个错的所以说我直接错到天上去了....

但是吴清月本人这个代码并没有问题,也就是说他的式子可能是等价的?

问题主要在于对于f'的求导,我们是对于f'而不是f(g(x))f(g(x))求导!

后来才检查出来一个巨大的锅:+1是常数项,不能每一项都+1

同时对于getln函数有个巨大的要求就是每次调用时要清空辅助数组,而且要清空到2倍,同时我们里面init数组的写法就要是那样子的就是<号而不是<=号,否则会getinv的时候循环卷积出锅再就是求导的时候n-1项清空而不是第n项/xk

以后都写<号吧,也不会有错,还能节省一点常数