题解 P7358 【窗花】
do_while_true
2021-02-15 21:30:01
[$\text{获得另一种阅读体验}$](https://www.cnblogs.com/do-while-true/p/14405475.html)
设赢，平局，输的概率分别为 $P_1,P_2,P_3$。
设 $E(i)$ 为从计数器为 $i$ 开始走，走到 $m$ 的步数。根据定义有 $E(m)=0$。
据题意，得出
$$
E(i)=
\begin{cases}
P_1E(i+1)+P_2E(i)+P_3E(i-1)+1,(i\geq 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)$
$$
\because E(0)=P_1E(1)+(1-P_1)E(0)+1
\\
\therefore E(1)=E(0)-\frac{1}{P_1}
\\
\therefore 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:$
$$
\\
\text{初始矩阵：}\begin{bmatrix}1&1\end{bmatrix}
\\
\text{转移矩阵：}\begin{bmatrix}0&t_1\\1&t_2\end{bmatrix}
$$
$b:$
$$
\\
\text{初始矩阵：}\begin{bmatrix}0&t_3&1\end{bmatrix}
\\
\text{转移矩阵：}
\begin{bmatrix}
0&t_2&0
\\
1&t_1&0
\\
0&t_3&1
\end{bmatrix}
$$
$m$ 可以高精除，模拟短除法得到 $m$ 的二进制，遍历一遍 $m$ 的二进制上各个位即可完成快速幂。
时间复杂度加上矩阵的常数为 $\mathcal{O}(\lg m\times\log_2 m+2^3\times\log_2 m+3^3\times\log_2 m)$，前面是求 $m$ 的二进制，后面两个分别为 $a,b$ 的矩阵快速幂求解。
即使本文略去简单的计算及推导，也难免会出现错误，欢迎指正。
$\mathcal{Code}$
```cpp
#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\n", 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;
}
```