自闭の多项式

多项式真有意思,我超爱多项式。。。

去nm的多项式,然而这个东西还是得学。

多项式位运算(FWT)

FWT是用来解决多项式位运算卷积的工具,支持O(nlogn)时间内解决形如\sum_{i|j=k}a_i*b_j之类的运算。

FWT的想法

我们无法快速得知卷积,于是考虑构造一个变换,使得A的变换点乘B的变换就等于卷积的变换,然后用卷积的变换倒推出卷积。

思路有点类似FFT,只是这里我们利用的是二进制的运算性质,而并非单位根。

or卷积和and卷积

这两个卷积的构造是比较好理解的。

or为例子,我们FWT(A)得到的多项式每一项A^{'}[i]=\sum_{j|i=i}A[j],所以这个运算就是一个迭代的过程,相当于每次枚举最高位,然后把下标这一位为0的迭代到这一位为1的上面,因为0|1=1

然后考虑一个事实:

x|i=k,y|i=k,(x|y)|k=k

设最后要求的多项式为C,我们将FWT(A)FWT(B)点乘,得到的就是:对于一个下标i,他的任意两个二进制子集x,y都会进行一次相乘。

考虑这个东西与FWT(C)的关系,x,y会在算(x|y)位的时候被贡献到C中,而(x|y)也是i的子集,所以

又会在FWT(C)的时候被贡献到i位。

所以我们得出:

FWT(A)*FWT(B)=FWT(A|B)

and卷积构造是一致的,or基于0|1=1and基于1&0=0,所以迭代的时候把右半部分迭代到左边即可。

xor卷积

这个卷积的构造就很神奇了,好像找不到除了数学归纳法以外的证明。

抛出结论:

FWT(A)=(FWT(A_0+A_1),FWT(A_0-A_1))

这样构造就能数学归纳证明,此时

FWT(A)*FWT(B)=FWT(AxorB)

IFWT

这里相当于是给你了X,Y关于x,y的解析式,让你通过X,Y算出x,y

IFWT_{or}(A)=(FWT_{or}(A_0),FWT_{or}(A_1-A_0))

IFWT_{and}(A)=(IFWT_{and}(A_0-A_1),IFWT_{and}(A_1))

IFWT_{xor}(A)=(IFWT(\frac{A_0+A_1}{2}),IFWT(\frac{A_0-A_1}{2}))

代码实现

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
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#define int long long
using namespace std;
#define gc getchar()
inline int read(){
register int x(0),f(1);register char c(gc);
while(c>'9'||c<'0')f=c=='-'?-1:1,c=gc;
while(c>='0'&&c<='9')x=x*10+c-48,c=gc;
return f*x;
}
const int N=210000,P=998244353,inv2=499122177;
int A[N],B[N],a[N],b[N],n;
void cp(){
memcpy(a,A,sizeof(a));
memcpy(b,B,sizeof(b));
}
void fwtor(int *A,int kind){
register int j,mid,r,k;
for(mid=1;mid<n;mid<<=1)
for(r=(mid<<1),j=0;j<n;j+=r)
for(k=0;k<mid;k++)
(A[j+k+mid]+=A[j+k]*kind%P+P)%=P;
}
void fwtand(int *A,int kind){
register int j,k,mid,r;
for(mid=1;mid<n;mid<<=1)
for(r=(mid<<1),j=0;j<n;j+=r)
for(k=0;k<mid;k++)
(A[j+k]+=A[j+k+mid]*kind%P+P)%=P;
}
void fwtxor(int *A,int kind){
register int j,k,mid,r;
for(mid=1;mid<n;mid<<=1)
for(r=(mid<<1),j=0;j<n;j+=r)
for(k=0;k<mid;k++){
int x=A[j+k],y=A[j+k+mid];
if(kind==1)A[j+k]=(x+y)%P,A[j+k+mid]=(x-y+P)%P;
else A[j+k]=(x+y)*inv2%P,A[j+k+mid]=(x-y+P)*inv2%P;
}
}
signed main(){
register int i;
n=(1<<read());
for(i=0;i<n;i++)A[i]=read();
for(i=0;i<n;i++)B[i]=read();
cp();fwtor(a,1);fwtor(b,1);for(i=0;i<n;i++)(a[i]*=b[i])%=P;fwtor(a,-1);for(i=0;i<n;i++)cout<<a[i]<<' ';cout<<'\n';
cp();fwtand(a,1);fwtand(b,1);for(i=0;i<n;i++)(a[i]*=b[i])%=P;fwtand(a,-1);for(i=0;i<n;i++)cout<<a[i]<<' ';cout<<'\n';
cp();fwtxor(a,1);fwtxor(b,1);for(i=0;i<n;i++)(a[i]*=b[i])%=P;fwtxor(a,-1);for(i=0;i<n;i++)cout<<a[i]<<' ';
return 0;
}

原根

定义

定义一个数P的原根g为一个最小的正整数,满足gP互质的同时,满足g^{\phi(P)}\equiv1(modP),且对于

任意k<=\phi(P)g^{\phi(P)}\not \equiv1(modP)

性质

对于P的原根gg^{[1,\phi(P)]}的取值互不相同,这是用来NTT的重要标准,另外NTT模数为了能够把长度每次平分成原根的次幂,要求NTT模数的\phi必须含有很多的质因子2(至少2^{p_1}>=n),如\phi(998244353)=998244352=2^{23}*7*17

求原根

考虑一个很简单的式子,若x^k\equiv1(modP),那么显然(x^{k})^z\equiv1(modP)。所以我们没有必要对于i=[1,\phi(P)]检验g^i是否为1,只要检验Y=g^{\frac{\phi(P)}{p_i}}(p_i\phi(p)的因子)是否为1即可,因为任意一个<\phi(p)的数必定都是某个Y的因子,若它做幂是1Y显然也是1

这样复杂度就是O(g*\sqrt{P})了。

一般因为模数是质数,\phi(P)=P-1

题目

SDOI2015序列统计

题目链接

solution:题目相当于求从集合里选n个数,使其乘积在模m意义下等于x。如果不是乘积而是和的话是很好做的,我们考虑一个多项式,若数i在集合中则第i项为1,否则为0,对这个多项式自己卷n次,所有x+km项之和就是答案。但是这是乘积,我们考虑变乘法为加法,在模意义下一个很好用的方法就是原根。考虑一个简单的式子:

x^a*x^b=x^{ab}

只要我们找到一个合适的x,将S中所有数对x取对数,然后就可以用加法了,这个数必须满足在模m意义下x^{k}各不相同,显然原根就具有这个性质,并且题目保证m是质数且很小,所以原根存在且易求,于是只要转成次幂之后做多项式快速幂即可,注意因为下标实在模\phi(m)意义下进行的,每次乘完之后要把\phi(m)(即m-1)次项及之后的部分并到前面去。

总复杂度O(nlog^2n)

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
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
#define int long long
#define gc getchar()
inline int read(){
register int x(0),f(1);register char c(gc);
while(c>'9'||c<'0')f=c=='-'?-1:1,c=gc;
while(c>='0'&&c<='9')x=x*10+c-48,c=gc;
return f*x;
}
const int N=210000,P=1004535809,G=3;
int rev[N],n,m,id,S,g,a[N],prim[N],cnt,Gi;
void pre(int x){
int u=x,i;
for(i=2;i*i<=x;i++)
if(x%i==0){
prim[++cnt]=i;
while(x%i==0)x/=i;
}
}
int quickpow(int a,int b,int p){
int ans=1,base=a;
while(b){
if(b&1)(ans*=base)%=p;
(base*=base)%=p;
b>>=1;
}
return ans;
}
int check(int x){
for(int i=1;i<=cnt;i++)if(quickpow(x,(m-1)/prim[i],m)==1)return 0;
return 1;
}
int get(int x){
pre(x);
for(int i=2;i<=x;i++)
if(check(i))return i;
}
int limit=1,l,ans[N],base[N],bt[N];
void NTT(int *A,int kind){
register int i,j,k,mid,r;
for(i=0;i<limit;i++)if(i<bt[i])swap(A[i],A[bt[i]]);
for(mid=1;mid<limit;mid<<=1){
int W=quickpow(kind==1?G:Gi,(P-1)/(mid<<1),P);
for(r=mid<<1,j=0;j<limit;j+=r){
int w=1;
for(k=0;k<mid;k++,(w*=W)%=P){
int x=A[j+k],y=A[j+k+mid]*w%P;
A[j+k]=(x+y)%P;
A[j+k+mid]=(x-y+P)%P;
}
}
}
}
void Rev(int *A){
register int i;int inv=quickpow(limit,P-2,P);
for(i=m-1;i<limit;i++)(A[i%(m-1)]+=A[i])%=P,A[i]=0;
for(i=0;i<m-1;i++)A[i]=A[i]*inv%P;
}
void Quickpow(int *A,int b){
register int i;
while(limit<=2*(m-1))l++,limit<<=1;
for(i=0;i<limit;i++)base[i]=A[i];ans[0]=1;
for(i=0;i<limit;i++)bt[i]=(bt[i>>1]>>1)|((i&1)<<(l-1));
while(b){
if(b&1){NTT(ans,1);NTT(base,1);for(i=0;i<limit;i++)(ans[i]*=base[i])%=P;NTT(ans,-1);NTT(base,-1);Rev(ans);Rev(base);}
NTT(base,1);
for(i=0;i<limit;i++)(base[i]*=base[i])%=P;
NTT(base,-1);
Rev(base);
b>>=1;
}
}
int A[N];
signed main(){
register int i;Gi=quickpow(G,P-2,P);
n=read();m=read();id=read();S=read();
for(i=1;i<=S;i++)a[i]=read();
g=get(m-1);
for(i=0;i<=m-2;i++)rev[quickpow(g,i,m)]=i;
for(i=1;i<=S;i++)if(a[i]%m!=0)A[rev[a[i]%m]]++;
Quickpow(A,n);
cout<<ans[rev[id]];
return 0;
}

还有一点就是当a[i]%m等于0的时候我们不能取对数,但是正好也无法对答案产生贡献,直接判掉即可。

多项式求逆

式子推导

若求出G_0(x)表示模x^{\frac{2}{n}}意义下的结果,那么因为显然F(x)G(x)在模x^\frac{n}{2}意义下也是1,我们就有F(x)(G(x)-G_0(x))\equiv0

于是

G(x)-G_0(x)=0(modx^{\frac{n}{2}})

(G(x)-G_0(x))^2=0(mod x^n)

G(x)^2-2G(x)G_0(x)+G_0(x)^2=0(modx^n)

乘上F(x)

G(x)-2G_0(x)+F(x)G_0(x)^2=0

G(x)=2G_0(x)-F(x)G_0(x)^2

多项式开根

推导与上面类似,依旧是由前半部分推,应该能现场yy就不写了。

G(x)=\frac{G_0^2(x)+f(x)}{2G_0(x)}

多项式对数函数(ln)

因为求导和积分水平过菜(反正式子也不难背),所以证明咕了。

组合意义

F(x)每一项表示若干个集合加起来凑出的大小为i的方案树,那么G(x)=ln(F(x)),且表示集合大小为i的方案数,一般来说F(x)常数项为1G(x)常数项为0

结论

ln(F(x))=∫F(x)F′(x)dx

于是就把积分后的多项式和求导后的多项式卷起来再积回去即可。

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
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
#define gc getchar()
inline int read(){
register int x(0),f(1);register char c(gc);
while(c>'9'||c<'0')f=c=='-'?-1:1,c=gc;
while(c>='0'&&c<='9')x=x*10+c-48,c=gc;
return f*x;
}
const int N=410000,P=998244353,G=3;
int quickpow(int a,int b){
int ans=1,base=a;
while(b){
if(b&1)ans=1ll*ans*base%P;
base=1ll*base*base%P;
b>>=1;
}
return ans;
}
int limit,l,bt[N],Gi;
void getbt(int len){
limit=1,l=0;
while(limit<=len)limit<<=1,l++;
for(int i=0;i<limit;i++)bt[i]=((bt[i>>1]>>1)|((i&1)<<(l-1)));
}
void NTT(int *A,int kind,int limit){
register int i,j,k,mid,r;
for(i=0;i<limit;i++)if(i<bt[i])swap(A[i],A[bt[i]]);
for(mid=1;mid<limit;mid<<=1){
int W=quickpow(kind==1?G:Gi,(P-1)/(mid<<1));
for(r=mid<<1,j=0;j<limit;j+=r){
int w=1;
for(k=0;k<mid;k++,w=1ll*w*W%P){
int x=A[j+k],y=1ll*w*A[j+k+mid]%P;
A[j+k]=(x+y)%P;
A[j+k+mid]=(x-y+P)%P;
}
}
}
}
void mul(int *A,int *B,int len){
register int i;
getbt(len);
NTT(A,1,limit);NTT(B,1,limit);
for(i=0;i<limit;i++)A[i]=1ll*A[i]*B[i]%P;
NTT(A,-1,limit);
int inv=quickpow(limit,P-2);
for(i=0;i<limit;i++)A[i]=1ll*A[i]*inv%P;
}
void getinv(int *A,int *B,int len){
static int tmp[N];
register int i;
if(len==1)B[0]=quickpow(A[0],P-2);
else{
getinv(A,B,len>>1);
for(i=0;i<len;i++)tmp[i]=A[i];
getbt(len);
NTT(tmp,1,limit);NTT(B,1,limit);
for(i=0;i<limit;i++)B[i]=1ll*B[i]*(2ll-1ll*B[i]*tmp[i]%P+P)%P;
NTT(B,-1,limit);
int inv=quickpow(limit,P-2);
for(i=0;i<len;i++)B[i]=1ll*B[i]*inv%P;
for(i=len;i<limit;i++)B[i]=0;
for(i=0;i<limit;i++)tmp[i]=0;
}
}
void Dao(int *A,int *B,int len){for(int i=0;i<len-1;i++)B[i]=1ll*A[i+1]*(i+1)%P;B[len-1]=0;}
void Fen(int *A,int *B,int len){for(int i=1;i<len;i++)B[i]=1ll*A[i-1]*quickpow(i,P-2)%P;B[0]=0;}
void getln(int *A,int *B,int len){
register int i;
static int tmp1[N],tmp2[N];
Dao(A,tmp1,len);
getinv(A,tmp2,len);
mul(tmp1,tmp2,len);
Fen(tmp1,B,len);
for(i=0;i<limit;i++)tmp1[i]=tmp2[i]=0;
}
int n,f[N],g[N];
int main(){
register int i;Gi=quickpow(G,P-2);
n=read();
for(i=0;i<n;i++)f[i]=read();
getbt(n);
getln(f,g,limit);
for(i=0;i<n;i++)cout<<g[i]<<' ';
return 0;
}

多项式exp

组合意义

多项式ln的逆运算,即知道集合大小为i的方案数,求取任意个组合成大小为j的方案数。

结论

G(x)=exp(F(x))=G_0(x)(1+F(x)-ln(G_0))

其中G_0表示前半部分做完后的结果,边界G(0)=1

注意此处是modx^n意义下的,所以尽管G_0只有\frac{n}{2}项,ln也要在modx^n意义下做。

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
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
#define gc getchar()
inline int read(){
register int x(0),f(1);register char c(gc);
while(c>'9'||c<'0')f=c=='-'?-1:1,c=gc;
while(c>='0'&&c<='9')x=x*10+c-48,c=gc;
return f*x;
}
const int N=410000,P=998244353,G=3;
int limit,l,bt[N],Gi;
void getbt(int n){
limit=1,l=0;
while(limit<=n)limit<<=1,l++;
for(int i=0;i<limit;i++)bt[i]=(bt[i>>1]>>1)|((i&1)<<(l-1));
}
int quickpow(int a,int b){
int ans=1,base=a;
while(b){
if(b&1)ans=1ll*ans*base%P;
base=1ll*base*base%P;
b>>=1;
}
return ans;
}
void NTT(int *A,int kind,int limit){
register int i,j,k,mid,r;
for(i=0;i<limit;i++)if(i<bt[i])swap(A[i],A[bt[i]]);
for(mid=1;mid<limit;mid<<=1){
int W=quickpow(kind==1?G:Gi,(P-1)/(mid<<1));
for(r=mid<<1,j=0;j<limit;j+=r){
int w=1;
for(k=0;k<mid;k++,w=1ll*w*W%P){
int x=A[j+k],y=1ll*A[j+k+mid]*w%P;
A[j+k]=(x+y)%P;
A[j+k+mid]=(x-y+P)%P;
}
}
}
}
void mul(int *A,int *B,int len){
getbt(len);
register int i;
NTT(A,1,limit);NTT(B,1,limit);
for(i=0;i<limit;i++)A[i]=1ll*A[i]*B[i]%P;
NTT(A,-1,limit);
int inv=quickpow(limit,P-2);
for(i=0;i<limit;i++)A[i]=1ll*A[i]*inv%P;
}
void Dao(int *A,int *B,int len){for(int i=0;i<len;i++)B[i]=1ll*A[i+1]*(i+1)%P;}
void Fen(int *A,int *B,int len){for(int i=0;i<len;i++)B[i]=1ll*A[i-1]*quickpow(i,P-2)%P;}
void getinv(int *A,int *B,int len){
static int tmp[N];
register int i;
if(len==1)B[0]=quickpow(A[0],P-2);
else{
getinv(A,B,len>>1);
getbt(len);
for(i=0;i<len;i++)tmp[i]=A[i];
NTT(tmp,1,limit);NTT(B,1,limit);
for(i=0;i<limit;i++)B[i]=1ll*B[i]*(2ll-1ll*tmp[i]*B[i]%P+P)%P;
NTT(B,-1,limit);
int inv=quickpow(limit,P-2);
for(i=0;i<len;i++)B[i]=1ll*B[i]*inv%P;
for(i=len;i<limit;i++)B[i]=0;
for(i=0;i<limit;i++)tmp[i]=0;
}
}
void getln(int *A,int *B,int len){
static int tmp1[N],tmp2[N];
Dao(A,tmp1,len);getinv(A,tmp2,len);
mul(tmp1,tmp2,len);
Fen(tmp1,B,len);
for(int i=0;i<limit;i++)tmp1[i]=tmp2[i]=0;
}
void getexp(int *A,int *B,int len){
register int i;
static int tmp[N];
if(len==1)B[0]=1;
else{
getexp(A,B,len>>1);
getln(B,tmp,len);
for(i=0;i<len;i++)tmp[i]=(A[i]-tmp[i]+P)%P;
tmp[0]++;
mul(B,tmp,len);
for(i=len;i<limit;i++)B[i]=0;
for(i=0;i<limit;i++)tmp[i]=0;
}
}
int n,f[N],g[N];
int main(){
register int i;
n=read();Gi=quickpow(G,P-2);
for(i=0;i<n;i++)f[i]=read();
getbt(n);
getexp(f,g,limit);
for(i=0;i<n;i++)cout<<g[i]<<' ';
return 0;
}

我的习惯在最开始把多项式化成2的次幂,这样可以避免一些特判,另外要注意两个清零的地方,一个是辅助数组,另一个是多项式乘法之后多出来的后一半,写的时候一定要想一想当前多项式的项数是多少以及每一个多项式运算在模多少的意义下进行

分治FFT/NTT

用途

处理形如这样的卷积:

f(i)=\sum_{j+k=i}f(j)*g(k)

因为f每一项是和前面的系数相关的,可以用一个类似CDQ分治的想法,每次只考虑左边对右边的贡献,然后累加。

考虑一段区间[l,mid][mid+1,r]贡献,他们的下标相差在[1,r-l]之间,所以只要把gr-l+1项拎出来与f前半部分做卷积,就可以得到对后半部分的贡献。

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
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
#define gc getchar()
inline int read(){
register int x(0),f(1);register char c(gc);
while(c>'9'||c<'0')f=c=='-'?-1:1,c=gc;
while(c>='0'&&c<='9')x=x*10+c-48,c=gc;
return f*x;
}
const int N=410000,P=998244353,G=3;
int limit,l,bt[N],Gi;
void getbt(int len){
limit=1,l=0;
while(limit<=len)limit<<=1,l++;
for(int i=0;i<limit;i++)bt[i]=(bt[i>>1]>>1)|((i&1)<<(l-1));
}
int quickpow(int a,int b){
int ans=1,base=a;
while(b){
if(b&1)ans=1ll*ans*base%P;
base=1ll*base*base%P;
b>>=1;
}
return ans;
}
void NTT(int *A,int kind,int limit){
register int i,j,k,mid,r;
for(i=0;i<limit;i++)if(i<bt[i])swap(A[i],A[bt[i]]);
for(mid=1;mid<limit;mid<<=1){
int W=quickpow(kind==1?G:Gi,(P-1)/(mid<<1));
for(r=mid<<1,j=0;j<limit;j+=r){
int w=1;
for(k=0;k<mid;k++,w=1ll*w*W%P){
int x=A[j+k],y=1ll*w*A[j+k+mid]%P;
A[j+k]=(x+y)%P;
A[j+k+mid]=(x-y+P)%P;
}
}
}
}
void mul(int *A,int *B,int *C,int len){
static int tmp1[N],tmp2[N];
getbt(len);
register int i;
for(i=0;i<len;i++)tmp2[i]=B[i];
for(i=0;i<(len>>1);i++)tmp1[i]=A[i];
NTT(tmp1,1,limit);NTT(tmp2,1,limit);
for(i=0;i<limit;i++)tmp1[i]=1ll*tmp1[i]*tmp2[i]%P;
NTT(tmp1,-1,limit);
int inv=quickpow(limit,P-2);
for(i=0;i<len;i++)C[i]=1ll*tmp1[i]*inv%P;;
for(i=0;i<limit;i++)tmp1[i]=tmp2[i]=0;
}
void solve(int *A,int *B,int l,int r){
static int tmp[N];
if(l==r)return;
else{
int mid=(l+r)>>1;
solve(A,B,l,mid);
mul(A+l,B,tmp,r-l+1);
for(int i=mid+1;i<=r;i++)A[i]=(A[i]+tmp[i-l])%P;
for(int i=0;i<=r-l+1;i++)tmp[i]=0;
solve(A,B,mid+1,r);
}
}
int n,f[N],g[N];
int main(){
register int i;
n=read();Gi=quickpow(G,P-2);
for(i=1;i<n;i++)g[i]=read();
getbt(n);f[0]=1;
solve(f,g,0,limit-1);
for(i=0;i<n;i++)cout<<f[i]<<' ';
return 0;
}
0%