点这里看题目。
在此之前,约定 \(x=(x_{m-1}x_{m-2}x_{m-3}\dots x_0)_k\),后者为 \(x\) 的 \(k\) 进制表示(从高位到低位),且根据题意至多 \(m\) 位。
约定 \(\omega_k\) 为 \(k\) 次单位根。
确实是大大地考验了我的算力。
显然,我们需要计算如下形式的卷积:
\[\bigoplus_{j=1}^n(1+x^{A_j}) \]在这里,我们不妨认为 \(\bigoplus\) 表示的是对于指数的 K 进制异或卷积。
如果你做过类似的题目,比如「UNR #2」黎明前的巧克力,那么应该记得,这种题目有一个普适的方法,就是直接计算答案序列 FWT 之后的结果,然后一遍 IFWT 回去。
这道题也适用于这个思路。首先考察一个序列的 FWT:
\[\operatorname{FWT}(1+x^{a})[x]=1+\prod_{j=0}^{m-1}\omega_k^{a_jx_j} \]所以,一个序列 FWT 之后的结果应当只有 \(k\) 种,形如 \(1+\omega_k^j\)。显然,我们可以尝试计算在第 \(x\) 位上,\(1+\omega_k^j\) 这样的因子出现了多少次,则最后可以直接算幂次得到结果。所以,不妨设 \(c_{x,j}\) 表示在第 \(x\) 位上,形如 \(1+\omega_k^j\) 因子出现的次数。
注意到,我们已然知道 \(\sum_j c_{x,j}=n\),这暗示我们可以通过解方程得到各个 \(c\)。尝试利用 FWT 本身来建立一个方程:注意到,根据 FWT 运算的线性性,我们很容易得到 \(\sum_{j} c_{x,j}\omega_k^j\),只需要考察 \(\operatorname{FWT}(\sum_{j=1}^nx^{A_j})[x]\) 即可。
再看一眼已有的两个方程:
\[\begin{cases} \sum_{j=0}^{k-1}c_{x,j}&=n\\ \sum_{j=0}^{k-1}c_{x,j}\omega_k^j&=\dots \end{cases} \]注意到,方程的系数已经有点像 DFT 矩阵的系数了。因此,我们尝试补全这个矩阵。由于 FWT 实际上就是高维 DFT,所以我们不难修改它,从而得到其余的几条方程。实际上,只需要改一下插入的值——无非就是将 \(\omega_k\) 改成 \(\omega_k^j\)——即可生成一条新的方程。
下面先考虑如何实现这个算法。由于 998244353 性质(在这里)不够好,或者 \(k=5,6\) 实在是两个精心构造的值,反正我们不能直接用原根的幂次替代 \(\omega_k\)。但是,在这种背景下,绝对不能忘记扩域这样一个选项——直接引入 \(\omega_k\) 及其幂次即可。
但是,又会有新的问题出现:由于单位根的某些性质,比如 \(\sum_{j=0}^{k-1}\omega_k^j=0\),一个整数在扩张后的域内不能唯一表示。当然,放到几何背景下,这其实就说明某些(我们认为是基向量的)向量可以被其它向量线性表示。解决办法也比较简单,消去多余的元(在复平面上画图观察)即可。
最终的结果是:在 \(k=5\) 的时候,我们只需要留下 \(\omega_5^0,\omega_5^1,\omega_5^2,\omega_5^3\);而在 \(k=6\) 的时候,我们只需要留下 \(\omega_k^0,\omega_k^2\)。
最后是一个观察——当我们试图生成下式:
\[\sum_{j=0}^{k-1}c_{x,j}\omega_k^{vj} \]的值的时候,我们实际上是在计算:
\[\sum_{i=1}^n\prod_{j=0}^{m-1}(\omega_k^v)^{x_j(a_i)_j} \]如果每次都换一下插入的值实在是太浪费了。注意到 \((\omega_k^v)^{x_j}=\omega_k^{vx_j}\),所以我们可以只需要插入 \(\omega_k\) 计算一次,就可以通过下标移位找到需要的值,从而降一个 \(k\)。
注意,这里不能把插入 \(\omega_k^v\) 看作是“复合”上一个新的元,因为外部的运算并不会满足真正的“复合”的逻辑。
最终复杂度大约为 \(O(k^{m+2}m)\)。
小结:
扩域的想法,或者说是解决问题的手段,要积累一下。
简单序列 FWT 的结果可以直接推导,这个是一个常见的思路;
考察不同的项的指数,通过解方程得到指数,也算是一个常用的思路,当然关键在于如何建立方程。
最后一步的观察,思路比较朴素,也就是要充分利用计算出来了的信息。支撑这个方法的还有一点,就是整个下标集在按位运算下是封闭的。
#include <cmath> #include <cstdio> #define rep( i, a, b ) for( int i = (a) ; i <= (b) ; i ++ ) #define per( i, a, b ) for( int i = (a) ; i >= (b) ; i -- ) const int mod = 998244353, g = 3, phi = mod - 1; const int MAXN = 1e5 + 5, MAXRT = 1e3 + 5; template<typename _T> void read( _T &x ) { x = 0; char s = getchar(); bool f = false; while( ! ( '0' <= s && s <= '9' ) ) { f = s == '-', s = getchar(); } while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); } if( f ) x = -x; } template<typename _T> void write( _T x ) { if( x < 0 ) putchar( '-' ), x = -x; if( 9 < x ) write( x / 10 ); putchar( x % 10 + '0' ); } template<typename _T> struct SwiftPow { _T sml[MAXRT], lrg[MAXRT]; int B; SwiftPow(): B( 1 ) {} inline void Init( const _T &x, const int &n ) { B = ceil( sqrt( n ) ); lrg[0] = sml[0] = _T( 1 ); rep( i, 1, B ) sml[i] = sml[i - 1] * x; lrg[1] = sml[B]; rep( i, 2, B ) lrg[i] = lrg[i - 1] * lrg[1]; } inline _T operator () ( const int &n ) const { return lrg[n / B] * sml[n % B]; } }; int app[MAXN]; int N, K, M; inline int Qkpow( int, int ); inline int Inv( const int &a ) { return Qkpow( a, mod - 2 ); } inline int Mul( int x, const int &v ) { return 1ll * x * v % mod; } inline int Sub( int x, const int &v ) { return ( x -= v ) < 0 ? x + mod : x; } inline int Add( int x, const int &v ) { return ( x += v ) >= mod ? x - mod : x; } inline int& MulEq( int &x, const int &v ) { return x = 1ll * x * v % mod; } inline int& SubEq( int &x, const int &v ) { return ( x -= v ) < 0 ? ( x += mod ) : x; } inline int& AddEq( int &x, const int &v ) { return ( x += v ) >= mod ? ( x -= mod ) : x; } inline int Qkpow( int base, int indx ) { int ret = 1; while( indx ) { if( indx & 1 ) MulEq( ret, base ); MulEq( base, base ), indx >>= 1; } return ret; } inline int readK() { static char buf[20]; int ret = 0; scanf( "%s", buf ); for( int i = 0 ; buf[i] ; i ++ ) ret = ret * K + ( int ) ( buf[i] - '0' ); return ret; } namespace KEq5 { struct Complex { int a[4]; Complex(): a{} {} Complex( int R ): a{} { a[0] = R; } inline int& operator [] ( const int &idx ) { return a[idx]; } inline const int operator [] ( const int &idx ) const { return a[idx]; } inline Complex& operator += ( const Complex &b ) { rep( i, 0, 3 ) AddEq( a[i], b[i] ); return *this; } inline Complex& operator -= ( const Complex &b ) { rep( i, 0, 3 ) SubEq( a[i], b[i] ); return *this; } inline Complex& operator /= ( int b ) { b = Inv( b ); rep( i, 0, 3 ) MulEq( a[i], b ); return *this; } inline Complex& operator *= ( const Complex &b ) { return *this = *this * b; } inline Complex operator + ( const Complex &b ) const { Complex ret; rep( i, 0, 3 ) ret[i] = Add( a[i], b[i] ); return ret; } inline Complex operator - ( const Complex &b ) const { Complex ret; rep( i, 0, 3 ) ret[i] = Sub( a[i], b[i] ); return ret; } inline Complex operator * ( const Complex &b ) const { Complex ret; int out = 0; rep( i, 0, 3 ) if( a[i] ) rep( j, 0, 3 ) if( b[j] ) { if( ( i + j ) % 5 == 4 ) AddEq( out, Mul( a[i], b[j] ) ); else AddEq( ret[( i + j ) % 5], Mul( a[i], b[j] ) ); } rep( i, 0, 3 ) SubEq( ret[i], out ); return ret; } inline Complex operator / ( int b ) const { Complex ret; b = Inv( b ); rep( i, 0, 3 ) ret[i] = Mul( a[i], b ); return ret; } }; Complex coe[MAXN], nw[MAXN], w[5], res[5]; SwiftPow<Complex> pw[5]; int lim; Complex Work( const int x ) { static int val[5] = {}; static int inv5 = Inv( 5 ); Complex f = 0, ret = 1; rep( i, 0, K - 1 ) { val[i] = 0; for( int j = 0, s = 1 ; j < M ; j ++, s *= K ) val[i] += s * ( ( x / s % 5 ) * i % 5 ); } rep( i, 0, K - 1 ) { f = 0; rep( j, 0, K - 1 ) f += coe[val[j]] * w[( K - i * j % K ) % K]; ret *= pw[i]( Mul( inv5, f.a[0] ) ); } return ret; } void Init() { lim = 1; rep( i, 1, M ) lim *= K; rep( i, 0, lim - 1 ) coe[i] = app[i]; rep( i, 0, 3 ) w[i].a[i] = 1, w[4].a[i] = mod - 1; rep( i, 0, 4 ) pw[i].Init( w[i] + 1, N ); } void Solve() { Init(); Complex tmp; for( int s = 1 ; s < lim ; s *= K ) for( int i = 0 ; i < lim ; i += s * K ) for( int j = i ; j < i + s ; j ++ ) { for( int k = 0 ; k < 5 ; k ++ ) { tmp = 0; for( int t = 0 ; t < 5 ; t ++ ) tmp += coe[j + t * s] * w[k * t % 5]; res[k] = tmp; } for( int k = 0 ; k < 5 ; k ++ ) coe[j + k * s] = res[k]; } rep( i, 0, lim - 1 ) nw[i] = Work( i ); rep( i, 0, lim - 1 ) coe[i] = nw[i]; for( int s = 1 ; s < lim ; s *= K ) for( int i = 0 ; i < lim ; i += s * K ) for( int j = i ; j < i + s ; j ++ ) { for( int k = 0 ; k < 5 ; k ++ ) { tmp = 0; for( int t = 0 ; t < 5 ; t ++ ) tmp += coe[j + t * s] * w[( 5 - k * t % 5 ) % 5]; res[k] = tmp / 5; } for( int k = 0 ; k < 5 ; k ++ ) coe[j + k * s] = res[k]; } rep( i, 0, lim - 1 ) write( coe[i].a[0] ), putchar( '\n' ); } } namespace KEq6 { struct Complex { int w0, w2; Complex(): w0( 0 ), w2( 0 ) {} Complex( int R ): w0( R ), w2( 0 ) {} Complex( int W0, int W2 ): w0( W0 ), w2( W2 ) {} inline Complex& operator += ( const Complex &b ) { AddEq( w0, b.w0 ); AddEq( w2, b.w2 ); return *this; } inline Complex& operator -= ( const Complex &b ) { SubEq( w0, b.w0 ); SubEq( w2, b.w2 ); return *this; } inline Complex& operator /= ( int b ) { b = Inv( b ); MulEq( w0, b ); MulEq( w2, b ); return *this; } inline Complex& operator *= ( const Complex &b ) { return *this = *this * b; } inline Complex operator + ( const Complex &b ) const { return Complex( Add( w0, b.w0 ), Add( w2, b.w2 ) ); } inline Complex operator - ( const Complex &b ) const { return Complex( Sub( w0, b.w0 ), Sub( w2, b.w2 ) ); } inline Complex operator * ( const Complex &b ) const { return Complex( Sub( Mul( w0, b.w0 ), Mul( w2, b.w2 ) ), Sub( Add( Mul( w2, b.w0 ), Mul( w0, b.w2 ) ), Mul( w2, b.w2 ) ) ); } inline Complex operator / ( int b ) const { b = Inv( b ); return Complex( Mul( w0, b ), Mul( w2, b ) ); } }; Complex nw[MAXN], coe[MAXN], w[6], res[6]; SwiftPow<Complex> pw[6]; int lim; Complex Work( const int x ) { static int val[6] = {}; static int inv6 = Inv( 6 ); Complex f = 0, ret = 1; rep( i, 0, K - 1 ) { val[i] = 0; for( int j = 0, s = 1 ; j < M ; j ++, s *= K ) val[i] += s * ( ( x / s % K ) * i % K ); } rep( i, 0, K - 1 ) { f = 0; rep( j, 0, K - 1 ) f += coe[val[j]] * w[( K - i * j % K ) % K]; ret *= pw[i]( Mul( inv6, f.w0 ) ); } return ret; } void Init() { lim = 1; rep( i, 1, M ) lim *= K; rep( i, 0, lim - 1 ) coe[i] = app[i]; w[0] = Complex( 1, 0 ); w[1] = Complex( 1, 1 ); w[2] = Complex( 0, 1 ); w[3] = Complex( mod - 1, 0 ); w[4] = Complex( mod - 1, mod - 1 ); w[5] = Complex( 0, mod - 1 ); rep( i, 0, 5 ) pw[i].Init( w[i] + 1, N ); } void Solve() { Init(); Complex tmp; for( int s = 1 ; s < lim ; s *= K ) for( int i = 0 ; i < lim ; i += s * K ) for( int j = i ; j < i + s ; j ++ ) { for( int k = 0 ; k < 6 ; k ++ ) { tmp = 0; for( int t = 0 ; t < 6 ; t ++ ) tmp += coe[j + t * s] * w[k * t % 6]; res[k] = tmp; } for( int k = 0 ; k < 6 ; k ++ ) coe[j + k * s] = res[k]; } rep( i, 0, lim - 1 ) nw[i] = Work( i ); rep( i, 0, lim - 1 ) coe[i] = nw[i]; for( int s = 1 ; s < lim ; s *= K ) for( int i = 0 ; i < lim ; i += s * K ) for( int j = i ; j < i + s ; j ++ ) { for( int k = 0 ; k < 6 ; k ++ ) { tmp = 0; for( int t = 0 ; t < 6 ; t ++ ) tmp += coe[j + t * s] * w[( 6 - k * t % 6 ) % 6]; res[k] = tmp / 6; } for( int k = 0 ; k < 6 ; k ++ ) coe[j + k * s] = res[k]; } rep( i, 0, lim - 1 ) write( coe[i].w0 ), putchar( '\n' ); } } int main() { read( N ), read( K ), read( M ); rep( i, 1, N ) AddEq( app[readK()], 1 ); if( K == 5 ) KEq5 :: Solve(); if( K == 6 ) KEq6 :: Solve(); return 0; }