多项式全家桶(纯模板)

已经完成的：

多项式乘法
多项式乘法逆
多项式 $\ln$
多项式 $\exp$
多项式开根( $a_0$ 为二次剩余)
多项式快速幂( $a_0$ 可以不为1)
多项式三角函数( $\sin,\cos$ )
多项式除法、取模
分治 $FFT$
任意模数 $NTT$ ( $MTT$ )
多项式多点求值
多项式快速插值
常系数齐次线性递推
BM

还未完成的

多项式非齐次线性递推
普通幂转下降幂/下降幂转普通幂
多项式复合逆

总结

一般多项式题，都是通过一些奇妙的操作（DP、生成函数）后化成一个能用多项式快速完成的东西
一些常见套路：

f_{n} = \sum_{i = 1}^{n} g_{i} h_{n - i}

$f_n=\sum_{i=1}^{n}g_ih_{n-i}$

就是卷积的基础柿子

f_{n} = \sum_{i = 1}^{n - 1} f_{i} g_{n - i}

$f_n=\sum_{i=1}^{n-1}f_ig_{n-i}$

分治 $FFT$ ：直接 $cdq$ 分治，每次计算 $f_{l\sim mid}$ 与 $g_{1\sim r-l}$ 卷积对 $f_{mid+1\sim r}$ 的贡献，复杂度 $\mathcal O(n\log^2 n)$ 。（当然可以多项式求逆，但分治 $FFT$ 可用于更广泛的情况）

特殊情况：自己卷自己：

f_{n} = \sum_{i = 1}^{n - 1} g_{i} h_{n - i}

$f_n=\sum_{i=1}^{n-1}g_ih_{n-i}$

其中 $g,h$ 都是与 $f$ 有关的函数，需要先求出 $f_x$ 才能得到 $g_x,h_x$ 。

还是考虑之前的做法，但在 $solve(l,r)$ 时，可能会出现 $r-l<mid$ 的情况导致需要的 $h$ 没有被算出来。注意到情况只会在 $l=1$ 时出现，因此在 $l=1$ 时，先只计算 $g_{1\sim mid}$ 与 $h_{1\sim mid}$ 卷积对 $f_{mid+1\sim r}$ 的贡献。在 $l\not=1$ 时，在上一条的卷积之外，额外增加 $g_{1\sim r-l}$ 与 $h_{l\sim mid}$ 卷积对 $f_{mid+1\sim r}$ 的贡献，弥补之前缺失的部分即可。

减法卷积：

f_{n} = \sum_{i = 1}^{d} g_{i} h_{n + i}

$f_n=\sum_{i=1}^{d}g_ih_{n+i}$

将 $g(x)$ $reverse$ 一下，于是就变成了：

f_{n} = \sum_{i = 1}^{d} g_{d - i} h_{n + i}

$f_n=\sum_{i=1}^{d}g_{d-i}h_{n+i}$

于是又可以卷积了，设 $g*h=A$ , $f_n=A_{n+d}$

f_{n} = \sum_{i = 1}^{n} (\binom{n}{i}) g_{i} h_{n - i}

$f_n=\sum_{i=1}^{n}\binom nig_{i}h_{n-i}$

设 $F,G,H$ 分别是 $f(x),g(x),h(x)$ 的 $EGF$ ,于是有 $F=GH$

5.指数公式定理：

如果存在两个EGF $F(x)$ 与 $G(x)$ ，满足 $e^{F(x)}=G(x)$ ， $F(x)$ 是 $f_i$ 的 EGF，那么 $G(x)$ 是

g_{n} = \sum_{π = {S_{1}, S_{2}, \dots, S_{n}}} \prod_{i = 1}^{n} f_{| S_{i} |}

$g_n=\sum_{\pi=\{S_1,S_2,\dots,S_n\}}\prod_{i=1}^n f_{|S_i|}$

的EGF，其中 $\pi$ 是 $[n]$ 的划分。

Update on 2022.3.22 Berlekamp–Massey 算法

BM求解最短递推式

```cpp namespace BM{ vector f,mn; int fail[N],id,tot,del[N]; inline vector work(int n,int *a){ tot=0; f.clear();mn.clear(); for(int i=0;ig; int tmp=1ll*del[i]*ksm(del[id],mod-2)%mod; for(int j=0;j

Update on 2022.2.26 常系数齐次线性递推

%x^n\bmod f(x)$

```cpp int main(){ // freopen("poly.in","r",stdin); // freopen("1.out","w",stdout); n=Rint;k=Rint; init_poly((k+1)<<1); for(int i=1;i<=k;++i) f[k-i]=-Rint,f[k-i]=(f[k-i]%mod+mod)%mod; f[k]=1;GINV=f.rev().inv(); poly x,ret; x[1]=ret[0]=1; for(;n;n>>=1,x=(x*x)%f) if(n&1) ret=(ret*x)%f; int ans=0; for(int i=0,x;i

波斯坦茉莉算法

```cpp inline int solve(poly f,poly g,ll n){ for(;n;n>>=1){ poly h=g; for(int i=1;i

Update on 2021.3.29 数组版本新鲜出炉，跑的比较快但难用的一匹。

Update on 2021.3.26 多点求值与快速插值

多点求值与快速插值

复制代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
#define lc (p<<1)
#define rc (p<<1|1)
#define mid ((l+r)>>1)


namespace Multiask{
	int a[N],n,m,ans[N];
	poly T[100010<<2],F[100010<<2],f;
	inline void build(int p,int l,int r){
		if(l==r){T[p][1]=mod-a[l];T[p][0]=1;return ;}
		build(lc,l,mid);build(rc,mid+1,r);
		T[p]=T[lc]*T[rc];
	}
	
	inline poly mulT(const poly &F,const poly &G,int siz,int tp){
		poly f=F,g=G;
		if(f.size()<=100){
			poly ret;
			for(int i=0;i<f.size();++i){
				int fr=0;
				if(tp==1) fr=max(fr,i-(int)(f.size()-g.size()));
				for(int j=fr;j<=i&&j<g.size();++j) inc(ret[i-j],1ll*f[i]*g[j]%mod);
			} 
			return ret;
		}
		int len=lg[f.size()*tp],lim=1<<len+1,gg=g.size();
		init(lim);
		reverse(g.v.begin(),g.v.end());
		NTT(f,lim,1);NTT(g,lim,1);
		for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
		NTT(f,lim,0);
		return vec(f.v.begin()+gg-1,f.v.begin()+gg+siz-1);
	}
	
	inline void ask(int p,int l,int r){
//		cout<<l<<" "<<r<<" ";F[p].print(r-l+1);
		if(l==r){
			ans[l]=F[p][0];
			return ;
		}
		F[lc]=mulT(F[p],T[rc],mid-l+1,1);F[rc]=mulT(F[p],T[lc],r-mid,1);
		ask(lc,l,mid);
		ask(rc,mid+1,r);
	}
	inline void query(int n){
		build(1,1,n);
		T[1]=T[1].inv();T[1].resize(n);f.resize(n);
		F[1]=mulT(f,T[1],n,2);
		ask(1,1,n);
	}
	
}


namespace lagrange{
	int x[N],y[N];
	poly g[100010<<2];
	inline void build(int p,int l,int r){
		if(l==r){g[p][0]=dec(0,x[l]);g[p][1]=1;return ;}
		build(lc,l,mid);build(rc,mid+1,r);
		g[p]=g[lc]*g[rc];
	}
	inline poly query(int p,int l,int r){
		if(l==r) return 1ll*y[l]*ksm(Multiask::ans[l],mod-2)%mod;
		return (query(lc,l,mid)*g[rc])+(query(rc,mid+1,r)*g[lc]);
	}
	inline void work(int n){
		build(1,1,n);	
		Multiask::f=g[1].der();
		for(int i=1;i<=n;++i) Multiask::a[i]=x[i];
		Multiask::query(n);
		query(1,1,n).print(n);
	}
}
#undef lc
#undef rc
#undef mid

Update on 2021.3.25 出现了第二代多项式板子，方便使用且常数较优，但考场上基本上写不出来

Update on 2021.1.13 大幅优化了下方模板的常数

Update on 2021.1.12

任意模数NTT

复制代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=2e7+10;
const int M=6;
int n,v,mod,q;
int k[M][M],b[M][M];

//--------iobuff--------

namespace iobuff{
	const int LEN=1000000;
	char in[LEN+5], out[LEN+5];
	char *pin=in, *pout=out, *ed=in, *eout=out+LEN;
	inline char gc(void)
	{
		#ifdef LOCAL
		return getchar();
		#endif
		return pin==ed&&(ed=(pin=in)+fread(in, 1, LEN, stdin), ed==in)?EOF:*pin++;
	}
	inline void pc(char c)
	{
		pout==eout&&(fwrite(out, 1, LEN, stdout), pout=out);
		(*pout++)=c;
	}
	inline void flush()
	{ fwrite(out, 1, pout-out, stdout), pout=out; }
	template<typename T> inline void read(T &x)
	{
		static int f;
		static char c;
		c=gc(), f=1, x=0;
		while(c<'0'||c>'9') f=(c=='-'?-1:1), c=gc();
		while(c>='0'&&c<='9') x=10*x+c-'0', c=gc();
		x*=f;
	}
	template<typename T> inline void putint(T x, char div)
	{
		static char s[15];
		static int top;
		top=0;
		x<0?pc('-'), x=-x:0;
		while(x) s[top++]=x%10, x/=10;
		!top?pc('0'), 0:0;
		while(top--) pc(s[top]+'0');
		pc(div);
	}
}
using namespace iobuff;

//---------math----------

int base[N],inv[N],jc[N];
inline ll add(ll x,ll y,ll m=mod){return (x+y>=m)?x+y-m:x+y;}
inline ll dec(ll x,ll y,ll m=mod){return (x-y<0)?x-y+m:x-y;}
inline void init(){
	base[0]=base[1]=inv[0]=inv[1]=jc[0]=jc[1]=1;
	for(int i=2;i<N;++i){
		base[i]=1ll*base[i-1]*i%mod;
		inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
		jc[i]=1ll*jc[i-1]*inv[i]%mod;
	}
}
inline int binom(int x,int y){
	return 1ll*base[x]*jc[y]%mod*jc[x-y]%mod;
}
inline int calc(int l,int r,int j,int u){
	return dec(binom(r+j+1,j+u+1),binom(l+j,j+u+1));
}
inline ll ksm(ll x,ll y,ll mod){
	ll ret=1ll;
	for(;y;y>>=1,x=x*x%mod) (y&1)&&(ret=ret*x%mod);
	return ret;
}

//--------poly----------
const ll mod1=998244353;
const ll mod2=1004535809;
const ll mod3=469762049;
const ll mod12=mod1*mod2;
const ll iv1=669690699;
const ll iv2=354521948;
const ll gg=3;
namespace Container{
	ll Wn[2][N];
	int lg[N];
	struct poly{
		vector<ll>v;
		inline ll& operator[](int x){while(x>=v.size())v.push_back(0);return v[x];}
		inline poly(int x=0):v(1){v[0]=x;}
		inline int size(){return v.size();}
		inline void resize(int x){v.resize(x);}
		inline void mem(int l,int lim){
			int t=min(lim,(int)v.size());
			for(int i=l;i<t;++i) v[i]=0;
		}	
	};
	int p1,p2;
	inline void init_poly(int n,ll mod){
		p1=1;
		while((p1<<1)<=n) p1<<=1;p2=p1<<1;
		ll giv=(mod+1)/3;
		ll wn1=ksm(3,(mod-1)/p2,mod),wn2=ksm(giv,(mod-1)/p2,mod);
		Wn[0][p1]=Wn[1][p1]=1;
		for(int i=p1+1;i<=p2;++i) 
			Wn[0][i]=Wn[0][i-1]*wn1%mod,Wn[1][i]=Wn[1][i-1]*wn2%mod;
		for(int i=p1-1;i>=1;--i) Wn[0][i]=Wn[0][i<<1],Wn[1][i]=Wn[1][i<<1];
		for(int i=2;i<n;++i)
			lg[i]=lg[i>>1]+1;
	}
}
using namespace Container;
namespace Poly{	
	int r[N];
	ll fr[N];
	inline void NTT(int lim,poly& f,int tp,ll mod){
		int kk=tp==1?0:1;
		for(int i=0;i<lim;++i) fr[i]=f[r[i]];
		for(int mid=1;mid<lim;mid<<=1){
			int len=mid<<1;
			for(int l=0;l+len-1<lim;l+=len){
				for(int k=l;k<=l+mid-1;++k){
					ll w1=fr[k],w2=Wn[kk][k-l+mid]*fr[k+mid]%mod;
					fr[k]=add(w1,w2,mod);fr[k+mid]=dec(w1,w2,mod);
				}
			}
		}
		if(tp==-1){
			ll t=ksm(lim,mod-2,mod);
			for(int i=0;i<lim;++i) fr[i]=fr[i]*t%mod;
		}
		for(int i=0;i<lim;++i) f[i]=fr[i];
	}
	inline ll merge(ll x1,ll x2,ll x3){
		ll k1=(x2+mod2-x1)%mod2*iv1%mod2*mod1+x1;
		return (((x3-k1%mod3+mod3)%mod3*iv2%mod3*(mod12%mod))%mod+k1)%mod;
	}
	inline poly poly_mul(int n,int m,poly f,poly g){
		int lim=1,len=0;
		while(lim<(n+m)) lim<<=1,len++;
		for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(len-1));
		poly f1=f,f2=f,f3=f,g1=g,g2=g,g3=g; 
		init_poly(lim,mod1);
		NTT(lim,f1,1,mod1);NTT(lim,g1,1,mod1);
		for(int i=0;i<lim;++i) f1[i]=f1[i]*g1[i]%mod1;
		NTT(lim,f1,-1,mod1);
		
		init_poly(lim,mod2);
		NTT(lim,f2,1,mod2);NTT(lim,g2,1,mod2);
		for(int i=0;i<lim;++i) f2[i]=f2[i]*g2[i]%mod2;
		NTT(lim,f2,-1,mod2);
		
		init_poly(lim,mod3);
		NTT(lim,f3,1,mod3);NTT(lim,g3,1,mod3);
		for(int i=0;i<lim;++i) f3[i]=f3[i]*g3[i]%mod3;
		NTT(lim,f3,-1,mod3);
		for(int i=0;i<lim;++i) f[i]=merge(f1[i],f2[i],f3[i]);
		return f;
	}
}
using namespace Poly;

//----------main------------

poly f,g;
int m;
int main(){
	read(n);read(m);read(mod);++n;++m;
	for(int i=0;i<n;++i) read(f[i]);
	for(int i=0;i<m;++i) read(g[i]);
	f=poly_mul(n,m,f,g);
	for(int i=0;i<n+m-1;++i) putint(f[i],' ');flush();
	return 0;
}

Update on ???

多项式求牛顿级数，也算是个比较有用的东西吧：

求 $g(x)$ 的牛顿级数 $c_k$ ：

c_{k} = △^{k} g (0) = \sum_{k} (\binom{n}{k}) (- 1)^{n - k} g (k)

$c_k=\triangle^kg(0)=\sum_k\binom{n}{k}(-1)^{n-k}g(k)$

于是令 $f(x)=(-1)^x$ ，故 $C=F*G$ ，其中 $C,F,G$ 都是 $EGF$

多项式牛顿级数

复制代码
1
2
3
4
5
6
7
8
9
10
11
12
inline int mi(int x){return x&1?mod-1:1;}
inline poly Newton(int n,poly g){
	poly c;
	for(int i=0;i<n;++i)
		g[i]=1ll*g[i]*jc[i]%mod;
	for(int i=0;i<n;++i)
		c[i]=1ll*mi(i)*jc[i]%mod;
	c=poly_mul(n,n,c,g);
	for(int i=0;i<n;++i)
		c[i]=1ll*c[i]*base[i]%mod;
	return c;
}//这里我们假设已经求出了$g(i)$的值，存在$g$中

第二代模板

view code

复制代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
#include<bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
const int N=(1<<21)+20;
const int mod=998244353;
typedef vector<int> vec;
 
inline int add(int x,int y){return (x+y>=mod)?x+y-mod:x+y;}
inline int dec(int x,int y){return (x-y<0)?x-y+mod:x-y;}
inline void inc(int &x,int y){x=add(x,y);}
inline void rec(int &x,int y){x=dec(x,y);}
inline int ksm(int x,int y){
	int ret=1;
	for(;y;y>>=1,x=1ll*x*x%mod) if(y&1) ret=1ll*ret*x%mod;
	return ret; 
}
 
namespace IO {
	inline char nc(){
//		return getchar();
		static char buf[500005],*p1=buf,*p2=buf;
		return p1==p2&&(p2=(p1=buf)+fread(buf,1,500000,stdin),p1==p2)?EOF:*p1++;
	}
	char out[500005],*pout=out,*eout=out+500000;
	inline void flush() { fwrite(out,1,pout-out,stdout),pout=out; }
	inline void pc(char c) { pout==eout&&(fwrite(out,1,500000,stdout),pout=out); (*pout++)=c; }
	template<typename T> inline void put(T x,char suf) {
		static char stk[65];int top=0; while(x) stk[top++]=x%10,x/=10;
		!top?pc('0'),0:0; while(top--) pc(stk[top]+'0'); pc(suf);
	}
	inline int read(){
		char ch=nc(); int sum=0; for(;ch<'0'||ch>'9';ch=nc());
		while(ch>='0'&&ch<='9')sum=(sum<<1)+(sum<<3)+ch-48,ch=nc();
		return sum;
	}
}
#define Rint IO::read()
using IO::put;
using IO::nc;
namespace Prep{
	int Wn[N<<1],lg[N],p2,mx=1,r[N],tot;
	inline void init_poly(int n){
		int p=1;while(p<=n)p<<=1;
		for(int i=2;i<=p;++i) lg[i]=lg[i>>1]+1;
		for(int i=1;i<p;i<<=1){
			int wn=ksm(3,(mod-1)/(i<<1));
			Wn[++tot]=1;
			for(int j=1;j<i;++j) ++tot,Wn[tot]=1ll*Wn[tot-1]*wn%mod;
		}
	}
	inline void init(int lim){
//		cout<<"init"<<lim<<endl;
		int len=lg[lim]-1;
		for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<len);
	}
	int iv[N],tp;
	inline void init_inv(int n){
		if(!tp){tp=2;iv[0]=iv[1]=1;}
		for(;tp<=n;++tp) iv[tp]=1ll*(mod-mod/tp)*iv[mod%tp]%mod;
	}
}
using namespace Prep;
 
namespace Cipolla{
	int I,fl=0;
	mt19937 rnd(time(0));
	struct pt{
		int a,b;
		pt(int _a=0,int _b=0){a=_a;b=_b;}
	};
	inline pt operator *(pt x,pt y){
		pt ret;
		ret.a=add(1ll*x.a*y.a%mod,1ll*x.b*y.b%mod*I%mod);
		ret.b=add(1ll*x.a*y.b%mod,1ll*x.b*y.a%mod);
		return ret;
	}
	inline bool check(int x){
		return ksm(x,(mod-1)/2)==1;
	}
	inline int random(){return rnd()%mod;}
	inline pt qpow(pt a,int b){
		pt ret=pt(1,0);
		for(;b;a=a*a,b>>=1) if(b&1) ret=ret*a;
		return ret;
	}
	inline int cipolla(int n){
		if(!check(n)) return 0;
		int a=random();
		while(!a||check(dec(1ll*a*a%mod,n))) a=random();
		I=dec(1ll*a*a%mod,n);
		int ans=qpow(pt(a,1),(mod+1)/2).a;
		return min(ans,(int)mod-ans);
	}
}
using namespace Cipolla;
 
int ddddd=0;
bool flag=0;
namespace Poly{
	struct poly{
		vec v;
		
		inline poly(int w=0):v(1){v[0]=w;}
		inline poly(const vec&w):v(w){}
		
		inline int operator [](int x)const{return x>=v.size()?0:v[x];}
		inline int& operator [](int x){
			if(x>=v.size()) v.resize(x+1);
			return v[x];
		}
		inline int size(){return v.size();}
		inline void resize(int x){v.resize(x);}
		inline void read(int n){v.resize(n);for(int i=0;i<n;++i) v[i]=Rint;}
		inline void print(int n)const{for(int i=0;i<n-1;++i) put(operator[](i),' ');put(operator[](n-1),'\n');}
		
		inline poly slice(int len)const{
			if(len<=v.size()) return vec(v.begin(),v.begin()+len);
			vec ret(v);ret.resize(len);
			return ret;
		}
		
		inline poly operator *(const int &x)const{
			poly ret(v);
			for(int i=0;i<v.size();++i) ret[i]=1ll*ret[i]*x%mod; 
			return ret;
		}
		inline poly operator -()const{
			poly ret(v);
			for(int i=0;i<v.size();++i) ret[i]=dec(0,ret[i]);
			return ret;
		}
		
		inline poly operator*(const poly &g)const;
		inline poly operator/(const poly &g)const;
		inline poly operator%(const poly &g)const;
		inline poly der()const{
			vec ret(v);
			for(int i=0;i<ret.size()-1;++i) ret[i]=1ll*ret[i+1]*(i+1)%mod;
			ret.pop_back();
			return ret; 
		}
		inline poly jifen()const{
			vec ret(v);
			init_inv(ret.size());ret.push_back(0);
			for(int i=ret.size()-1;i;--i) ret[i]=1ll*ret[i-1]*iv[i]%mod;ret[0]=0;
			return ret;
		}
		inline poly rev()const{
			vec ret(v);
			reverse(ret.begin(),ret.end());
			return ret;
		}
		
		inline poly inv()const;
		inline poly div(const poly &FF)const;
		inline poly ln()const;
		inline poly exp()const;
		inline poly pow(int k)const;
		inline poly sqrt()const;
		inline poly mulT(const poly &g,int siz,int tp)const;
	};
	
	inline poly operator +(const poly &x,const poly &y){
		vec v(max(x.v.size(),y.v.size()));
		for(int i=0;i<v.size();++i) v[i]=add(x[i],y[i]);
		return v;
	}
	inline poly operator -(const poly &x,const poly &y){
		vec v(max(x.v.size(),y.v.size()));
		for(int i=0;i<v.size();++i) v[i]=dec(x[i],y[i]);
		return v;
	}
	
	ull fr[N];
	const ull Mod=998244353;
	inline void NTT(poly& f,int lim,int tp){
//		ddddd+=lim*__builtin_ctz(lim);
		for(int i=0;i<lim;++i) fr[i]=f[r[i]];
		for(int mid=1;mid<lim;mid<<=1){
			for(int len=mid<<1,l=0;l+len-1<lim;l+=len){
				for(int k=l;k<l+mid;++k){
					ull w1=fr[k],w2=fr[k+mid]*Wn[mid+k-l]%Mod;
					fr[k]=w1+w2;fr[k+mid]=w1+Mod-w2; 
				}
			}
		}
		for(int i=0;i<lim;++i) fr[i]>=Mod?fr[i]%=Mod:0;
		if(!tp){
			reverse(fr+1,fr+lim);
			int iv=ksm(lim,mod-2);
			for(int i=0;i<lim;++i) fr[i]=fr[i]*iv%mod;
		}
		for(int i=0;i<lim;++i) f[i]=fr[i];
	}
	
	inline poly poly::operator *(const poly &G)const{
		poly f(v),g=G;
		int rec=f.size()+g.size()-1,d=max(f.size(),g.size());
		int len=lg[rec],lim=1<<len+1;
		init(lim);
		NTT(f,lim,1);NTT(g,lim,1);
		for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
		NTT(f,lim,0);
		return f.slice(rec);
	}
	
	inline poly poly::inv()const{
		poly g,g0,d;
		g[0]=ksm(v[0],mod-2);
		for(int lim=2;(lim>>1)<v.size();lim<<=1){
			g0=g;d=slice(lim);
			init(lim);
			NTT(g0,lim,1);NTT(d,lim,1);
			for(int i=0;i<lim;++i) d[i]=1ll*g0[i]*d[i]%mod;
			NTT(d,lim,0);
			fill(d.v.begin(),d.v.begin()+(lim>>1),0);
			NTT(d,lim,1);
			for(int i=0;i<lim;++i) d[i]=1ll*d[i]*g0[i]%mod;
			NTT(d,lim,0);
			for(int i=lim>>1;i<lim;++i) g[i]=dec(g[i],d[i]);
		}
		return g.slice(v.size());
	}//10E(n)
	
	inline poly poly::div(const poly &FF)const{
		if(v.size()==1) return 1ll*v[0]*ksm(FF[0],mod-2)%mod;	
		int len=lg[v.size()],lim=1<<len+1,nlim=lim>>1;
		poly F=FF,G0=FF.slice(nlim);	
		G0=G0.inv();
		poly H0=slice(nlim),Q0;
		
		init(lim);
		NTT(G0,lim,1);NTT(H0,lim,1);
		for(int i=0;i<lim;++i) Q0[i]=1ll*G0[i]*H0[i]%mod;
		NTT(Q0,lim,0);Q0.resize(nlim);
		
		poly ret=Q0;
		NTT(Q0,lim,1);NTT(F,lim,1);
		for(int i=0;i<lim;++i) Q0[i]=1ll*Q0[i]*F[i]%mod;
		NTT(Q0,lim,0);
		fill(Q0.v.begin(),Q0.v.begin()+nlim,0);
		for(int i=nlim;i<lim&&i<v.size();++i) Q0[i]=dec(Q0[i],v[i]);
		NTT(Q0,lim,1);
		for(int i=0;i<lim;++i) Q0[i]=1ll*Q0[i]*G0[i]%mod;
		NTT(Q0,lim,0);
		for(int i=nlim;i<lim;++i) ret[i]=dec(ret[i],Q0[i]); 
		return ret.slice(v.size());
	}
	
	inline poly poly::ln()const{
		return der().div(*this).jifen(); 
	}//13E
	
	namespace EXP{
		const int logB=4;
		const int B=16;
		poly f,ret,g[30][B];
		inline void exp(int lim,int l,int r){
//			cout<<lim<<" "<<l<<" "<<r<<endl;
			if(r-l<=64){
				for(int i=l;i<r;++i){
					ret[i]=(!i)?1:1ll*ret[i]*iv[i]%mod;
					for(int j=i+1;j<r;++j) inc(ret[j],1ll*ret[i]*f[j-i]%mod);
				}
				return ;
			}
			int k=(r-l)/B;
			poly bl[B];
			for(int i=0;i<B;++i) bl[i].resize(k<<1);
			int len=1<<lim-logB+1;
			for(int i=0;i<B;++i){
				if(i>0){
					init(len);NTT(bl[i],len,0);
					for(int j=0;j<k;++j)
						inc(ret[l+i*k+j],bl[i][j+k]);
				}
				exp(lim-logB,l+i*k,l+(i+1)*k);
				if(i<B-1){
					poly H;H.resize(k<<1);
					for(int j=0;j<k;++j) H[j]=ret[j+l+i*k];
					init(len);NTT(H,len,1);
					for(int j=i+1;j<B;++j)
						for(int t=0;t<(k<<1);++t) inc(bl[j][t],1ll*H[t]*g[lim][j-i-1][t]%mod); 
				}
			}
		}
		inline void init_exp(){
			ret.resize(f.size());
			for(int i=0;i<f.size();++i) f[i]=1ll*f[i]*i%mod,ret[i]=0;
			int mx=lg[f.size()]+1;
			init_inv(1<<mx);
			for(int lim=mx;lim>=logB;lim-=logB){
				int bl=1<<(lim-logB),tot=0,ll=1<<(lim-logB+1);
				init(ll);
				for(int i=0;i<B-1;++i){
					g[lim][i].resize(bl<<1);
					for(int j=0;j<(bl<<1);++j) g[lim][i][j]=f[j+bl*i];
					NTT(g[lim][i],ll,1);
				}
			}
		}
	}
	
	inline poly poly::exp()const{
		EXP::f=*this;
		EXP::init_exp();
		EXP::exp(lg[v.size()]+1,0,1<<lg[v.size()]+1);
		return EXP::ret.slice(v.size());
	}
	
	
	inline poly poly::pow(int k)const{
		return ((*this).ln()*k).exp();
	}
	
	inline poly poly::operator /(const poly &Q)const{
		if(v.size()<Q.v.size()) return 0;
		int p=v.size()-Q.v.size()+1;
		poly fr=rev(),qr=Q.rev();
		fr.resize(p);qr.resize(p);
		return fr.div(qr).rev();
	} 
	inline poly poly::operator %(const poly &Q)const{
		poly F(v);
		return 	(F-(Q*(F/Q))).slice(Q.v.size()-1);
	} 
	inline poly poly::sqrt()const{
		poly g,h,gf,F1,F2,F3,f(v);
		g[0]=cipolla(operator[](0));
		h[0]=ksm(g[0],mod-2);
		gf[0]=g[0];gf[1]=g[0];
		int iv=(mod+1)/2; 
		init(1);
		for(int lim=1;lim<v.size();lim<<=1){
			for(int i=0;i<lim;++i) F1[i]=1ll*gf[i]*gf[i]%mod;
			NTT(F1,lim,0);
			for(int i=0;i<lim;++i) F1[i+lim]=dec(F1[i],f[i]),F1[i]=0;
				
			int nlim=lim<<1;init(nlim);
			for(int i=lim;i<nlim;++i) rec(F1[i],f[i]);
			F2=h;F2.resize(lim);
			NTT(F1,nlim,1);NTT(F2,nlim,1);
			
			for(int i=0;i<nlim;++i) F1[i]=1ll*F1[i]*F2[i]%mod;
			NTT(F1,nlim,0);
			for(int i=lim;i<nlim;++i) g[i]=dec(0,1ll*F1[i]*iv%mod); 
			
			if(nlim<v.size()){
				gf=g;
				NTT(gf,nlim,1);
				for(int i=0;i<nlim;++i) F3[i]=1ll*gf[i]*F2[i]%mod;
				NTT(F3,nlim,0);
				fill(F3.v.begin(),F3.v.begin()+lim,0);
				NTT(F3,nlim,1);
				for(int i=0;i<nlim;++i) F3[i]=1ll*F3[i]*F2[i]%mod;
				NTT(F3,nlim,0);
				for(int i=lim;i<nlim;++i) rec(h[i],F3[i]);
			}
		}
		return g.slice(v.size());
	}
	inline poly poly::mulT(const poly &G,int siz,int tp)const{
		poly f(v),g=G;
		if(f.size()<=100){
			poly ret;
			for(int i=0;i<f.size();++i){
				int fr=0;
				if(tp==1) fr=max(fr,i-(int)(f.size()-g.size()));
				for(int j=fr;j<=i&&j<g.size();++j) inc(ret[i-j],1ll*f[i]*g[j]%mod);
			} 
			return ret;
		}
		int len=lg[f.size()*tp],lim=1<<len+1,gg=g.size();
		init(lim);
		reverse(g.v.begin(),g.v.end());
		NTT(f,lim,1);NTT(g,lim,1);
		for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
		NTT(f,lim,0);
		return vec(f.v.begin()+gg-1,f.v.begin()+gg+siz-1);
	}
}
using namespace Poly;
 
inline int sread(){
	long long m=mod,x=0;char ch=nc();
	while(!isdigit(ch)){ch=nc();}
	while(isdigit(ch)){x=(x<<3)+(x<<1)+(ch^48);x%=m;ch=nc();}
	return x;
}
//remember to init_poly(多项式最高可能次数) 
int n,k;poly f;
/*
int main(){
//	freopen("poly2.in","r",stdin);
//	freopen("1.out","w",stdout);
	n=Rint+1,k=Rint;
	f.read(n);
	init_poly(n<<1);
	poly g=f+poly(2)-poly(f[0]);
	g=g-(f.sqrt().inv().jifen().exp());
	(g.ln()+poly(1)).pow(k).der().print(n-1);
	IO::flush(); 
	return 0;
}
*/
int main(){
	n=Rint;
	f.read(n);
	init_poly(n);
	f.exp().print(n);
	IO::flush();
	return 0;
}

多点求值与快速插值

view code

复制代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
#define lc (p<<1)
#define rc (p<<1|1)
#define mid ((l+r)>>1)


namespace Multiask{
	int a[N],n,m,ans[N];
	poly T[100010<<2],F[100010<<2],f;
	inline void build(int p,int l,int r){
		if(l==r){T[p][1]=mod-a[l];T[p][0]=1;return ;}
		build(lc,l,mid);build(rc,mid+1,r);
		T[p]=T[lc]*T[rc];
	}
	
	inline poly mulT(const poly &F,const poly &G,int siz,int tp){
		poly f=F,g=G;
		if(f.size()<=100){
			poly ret;
			for(int i=0;i<f.size();++i){
				int fr=0;
				if(tp==1) fr=max(fr,i-(int)(f.size()-g.size()));
				for(int j=fr;j<=i&&j<g.size();++j) inc(ret[i-j],1ll*f[i]*g[j]%mod);
			} 
			return ret;
		}
		int len=lg[f.size()*tp],lim=1<<len+1,gg=g.size();
		init(lim);
		reverse(g.v.begin(),g.v.end());
		NTT(f,lim,1);NTT(g,lim,1);
		for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
		NTT(f,lim,0);
		return vec(f.v.begin()+gg-1,f.v.begin()+gg+siz-1);
	}
	
	inline void ask(int p,int l,int r){
//		cout<<l<<" "<<r<<" ";F[p].print(r-l+1);
		if(l==r){
			ans[l]=F[p][0];
			return ;
		}
		F[lc]=mulT(F[p],T[rc],mid-l+1,1);F[rc]=mulT(F[p],T[lc],r-mid,1);
		ask(lc,l,mid);
		ask(rc,mid+1,r);
	}
	inline void query(int n){
		build(1,1,n);
		T[1]=T[1].inv();T[1].resize(n);f.resize(n);
		F[1]=mulT(f,T[1],n,2);
		ask(1,1,n);
	}
	
}


namespace lagrange{
	int x[N],y[N];
	poly g[100010<<2];
	inline void build(int p,int l,int r){
		if(l==r){g[p][0]=dec(0,x[l]);g[p][1]=1;return ;}
		build(lc,l,mid);build(rc,mid+1,r);
		g[p]=g[lc]*g[rc];
	}
	inline poly query(int p,int l,int r){
		if(l==r) return 1ll*y[l]*ksm(Multiask::ans[l],mod-2)%mod;
		return (query(lc,l,mid)*g[rc])+(query(rc,mid+1,r)*g[lc]);
	}
	inline void work(int n){
		build(1,1,n);	
		Multiask::f=g[1].der();
		for(int i=1;i<=n;++i) Multiask::a[i]=x[i];
		Multiask::query(n);
		query(1,1,n).print(n);
	}
}

数组版poly：略快但难以使用

view code

复制代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
#include<bits/stdc++.h>
using namespace std;
const int N=(1<<19)+20;
const int mod=998244353;
typedef vector<int> vec;
typedef unsigned long long ull;

namespace IO{
	inline char nc(){
		static char buf[500005],*p1=buf,*p2=buf;
		return p1==p2&&(p2=(p1=buf)+fread(buf,1,500000,stdin),p1==p2)?EOF:*p1++;
	}
	char out[500005],*pout=out,*eout=out+500000;
	inline void flush() { fwrite(out,1,pout-out,stdout),pout=out; }
	inline void pc(char c) { pout==eout&&(fwrite(out,1,500000,stdout),pout=out); (*pout++)=c; }
	template<typename T> inline void put(T x,char suf) {
		static char stk[40];int top=0; while(x) stk[top++]=x%10,x/=10;
		!top?pc('0'),0:0; while(top--) pc(stk[top]+'0'); pc(suf);
	}
	inline int read(){
		char ch=nc(); int sum=0; for(;ch<'0'||ch>'9';ch=nc());
		while(ch>='0'&&ch<='9')sum=(sum<<1)+(sum<<3)+ch-48,ch=nc();
		return sum;
	}
}
#define Rint IO::read()
using IO::put;
using IO::nc;
//IObuff没判负数 

namespace Math{
	inline int add(int x,int y){return (x+y>=mod)?x+y-mod:x+y;}
	inline int dec(int x,int y){return (x-y<0)?x-y+mod:x-y;}
	inline void inc(int &x,int y){x=add(x,y);}
	inline void rec(int &x,int y){x=dec(x,y);}
	inline int ksm(int x,int y){
		int ret=1;
		for(;y;y>>=1,x=1ll*x*x%mod) if(y&1) ret=1ll*ret*x%mod;
		return ret; 
	}
	int iv[N],tp;
	inline void init_inv(int n){
		if(!tp){iv[0]=iv[1]=1;tp=2;}
		for(;tp<=n;++tp) iv[tp]=1ll*(mod-mod/tp)*iv[mod%tp]%mod;
	}
}
using namespace Math;

namespace Cipolla{
	int I,fl=0;
	mt19937 rnd(time(0));
	struct pt{int a,b;pt(int _a=0,int _b=0){a=_a;b=_b;}};
	inline pt operator *(pt x,pt y){
		pt ret;
		ret.a=add(1ll*x.a*y.a%mod,1ll*x.b*y.b%mod*I%mod);
		ret.b=add(1ll*x.a*y.b%mod,1ll*x.b*y.a%mod);
		return ret;
	}
	inline bool check(int x){return ksm(x,(mod-1)/2)==1;}
	inline int random(){return rnd()%mod;}
	inline pt qpow(pt a,int b){
		pt ret=pt(1,0);
		for(;b;a=a*a,b>>=1) if(b&1) ret=ret*a;
		return ret;
	}
	inline int cipolla(int n){
		if(!check(n)) return 0;
		int a=random();
		while(!a||check(dec(1ll*a*a%mod,n))) a=random();
		I=dec(1ll*a*a%mod,n);
		int ans=qpow(pt(a,1),(mod+1)/2).a;
		return min(ans,(int)mod-ans);
	}
}
using namespace Cipolla;

struct poly{
	vec v;
	inline poly(int w=0):v(1){v[0]=w;}
	inline poly(const vec&w):v(w){}
		
	inline int operator [](int x)const{return x>=v.size()?0:v[x];}
	inline int& operator [](int x){if(x>=v.size()) v.resize(x+1);return v[x];}
	inline int size(){return v.size();}
	inline void resize(int x){v.resize(x);}
	inline void read(int n){v.resize(n);for(int i=0;i<n;++i) v[i]=Rint;}
	inline void print(int n)const{for(int i=0;i<n-1;++i) put(operator[](i),' ');put(operator[](n-1),'\n');}
	
	inline poly slice(int len)const{
		if(len<=v.size()) return vec(v.begin(),v.begin()+len);
		vec ret(v);ret.resize(len);
		return ret;
	}
	inline poly operator *(const int &x)const{
		poly ret(v);
		for(int i=0;i<v.size();++i) ret[i]=1ll*ret[i]*x%mod; 
		return ret;
	}
	inline poly operator +(const poly &x)const{
		vec ret(max(v.size(),x.v.size()));
		for(int i=0;i<v.size();++i) ret[i]=add(v[i],x[i]);
		return ret;
	}
	inline poly operator -(const poly &x)const{
		vec ret(max(v.size(),x.v.size()));
		for(int i=0;i<v.size();++i) ret[i]=dec(v[i],x[i]);
		return ret;
	}
};


int Wn[N<<1],lg[N],r[N],tot;
inline void init_poly(int n){
	int p=1;while(p<=n)p<<=1;
	for(int i=2;i<=p;++i) lg[i]=lg[i>>1]+1;
	for(int i=1;i<p;i<<=1){
		int wn=ksm(3,(mod-1)/(i<<1));
		Wn[++tot]=1;
		for(int j=1;j<i;++j) ++tot,Wn[tot]=1ll*Wn[tot-1]*wn%mod;
	}
}
inline void init_pos(int lim){
	int len=lg[lim]-1;
	for(int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<len);
}

ull fr[N];
const ull Mod=998244353;
inline void NTT(int *f,int lim,int tp){
	for(int i=0;i<lim;++i) fr[i]=f[r[i]];
	for(int mid=1;mid<lim;mid<<=1){
		for(int len=mid<<1,l=0;l+len-1<lim;l+=len){
			for(int k=l;k<l+mid;++k){
				ull w1=fr[k],w2=fr[k+mid]*Wn[mid+k-l]%Mod;
				fr[k]=w1+w2;fr[k+mid]=w1+Mod-w2; 
			}
		}
	}
	for(int i=0;i<lim;++i) f[i]=fr[i]%Mod;
	if(!tp){
		reverse(f+1,f+lim);
		int iv=ksm(lim,mod-2);
		for(int i=0;i<lim;++i) f[i]=1ll*f[i]*iv%mod;
	}
}
inline poly to_poly(int *a,int n){
	poly ret;
	ret.resize(n);
	memcpy(ret.v.data(),a,n<<2);
	return ret;
}
inline poly der(poly F,int n){
	vec ret;ret.resize(n-1);
	for(int i=0;i<n-1;++i) ret[i]=1ll*F[i+1]*(i+1)%mod;
	return ret; 
}
inline poly jifen(poly F,int n){
	vec ret;ret.resize(n);
	init_inv(n);
	for(int i=n-1;i;--i) ret[i]=1ll*F[i-1]*iv[i]%mod;
	return ret;
}
		
namespace Basic{
	int f[N],g[N];
	inline poly mul(poly F,poly G,int n,int m){
		int rec=n+m-1,d=max(n,m);F.resize(n);G.resize(m);
		memcpy(f,F.v.data(),sizeof(int)*(n));
		memcpy(g,G.v.data(),sizeof(int)*(m));
		if(d<=150){
			poly ret;ret.resize(rec);
			for(int i=0;i<n;++i)
				for(int j=0;j<m;++j) inc(ret[i+j],1ll*f[i]*g[j]%mod);
			return ret;
		}
		int len=lg[rec],lim=1<<len+1;
		init_pos(lim);
		NTT(f,lim,1);NTT(g,lim,1);
		for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g[i]%mod;
		NTT(f,lim,0);
		return to_poly(f,rec);
	}
	
	int g0[N],f0[N];
	inline poly Inv(poly F,int n){
		int l=1<<lg[n]+1;
		memset(g,0,sizeof(int)*(l));
		memset(g0,0,sizeof(int)*(l));
		memset(f0,0,sizeof(int)*(l));
		memcpy(f0,F.v.data(),sizeof(int)*(n));
		g[0]=ksm(f0[0],mod-2);
		for(int lim=2,nlim=1;nlim<n;lim<<=1,nlim<<=1){
			memcpy(g0,g,sizeof(int)*(nlim));
			memcpy(f,f0,sizeof(int)*(lim));
			init_pos(lim);
			NTT(g0,lim,1);NTT(f,lim,1);
			for(int i=0;i<lim;++i) f[i]=1ll*g0[i]*f[i]%mod;
			NTT(f,lim,0);
			fill(f,f+nlim,0);
			NTT(f,lim,1);
			for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g0[i]%mod;
			NTT(f,lim,0);
			for(int i=nlim;i<lim;++i) rec(g[i],f[i]);
		}
		return to_poly(g,n);
	}
	inline void Inv(int *F,int *G,int n){
		int l=1<<lg[n]+1;
		memset(g,0,sizeof(int)*(l));
		memset(g0,0,sizeof(int)*(l));
		memset(f,0,sizeof(int)*(l));
		memcpy(f0,F,sizeof(int)*(n));
		memset(f0+n,0,sizeof(int)*(l-n));
		g[0]=ksm(f0[0],mod-2);
		
		for(int lim=2,nlim=1;nlim<n;lim<<=1,nlim<<=1){
			memcpy(g0,g,sizeof(int)*(nlim));
			memcpy(f,f0,sizeof(int)*(lim));
			init_pos(lim);
			NTT(g0,lim,1);NTT(f,lim,1);
			for(int i=0;i<lim;++i) f[i]=1ll*g0[i]*f[i]%mod;
			NTT(f,lim,0);
			fill(f,f+nlim,0);
			NTT(f,lim,1);
			for(int i=0;i<lim;++i) f[i]=1ll*f[i]*g0[i]%mod;
			NTT(f,lim,0);
			for(int i=nlim;i<lim;++i) rec(g[i],f[i]);
		}
		memcpy(G,g,sizeof(int)*(n));
	}
}

namespace Ln{
	int f[N],g0[N],h0[N],q0[N],ret[N],rc[N];
	inline poly div(poly H,poly F,int n,int m){
		if(n==1) return 1ll*H[0]*ksm(F[0],mod-2)%mod;	
		int len=lg[n],lim=1<<len+1,nlim=lim>>1;
		
		memcpy(f,F.v.data(),sizeof(int)*(m));
		memcpy(g0,F.v.data(),sizeof(int)*(nlim));
		memcpy(h0,H.v.data(),sizeof(int)*(nlim));
		memcpy(rc,H.v.data(),sizeof(int)*(min(lim,m)));
		memset(f+n,0,sizeof(int)*(lim-n));
		memset(g0+nlim,0,sizeof(int)*(nlim));
		memset(h0+nlim,0,sizeof(int)*(nlim)); 
		memset(ret+nlim,0,sizeof(int)*(nlim));
		
		Basic::Inv(g0,g0,nlim);
		init_pos(lim);
		NTT(g0,lim,1);NTT(h0,lim,1);
		for(int i=0;i<lim;++i) q0[i]=1ll*g0[i]*h0[i]%mod;
		NTT(q0,lim,0);memset(q0+nlim,0,sizeof(int)*(nlim));
		
		memcpy(ret,q0,sizeof(int)*(nlim));
		NTT(q0,lim,1);NTT(f,lim,1);
		for(int i=0;i<lim;++i) q0[i]=1ll*q0[i]*f[i]%mod;
		NTT(q0,lim,0);
		memset(q0,0,sizeof(int)*(nlim));
		for(int i=nlim;i<lim&&i<m;++i) rec(q0[i],rc[i]);
		NTT(q0,lim,1);
		for(int i=0;i<lim;++i) q0[i]=1ll*q0[i]*g0[i]%mod;
		NTT(q0,lim,0);
		for(int i=nlim;i<lim;++i) rec(ret[i],q0[i]); 
		return to_poly(ret,n);
	}
	
	inline poly getln(poly F,int n){
		return jifen(div(der(F,n),F,n,n),n);
	}//11E
}

namespace Exp{
	const int logB=4;
	const int B=16;
	int f[N],ret[N],H[N];
	poly g[4][B];
	inline void exp(int lim,int l,int r,int dep){
		if(r-l<=128){
			for(int i=l;i<r;++i){
				ret[i]=(!i)?1:1ll*ret[i]*iv[i]%mod;
				for(int j=i+1;j<r;++j)	
					inc(ret[j],1ll*ret[i]*f[j-i]%mod);
			}
			return ;
		}
		int k=(r-l)/B;
		int len=1<<lim-logB+1;
		vector<unsigned long long> bl[B];
		for(int i=0;i<B;++i) bl[i].resize(k<<1);
		for(int i=0;i<B;++i){
			if(i>0){
				init_pos(len);
				for(int j=0;j<(k<<1);++j) H[j]=bl[i][j]%mod; 
				NTT(H,len,0);
				for(int j=0;j<k;++j)
					inc(ret[l+i*k+j],H[j+k]);
			}
			exp(lim-logB,l+i*k,l+(i+1)*k,dep+1);
			if(i<B-1){
				memcpy(H,ret+l+i*k,sizeof(int)*(k));
				memset(H+k,0,sizeof(int)*(k));
				init_pos(len);NTT(H,len,1);
				for(int j=i+1;j<B;++j)
					for(int t=0;t<(k<<1);++t) 
						bl[j][t]+=1ll*H[t]*g[dep][j-i-1][t]; 
			}
		}
	}
	
	inline poly getexp(poly F,int n){
		memcpy(f,F.v.data(),sizeof(int)*(n));
		int mx=lg[n]+1;init_inv(1<<mx);
		for(int i=0;i<n;++i) f[i]=1ll*f[i]*i%mod;
		memset(ret,0,sizeof(int)*(1<<mx));
		for(int lim=mx,dep=0;lim>=8;lim-=logB,dep++){
			int len=1<<(lim-logB+1);
			init_pos(len);
			for(int i=0;i<B-1;++i){
				g[dep][i].resize(len);
				memcpy(g[dep][i].v.data(),f+(len>>1)*i,sizeof(int)*(len));
				NTT(g[dep][i].v.data(),len,1);
			}
		}
		exp(mx,0,1<<mx,0);
		return to_poly(ret,n);
	}
	inline poly getpow(poly F,int n,int k){
		F=Ln::getln(F,n);
		return getexp(F*k,n);
	}
}

namespace Sqrt{
	int f[N],g[N],h[N],gf[N],F1[N],F2[N],F3[N];
	inline void allmem(int n){
		memset(f,0,sizeof(int)*(n));
		memset(g,0,sizeof(int)*(n));
		memset(h,0,sizeof(int)*(n));
		memset(gf,0,sizeof(int)*(n));
		memset(F1,0,sizeof(int)*(n));
		memset(F2,0,sizeof(int)*(n));
		memset(F3,0,sizeof(int)*(n));
	}
	inline poly sqrt(poly F,int n){
//		allmem(1<<lg[n]+1);
		memcpy(f,F.v.data(),sizeof(int)*(n));
		g[0]=cipolla(f[0]);
		h[0]=ksm(g[0],mod-2);
		gf[0]=g[0];gf[1]=g[0];
		int iv=(mod+1)/2; 
		init_pos(1);
		for(int lim=1;lim<n;lim<<=1){
			for(int i=0;i<lim;++i) F1[i]=1ll*gf[i]*gf[i]%mod;
			NTT(F1,lim,0);
			for(int i=0;i<lim;++i) F1[i+lim]=dec(F1[i],f[i]),F1[i]=0;
				
			int nlim=lim<<1;init_pos(nlim);
			for(int i=lim;i<nlim;++i) rec(F1[i],f[i]);
			memcpy(F2,h,sizeof(int)*(lim));
			NTT(F1,nlim,1);NTT(F2,nlim,1);
			for(int i=0;i<nlim;++i) F1[i]=1ll*F1[i]*F2[i]%mod;
			NTT(F1,nlim,0);
			for(int i=lim;i<nlim;++i) g[i]=dec(0,1ll*F1[i]*iv%mod); 
			if(nlim<n){
				memcpy(gf,g,sizeof(int)*(nlim));
				NTT(gf,nlim,1);
				for(int i=0;i<nlim;++i) F3[i]=1ll*gf[i]*F2[i]%mod;
				NTT(F3,nlim,0);
				memset(F3,0,sizeof(int)*(lim));
				NTT(F3,nlim,1);
				for(int i=0;i<nlim;++i) F3[i]=1ll*F3[i]*F2[i]%mod;
				NTT(F3,nlim,0);
				for(int i=lim;i<nlim;++i) rec(h[i],F3[i]);
			}
		}
		return to_poly(g,n);
	}
}
int n,k;
poly f;

inline int sread(){
	long long m=mod,x=0;char ch=nc();
	while(!isdigit(ch)){ch=nc();}
	while(isdigit(ch)){x=(x<<3)+(x<<1)+(ch^48);x%=m;ch=nc();}
	return x;
}

int main(){
//	freopen("poly.in","r",stdin);
//	freopen("1.out","w",stdout);
	n=Rint+1,k=Rint;
	f.read(n);
	init_poly(n);
	poly g=f+poly(2)-poly(f[0]);
	f=Basic::Inv(Sqrt::sqrt(f,n),n);
	g=Ln::getln(g-Exp::getexp(jifen(f,n),n),n);
	der(Exp::getpow((g+poly(1)),n,k),n).print(n-1);
	IO::flush(); 
	return 0;
}

posted @ 2021-04-12 20:20 cjTQX 阅读(661) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

公告

昵称： cjTQX
园龄： 4年4个月
粉丝： 43
关注： 5

+加关注

2025年3月

日

一

二

三

四

五

六

cjTQX

多项式全家桶(纯模板)

已经完成的：

还未完成的

总结

Update on 2022.3.22 Berlekamp–Massey 算法

Update on 2022.2.26 常系数齐次线性递推

Update on 2021.3.29 数组版本新鲜出炉，跑的比较快但难用的一匹。

Update on 2021.3.26 多点求值与快速插值

Update on 2021.3.25 出现了第二代多项式板子，方便使用且常数较优，但考场上基本上写不出来

Update on 2021.1.13 大幅优化了下方模板的常数

Update on 2021.1.12

Update on ???

第二代模板

公告

搜索

常用链接

最新随笔

随笔档案

阅读排行榜