题解 P7358 【窗花】

· · 题解

\text{获得另一种阅读体验}

设赢，平局，输的概率分别为 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} 时间复杂度加上矩阵的常数为 $\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} 
#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;
}