洛谷5577

洛谷5577

Description

$n$ 个不超过 $m$ 位的 $k$ 进制数,每次随意取出一些数 (注意可以一个数都不选), 然后把这些数进行 $k$ 进制下不进位加法。

求最终 $0\sim k^m-1$ 的每个数字被表示出来的方案数对 $998244353$ 取模

Input

第一行三个整数 $n,k,m$ 意义为题目中所述。

第二行 $n$ 个整数,之间用空格隔开,表示纸条上的数字。

注意 , 纸条上的数均为 $k$ 进制数)

Output

共 $k^m$ 行,你需要在第 $i$ 行输出最终得到数字 $i−1$ 的方案数模 $998244353$ 的值。

Sample Input

1
2
3 5 1
1 2 3

Sample Output

1
2
3
4
5
2
2
1
2
1

Constraints

$n\leq 10^6$
$k=5$ 时 $m\leq 7$
$k=6$ 时 $m\leq 6$

Solution

答案是形如 $\prod_{i=1}^n(1+x^{a_i})$ 的所有项的系数并且此处卷积的定义是:两个下标乘起来得到的下标是这两个下标的 $k$ 进制下不进位加法得到的值。

于是就有了个 $O(nk^m)$ 的优质做法:每次加入一个数就枚举所有下标并更新它乘上 $1+x^{a_i}$ 到达的位置。可以获得 20pts 的好成绩

如果这是普通的多项式乘法那么直接对于每个数 DFT 一遍算出许多点值再对应乘起来再把这些点值 IDFT 回去就是答案。

然后你大概感性理解意会一下题目中的「 $k$ 进制下不进位加法」的卷积其实也可以类似的做:把每个下标看成一个 $m$ 维的向量 $\mathbf{V}=\{v_0,v_1,v_2,\cdots v_{m-1}\}$ 且其中每一位 $v_i<k$ (就是 $k$ 进制下的每一位),那么一个下标 $\mathbf{V}=\{v_0,v_1,v_2\cdots,v_m\}$ 和 $\mathbf{U}=\{u_0,u_1,u_2,\cdots u_{m-1}\}$ 的卷积 $\mathbf{V\cdot U}=\{(v_0+u_0)\mod k,(v_1+u_1)\mod k,\cdots,(v_{m-1}+u_{m-1})\mod k\}$ 就是向量中对应的每一位在模 $x^k$ 意义下的循环卷积。

循环卷积的例子:模 $x^5$ 意义下的 $x^3\cdot x^4=x^7$ 就变成了 $x^2$

我们尝试对于一个 $1+x^{a_i}$ 如何进行 $m$ 维下的 DFT

回想一下一般的一维 DFT 含义:对于第 $i$ 项系数它对 DFT 的结果中位置 $j$ 的贡献的系数是 $\omega^{ij}$

二维的话位置 $(Ax,Ay)$ 对位置 $(Bx,By)$ 的贡献就是 $\omega^{AxBx}\cdot\omega^{AyBy}=\omega^{AxBx+AyBy}$ 直接先对每一行 DFT 再对每一列 DFT 就行

那么似乎可以直接套到多维上去:对于一个向量 $\mathbf{U}=\{u_0,u_1,\cdots u_{m-1}\}$ ,其对结果中位置 $\mathbf{V}=\{v_0,v_1,\cdots,v_{m-1}\}$ 的贡献是 $\omega^{u_0v_0}\cdot\omega^{u_1v_1}\cdots\omega^{u_{m-1}v_{m-1}}=\omega^{u_0v_0+u_1v_1+\cdots+u_{m-1}v_{m-1}}$

回到题目。

因为模 $998244353$ 意义下不存在 $\omega_5$ 和 $\omega_6$ 的整数解,但我们又需要表示出 $\omega_5$ 和 $\omega_6$ ,以下内容均在扩域后的意义下进行

即把每个「数」看成是 $a_0\omega_k^0+a_1\omega_k^1+\cdots+a_{k-1}\omega_k^{k-1}$

易证这个域里的「数」满足加法和乘法的封闭性

结合 20pts 做法,注意到每次乘上的只有两项,其中 $1=x^0=x^{\{0,0,\cdots,0\}}$ 它对所有位置都会造成 $\omega_k^{0+0+\cdots+0}=1$ 的贡献;而 $x^{V}=x^{\{v_0,v_1,\cdots,v_{m-1}\}}$ 对于每个位置 $\mathbf{U}=\{u_0,u_1,\cdots u_{m-1}\}$ 会造成 $\omega_k^{u_0v_0}\cdot\omega_k^{u_1v_1}\cdots\omega_k^{u_{m-1}v_{m-1}}=\omega_k^{u_0v_0+u_1v_1+\cdots+u_{m-1}v_{m-1}}$ 的贡献

那我们对于 $n$ 个数可以每次 DFT 得到 $k^m$ 个形如 $1+\omega_k^{t}$ 的的点值,一共求 $n$ 遍对于每一位把它对应的 $n$ 个 $1+\omega_k^t$ 乘起来最后再 IDFT 回去就是答案

这样做复杂度可能是 $O(n(k^m\times k^2))$ 的无疑是过不了的。

但是发现 $n$ 的范围很大但 $k^m$ 却小得多似乎指引我们按照值域做。

仔细一想对于某一位来说,似乎并不需要求出所有的 $n$ 个 $1+\omega_k^t$ 分别是多少,只需要对于每个 $t\in[0,k-1]$ 知道有多少个 $1+\omega_k^t$ 就行了(只需要个数,至于它们具体属于 $n$ 个中的哪几个是无关紧要的)

那问题就转化为有 $n$ 个数 $\mathbf{V_1,V_2,\cdots,V_{n}}$ ,每一项对于位置 $\mathbf{U}$ 的贡献都是上面那个很长的式子。求对于每个位置 $\mathbf U$ 的每个 $t\in[0,k-1]$ ,在这 $n$ 个数中有多少个数对于他的贡献是 $\omega_k^t$

这东西看上去异常的难以下手,出题人的标程里写的好像是什么奇奇怪怪的折半引理消去引理的奇淫技巧看不懂

讨论区里又有一位 djq 大爷说是什么单位根反演,但(道理我都懂板子我都抄过,哪里有单位根反演啊/晕)

然后我不知道怎么搞直接把每个值对应下标加一然后做一遍 DFT ,令人惊奇的发现它!就!对!了!

(我写到这里一看上面我写的,就想问自己,这**为什么不对啊)

但这里给读者讲一下(也给将来的我自己)(看懂的就直接跳过好了):

大概是对应下标加一后位置 $\mathbf V$ 的系数就是这个 $\mathbf{V}$ 出现的次数,对于每一个 $\mathbf V$ 它对于每个下标 $\mathbf U$ 的影响都是个定值 $\omega_k^t$ 。那么有 $a_{\mathbf V}$ 个,贡献就是 $a_{\mathbf V}\omega_k^t$ 。那么我们直接把这 $a_V$ 个一起做 DFT 肯定是没有问题的(我们每次是 $a_{\mathbf V}$ 加一(而没有改变后面的 $\omega_k^1,\omega_k^2\cdots $ 的项的系数 ),所以最后乘出来对于每个位置的贡献确实就是 $\omega_k^t$ 有多少个算多少个)感觉跟没讲一样,有啥问题评论区问吧

现在我们知道了对于一个位置 $\mathbf U$ 和 $t\in[0,k-1]$ 有多少个数最后对其的贡献是 $\omega_k^t$ (记为 $a(\mathbf V,t)$ 。那么这些数点积时都会是 $1+\omega_k^t$ 的形式,那么最后点积出来位置 $\mathbf V$ 的值就是 $\prod_{t=0}^{k-1}(1+\omega_k^t)^{a(\mathbf V,t)}$

最后再做一遍 IDFT 就能得到最终的答案了,总复杂度大概是 $O(k^{m+2}\sim k^{m+3})$ 左右,如果快速幂没有预处理(见实现细节)可能还要多个 $\log n$

好你就这样高高兴兴的写完然后输出每个答案 ( $a_0\omega_k^0+a_1\omega_k^1+\cdots+a_{k-1}\omega_k^{k-1}$ ) 中的 $a_0$ 就发现样例都过不了???

你一看题解发现要输出 $a_0+a_1-a_2-a_3$ 才行???而且写了就直接能过???

如果你做过 cf1103E 的话可能也会对洛谷上唯一的题解里最终输出的是 a[i].a[0]-a[i].a[4] 而不是 a[i].a[0] 感到疑惑(然而题解并没有提到这个 magic number 是哪来的(以及他把 $\omega_{10}$ 写成 $-\omega_5^3$ 为啥是对的都没有解释,像我等蒟蒻看的真是难受))

这里就方便一下我等蒟蒻略作证明(不需要的可以跳过看代码了):

$k=5$

先说结论:可以证得:$a_1=a_2=a_3=a_4$

(可以像那篇题解一样「带进去检验」,然而对于我等来说真不容易)

我的证法是:你可以发现最终这个东西一定是个整数,并且无论我们把 $\omega_5^1,\omega_5^2,\omega_5^3,\omega_5^4$ 中的哪一个带入答案都会是整数(这个我确实不会证,应该和 FFT 最后会是不带虚部的实数是同一个原因)

那么我们在复平面中画出所有的 $\omega_5$ 的幂

pic1

我们最后要得到 $a_0\omega_5^0+a_1\omega_5^1+a_2\omega_5^2+a_3\omega_5^3+a_4\omega_5^4$ 是一个不带虚部的整数

而 $\sin 18^{\circ}=\frac{\sqrt 5-1}4,\cos 18^{\circ}=\frac{\sqrt{10+2\sqrt 5}}4,\sin 36^{\circ}=\frac{\sqrt{10-2\sqrt 5}}4,\cos 36^{\circ}=\frac{1+\sqrt 5}4$

所以得到

又因为 $a_0,a_1,a_2,a_3,a_4\in \mathbf Z$ ,由上面的第二个式子得到 $a_1-a_4=a_2-a_3=0$ ,又由第一个式子得到 $a_1+a_4-a_2-a_3=0$ 于是就有 $a_1=a_2=a_3=a_4$

实际上不严谨的证明大概想一下上下对应的系数一定相同才能使得最后 $y$ 被抵消完落在 $x$ 轴上成为不带虚部的整数,左右两边如果相同的话就是 $\cos 72^{\circ}-\cos 36^{\circ}=-\frac 12$ 是个有理数 感觉是对的

知道这个后求值就很简单了:

当然也可以直接用刚才三角函数的式子: $x=\frac {\sqrt 5 }4(a_1+a_4-a_2-a_3)+\frac 14(4a_0-a_1-a_4-a_2-a_3)=a_0-a_1$

$k=6$

也差不多,不知道其他人怎么证的我还是画了个图:

pic2

和刚才类似的,我们可以得到

$\sin 60^{\circ}=\frac{\sqrt 3}2,\cos 60^{\circ}=\frac 12$ ,我们得到 $a_1+a_2-a_4-a_5=0$ 以及 $x=a_0-a_3+\frac 12(a_1-a_2-a_4+a_5)\in \mathbf Z$

这里并不能保证上下对应位置的系数相等了(即不能保证 $a_1=a_2=a_4=a_5$ ),但由 $a_1+a_2-a_4-a_5=0$ 可以得到 $a_1-a_4=a_5-a_2$ 得到 $x=a_0-a_3+\frac 12(a_1-a_4+a_5-a_2)=a_0-a_3+a_1-a_4$

答案就是 $a_0+a_1-a_4-a_3$ (题解里好像是 $a_0+a_1-a_2-a_3$ 所以应该可以证明 $a_2=a_4,a_1=a_5$ ?感觉我还不会证,欢迎交流。反正打表出来是对的那我们就当它是对的

一些实现上的细节

  1. 求 $1+\omega_k^t$ 在 $t=0$ 的时候是 $1+1=2$ 注意不要对应系数直接赋值成 $1$ 了
  2. $(1+\omega_k^t)^{a(\mathbf V,t)}$ 的时候快速幂可能被卡常,可以预处理出 $(1+\omega_k^t)^{1000}$ 的幂来加速
  3. $k=5$ 也可以用 $a_0+a_1-a_2-a_3$ 表示( $=a_0-a_1$ )不用分类讨论
  4. DFT 时的单位根只有一项有值,可以把 $O(k^2)$ 的乘法优化成 $O(k)$ 我太懒了就没写

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
#include<bits/stdc++.h>
#define FIO "P5577"
#define ll long long

char buf[1<<20]; int bufl,bufr;
#define getch ((bufl^bufr||(bufl=0,bufr=fread(buf,1,1<<20,stdin)))?buf[bufl++]:EOF)
template <class T>inline void read (T &x,const int &bas=10) {
T f=1;
x=0;
char ch=getch;
for (; !isdigit (ch)&&ch!='-'; ch=getch);
if (ch=='-')f=-1,ch=getch;
for (; isdigit (ch); ch=getch)x=x*bas+ch-'0';
x*=f;
}

const int MOD=998244353,V=6,N=pow(5,7)+5;

inline int add(int a,const int &b){a+=b;return a>=MOD?a-MOD:a;}
inline int sub(int a,const int &b){a-=b;return a< 0?a+MOD:a;}
inline int mul(const int &a,const int &b){return 1ll*a*b%MOD;}
inline int& inc(int &a,const int &b){return a=add(a,b);}
inline int& dec(int &a,const int &b){return a=sub(a,b);}
inline int& pro(int &a,const int &b){return a=mul(a,b);}
inline int qpow(int a,int b){int c=1;for(;b;b>>=1,pro(a,a))if(b&1)pro(c,a);return c;}

int n,k,m;

struct Com{
int a[V];
inline void clear(){a[0]=a[1]=a[2]=a[3]=a[4]=a[5]=0;}
inline Com(){clear();}
inline int& operator [](int x){return a[x];}
inline const int& operator [](const int &x)const{return a[x];}
inline Com operator +(const Com &b)const{Com c;for(int i=0;i<k;i++)c[i]=add(a[i],b[i]);return c;}
inline Com operator -(const Com &b)const{Com c;for(int i=0;i<k;i++)c[i]=sub(a[i],b[i]);return c;}
inline Com operator *(const int &b)const{Com c;for(int i=0;i<k;i++)c[i]=mul(a[i],b);return c;}
inline Com& operator +=(const Com &b){return *this=*this+b;}
inline Com& operator *=(const int &b){return *this=*this*b;}
inline Com operator *(const Com &b)const{
static ll r[V<<1];
const ll Moc=8ll*MOD*MOD;
for(int i=0;i<k;i++)if(a[i])for(int j=0;j<k;j++)if(b[j]){
r[i+j]+=1ll*a[i]*b[j];
if(r[i+j]>=Moc)r[i+j]-=Moc;
}
Com c;
for(int i=0;i<k;i++)c[i]=(r[i]+r[i+k])%MOD,r[i]=r[i+k]=0;
return c;
}
inline Com& operator *=(const Com &b){return *this=*this*b;}
inline void out(){for(int i=0;i<k;i++)printf("%d%c",a[i],i^(k-1)?' ':'\n');}
inline int val(){
return sub(add(a[0],a[1]),add(a[2],a[3]));
}
}a[N],w[2][V],b[N];

inline void pre(){
w[0][0][0]=w[1][0][0]=1;
w[0][1][1]=w[1][1][k-1]=1;
for(int i=2;i<k;i++)
w[0][i]=w[0][i-1]*w[0][1],w[1][i]=w[1][i-1]*w[1][1];
}

//inline Com qpow(Com a,int b){Com c;c[0]=1;for(;b;b>>=1,a*=a)if(b&1)c*=a;return c;}
const int B=1000;
Com pw0[B+5],pw1[B+5];
inline Com qpow(Com a,int b){return pw1[b/B]*pw0[b%B];}

inline void getpow(){
for(int i=0;i<n;i++)b[i][0]=1;
for(int i=0;i<k;i++){
Com omega;
//!!! i=1 gg
omega[0]=1;omega[i]++;
pw0[0][0]=pw1[0][0]=1;
for(int i=1;i<=B;i++)pw0[i]=pw0[i-1]*omega;
pw1[1]=pw0[B];
for(int i=1;i<=B;i++)pw1[i]=pw1[i-1]*pw1[1];
for(int id=0;id<n;id++)b[id]*=qpow(omega,a[id][i]);
}
for(int i=0;i<n;i++)a[i]=b[i];
}

int pw[V];
inline void fft(int opt){
pw[0]=1;for(int i=1;i<m;i++)pw[i]=pw[i-1]*k;
for(int i=0;i<m;i++){
for(int j=0;j<n;j++)
if(!((j/pw[i])%k)){
static Com f[V];
for(int x=0;x<k;x++)
for(int y=0;y<k;y++)
f[y]+=a[j+pw[i]*x]*w[opt][x*y%k];
for(int x=0;x<k;x++)
a[j+pw[i]*x]=f[x],f[x].clear();
}
}
if(opt)for(int i=0,inv=qpow(n,MOD-2);i<n;i++)a[i]*=inv;
}

int main(){
freopen(FIO".in","r",stdin);
freopen(FIO".out","w",stdout);
read(n);read(k);read(m);
for(int i=1,x;i<=n;i++)read(x,k),a[x][0]++;
n=qpow(k,m);
pre();

fft(0);

getpow();

fft(1);

for(int i=0;i<n;i++)printf("%d\n",a[i].val());

return 0;
}