[ZJOI2019]开关

题目

点这里看题目。

分析

不妨设 (P=sum_i p_i)

简单的情况是,如果我们去掉第一次为目标状态的限制,那么可以设 (g_n) 为操作 (n) 次后达成目标状态的概率,就得到了一个很好计算的东西。不难发现这其实就是限制了按钮被操作次数,可以直接得出其指数型生成函数

[hat G(x)=prod_{i=1}^nfrac{e^{frac{p_i}{P{}}x}+(-1)^{s_i}e^{-frac{p_i}{P}x}}{2} ]

加上第一次的限制之后,我们可以知道两次目标状态之间操作的结果是 " 没有操作 " 。那操作结果为 " 没有操作 " 的指数型生成函数是什么?

[hat H(x)=prod_{i=1}^nfrac{e^{frac{p_i}{P}x}+e^{-frac{p_i}{P}x}}{2} ]

(h_n=n![x^n]hat H(x)) ,且 (f_n) 为操作 (n)第一次为目标状态的概率,那么可以得到:

[g_n=sum_{k=0}^nf_kh_{n-k} ]

(f,g,h)普通生成函数分别为 (F(x),G(x), H(x)) ,那么上式告诉我们 (G(x)=F(x)H(x)) ,所以有 (F(x)=frac{G(x)}{H(x)})


现在的推导出现了一点问题:怎么让 (hat G) 变成 (G)

注意 (hat G) 的形式,它必然为多项 (alpha e^{frac{eta}{P}x}) 之和。而一项 (alpha e^{frac{eta}Px}) 对应的普通生成函数为 (frac{alpha}{1-frac{eta}{P}x})

另一个限制是,显然 (-Ple etale P) ,所以我们可知:

[G(x)=sum_{eta=-P}^{P}frac{hat g_eta}{1-frac{eta}{P}x} ]

(hat H) 可做同样处理,且使用背包我们可以快速求出 (hat g)(hat h)

据说可以更快,不过在这里就没必要了

根据概率生成函数的知识,我们只需要求出 (F’(1)) 。根据除法求导法则有:

[F'(x)=left(frac{G(x)}{H(x)} ight)'=frac{G'(x)H(x)-H'(x)G(x)}{H^2(x)} ]

代入计算 (F'(1)) 就大功告成......了?

个锤子,出题人阴险得很

(G(x),H(x)) 里面居然都有 (frac{1}{1-x}) 的项,这代进去不就 GG 了?

这里的解决办法是,给 (G(x),H(x)) 同时乘上 ((1-x)) ,就可以规避这个问题。

写题解的人水平太低,解释不通。话说求导过程中, ((1-x)) 或许是可以直接约掉的。

对于新的函数而言,就有:

[egin{aligned} G(1)&=hat g_P\ H(1)&=hat h_P\ G'(1) &=sum_{eta=-P}^{P}hat g_etaleft.left(frac{1-x}{1-frac{eta}{P}x} ight)' ight|_{x=1}\ &=sum_{eta=-P}^{P}frac{hat g_{eta}}{frac{eta}{P}-1}\ H'(1)&=sum_{eta=-P}^Pfrac{hat h_eta}{frac{eta}{P}-1} end{aligned} ]

代入计算即可得到 (F'(1)) 。时间复杂度为 (O(nP))

小结:

  1. 容斥方法计算 (f) 比较常用,转成生成函数的表达更为精妙;
  2. 此类 OGF 转为 EGF 的方法好像很显
  3. 对于 (frac{1}{1-x}) 的处理很有意思;
  4. 求导法则一定好好背,微积分一定好好学

代码

#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, inv2 = 499122177;
const int MAXN = 105, MAXP = 5e4 + 5;

template<typename _T>
void read( _T &x )
{
	x = 0; char s = getchar(); int f = 1;
	while( s < '0' || '9' < s ) { f = 1; if( s == '-' ) f = -1; s = getchar(); }
	while( '0' <= s && s <= '9' ) { x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar(); }
	x *= f;
}

template<typename _T>
void write( _T x )
{
	if( x < 0 ) putchar( '-' ), x = - x;
	if( 9 < x ) write( x / 10 );
	putchar( x % 10 + '0' );
}

int BP[2][MAXP << 1], *B[2] = { BP[0] + MAXP, BP[1] + MAXP };
int AP[2][MAXP << 1], *A[2] = { AP[0] + MAXP, AP[1] + MAXP };

int S[MAXN], p[MAXN];
int N;

inline int Qkpow( int, int );
inline int Mul( int x, int v ) { return 1ll * x * v % mod; }
inline int Inv( const int a ) { return Qkpow( a, mod - 2 ); }
inline int Sub( int x, int v ) { return ( x -= v ) < 0 ? x + mod : x; }
inline int Add( int x, int v ) { return ( x += v ) >= mod ? x - mod : x; }
inline void Upt( int &x, const int v ) { x = Add( x, v ); }

int Qkpow( int base, int indx )
{
	int ret = 1;
	while( indx )
	{
		if( indx & 1 ) ret = Mul( ret, base );
		base = Mul( base, base ), indx >>= 1;
	}
	return ret;
}

int main()
{
	read( N ); int P = 0;
	rep( i, 1, N ) read( S[i] );
	rep( i, 1, N ) read( p[i] ), P = Add( P, p[i] );
	A[0][0] = B[0][0] = 1;
	int nxt = 0, pre = 1;
	rep( i, 1, N )
	{
		nxt ^= 1, pre ^= 1;
		int c = S[i] & 1 ? mod - inv2 : inv2;
		rep( j, - P, P ) A[nxt][j] = B[nxt][j] = 0;
		rep( j, - P, P )
		{
			if( A[pre][j] ) 
				Upt( A[nxt][j + p[i]], Mul( inv2, A[pre][j] ) ),
				Upt( A[nxt][j - p[i]], Mul( c, A[pre][j] ) );
			if( B[pre][j] )
				Upt( B[nxt][j + p[i]], Mul( inv2, B[pre][j] ) ),
				Upt( B[nxt][j - p[i]], Mul( inv2, B[pre][j] ) );
		}
	}
	int ans = 0, inv = Inv( P );
	rep( i, -P, P - 1 )
		ans = Add( ans, Mul( Sub( Mul( A[nxt][i], B[nxt][P] ), Mul( B[nxt][i], A[nxt][P] ) ),
			            Inv( Mul( Sub( Mul( i, inv ), 1 ), Mul( B[nxt][P], B[nxt][P] ) ) ) ) );
	write( ans ), putchar( '
' );
	return 0;
}
原文地址:https://www.cnblogs.com/crashed/p/14562456.html