其他分享
首页 > 其他分享> > 「CmdOI2019」算力训练

「CmdOI2019」算力训练

作者:互联网

题目

点这里看题目。

分析

在此之前,约定 \(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)\)。

小结:

  1. 扩域的想法,或者说是解决问题的手段,要积累一下。

  2. 简单序列 FWT 的结果可以直接推导,这个是一个常见的思路;

    考察不同的项的指数,通过解方程得到指数,也算是一个常用的思路,当然关键在于如何建立方程。

  3. 最后一步的观察,思路比较朴素,也就是要充分利用计算出来了的信息。支撑这个方法的还有一点,就是整个下标集在按位运算下是封闭的

代码

#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;
}

标签:const,CmdOI2019,训练,int,rep,Complex,inline,return,算力
来源: https://www.cnblogs.com/crashed/p/16226540.html