Luogu P7358 窗花

( ext{获得另一种阅读体验})

设赢,平局,输的概率分别为 (P_1,P_2,P_3)

(E(i)) 为从计数器为 (i) 开始走,走到 (m) 的步数。根据定义有 (E(m)=0)

据题意,得出

[E(i)= egin{cases} P_1E(i+1)+P_2E(i)+P_3E(i-1)+1,(igeq 1) \ P_1E(1)+(1-P_1)E(0)+1,(i==0) end{cases} ]

(E(0)=x)

(E(i)=a_ix+b_i,(i>1))

[ecause E(0)=P_1E(1)+(1-P_1)E(0)+1 \ herefore E(1)=E(0)-frac{1}{P_1} \ herefore a_1=1,b_1=-frac{1}{P_1} ]

同理可得:

[a_i=frac{1-P_2}{P_1}a_{i-1}-frac{P_3}{P_1}a_{i-2},(i>1) \ b_i=frac{1-P_2}{P_1}b_{i-1}-frac{P_3}{P_1}b_{i-2}-frac{1}{P_1},(i>1) ]

推出 (a_m,b_m) 后可得一个一元一次方程

[E(m)=a_m x+b_m=0 \ E(0)=x=-frac{b_m}{a_m} ]

由定义可知 (E(0)) 即为所求答案。

如何快速求得 (a_m,b_m)?可以矩阵快速幂。

这里就是朴素的矩阵快速幂了,这里给出我推得的初始矩阵和转移矩阵:

(t_1=frac{1-P_2}{P_1}a_{i-1},t_2=-frac{P_3}{P_1}b_{i-2},t_3=-frac{1}{P_1})

(a:)

[\ ext{初始矩阵:}egin{bmatrix}1&1end{bmatrix} \ ext{转移矩阵:}egin{bmatrix}0&t_1\1&t_2end{bmatrix} ]

(b:)

[\ ext{初始矩阵:}egin{bmatrix}0&t_3&1end{bmatrix} \ ext{转移矩阵:} egin{bmatrix} 0&t_2&0 \ 1&t_1&0 \ 0&t_3&1 end{bmatrix} ]

(m) 可以高精除,模拟短除法得到 (m) 的二进制,遍历一遍 (m) 的二进制上各个位即可完成快速幂。

时间复杂度加上矩阵的常数为 (mathcal{O}(lg m imeslog_2 m+2^3 imeslog_2 m+3^3 imeslog_2 m)),前面是求 (m) 的二进制,后面两个分别为 (a,b) 的矩阵快速幂求解。

即使本文略去简单的计算及推导,也难免会出现错误,欢迎指正。

(mathcal{Code})

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
namespace do_while_true {	
	#define ld double
	#define ll long long
	#define re register
	#define pb push_back
	#define fir first
	#define sec second
	#define pp std::pair
	#define mp std::make_pair
	const ll mod = 1000000007;
	template <typename T> inline T Max(T x, T y) { return x > y ? x : y; }
	template <typename T> inline T Min(T x, T y) { return x < y ? x : y; }
	template <typename T> inline T Abs(T x) { return x < 0 ? -x : x; }
	#define p mod 
	template <typename T> inline T Add(T x, T y) { return (x + y) % p; }
	template <typename T> inline T Mul(T x, T y) { return (x * y) % p; } 
	template <typename T> inline T Mod(T x) { return x % mod; }
	#undef p
	template <typename T>
	inline T& read(T& r) {
		r = 0; bool w = 0; char ch = getchar();
		while(ch < '0' || ch > '9') w = ch == '-' ? 1 : 0, ch = getchar();
		while(ch >= '0' && ch <= '9') r = r * 10 + (ch ^ 48), ch = getchar();
		return r = w ? -r : r;
	}
	template <typename T>
	inline T qpow(T x, T y) {
		re T sumq = 1; x %= mod;
		while(y) {
			if(y&1) sumq = sumq * x % mod;
			x = x * x % mod;
			y >>= 1;
		}
		return sumq;
	}
	template <typename T> inline T Inv(T x) { return qpow(x, mod-2); }
	char outch[110];
	int outct;
	template <typename T>
	inline void print(T x) {
		do {
			outch[++outct] = x % 10 + '0';
			x /= 10;
		} while(x);
		while(outct >= 1) putchar(outch[outct--]);
	}
}
using namespace do_while_true;

const int N = 1000100;
int n;
ll P1, P2, P3;
ll a[N], b[N], suma[N], sumb[N], isuma, isumb;
ll c[N], d[N];
char ch[N];
int _m[N], m[N], len, ct;

struct Matrix {
	int n, m;
	ll a[5][5];
	void mem() {
		n = m = 0;
		for(int i = 0; i <= 4; ++i)
			for(int j = 0; j <= 4; ++j)
				a[i][j] = 0;
	}
	Matrix operator * (const Matrix& y) {
		Matrix c; c.mem();
		c.n = n; c.m = y.m;
		for(int i = 1; i <= c.n; ++i)
			for(int j = 1; j <= c.m; ++j)
				for(int k = 1; k <= m; ++k)
					(c.a[i][j] += a[i][k] * y.a[k][j]) %= mod;
		return c;
	}
};

void solve() {
	scanf("%d %s", &n, ch+1); len = std::strlen(ch+1);
	for(int i = 1; i <= len; ++i) _m[i] = ch[len-i+1]-'0';
	while(len>=0 && ++ct) { 
		if(_m[1] & 1) m[ct] = 1;
		for(int i = len, x = 0; i >= 1; --i) _m[i] += x*10, x = _m[i]%2, _m[i] /= 2;
		while(_m[len] == 0 && len >= 0) --len;
	}
	for(int i = 1; i <= n; ++i) read(a[i]), suma[i] = Add(suma[i-1], a[i]);
	for(int i = 1; i <= n; ++i) read(b[i]), sumb[i] = Add(sumb[i-1], b[i]);
	isuma = qpow(suma[n], mod-2); isumb = qpow(sumb[n], mod-2);
	for(int i = 1; i <= n; ++i) {
		P1 = Add(P1, Mul(a[i] * isuma % mod, sumb[i-1] * isumb % mod));
		P2 = Add(P2, Mul(a[i] * isuma % mod, b[i] * isumb % mod));
		P3 = Add(P3, Mul(a[i] * isuma % mod, Mod(sumb[n] - sumb[i] + mod) * isumb % mod));
	}
	c[0] = 1; d[0] = 0; c[1] = 1; d[1] = Mod(-Inv(P1)+mod);
	ll t1 = Mod(1 - P2 + mod) * Inv(P1) % mod;
	ll t2 = Mod(-P3 + mod) * Inv(P1) % mod;
	ll t3 = d[1];
	Matrix basea, ansa, baseb, ansb;
	basea.mem(); ansa.mem(); baseb.mem(); ansb.mem();
	ansa.n = 1; ansa.m = 2;
	ansa.a[1][1] = 1; ansa.a[1][2] = 1;
	basea.n = 2; basea.m = 2;
	basea.a[1][1] = 0; basea.a[1][2] = t1;
	basea.a[2][1] = 1; basea.a[2][2] = t2;
	for(int i = 1; i <= ct; ++i, basea = basea * basea) if(m[i]) ansa = ansa * basea; 
	ansb.n = 1; ansb.m = 3;
	ansb.a[1][1] = 0; ansb.a[1][2] = t3; ansb.a[1][3] = 1;
	baseb.n = 3; baseb.m = 3;
	baseb.a[1][1] = 0; baseb.a[1][2] = t2; baseb.a[1][3] = 0;
	baseb.a[2][1] = 1; baseb.a[2][2] = t1; baseb.a[2][3] = 0;
	baseb.a[3][1] = 0; baseb.a[3][2] = t3; baseb.a[3][3] = 1;
	for(int i = 1; i <= ct; ++i, baseb = baseb * baseb) if(m[i]) ansb = ansb * baseb; 
	printf("%lld
", Mod(-ansb.a[1][1]+mod) * Inv(ansa.a[1][1]) % mod);
}

signed main() {
	#ifndef ONLINE_JUDGE
	freopen("in.txt", "r", stdin);
	#endif
	int T = 1;
//	read(T);
	while(T--) solve();
	fclose(stdin);
	return 0;
}
原文地址:https://www.cnblogs.com/do-while-true/p/14405475.html