标签: 概率期望

  • noi.ac 154 Curiosity

    $$ \text{设}f\left( i,j \right) \text{表示第}i\text{次扔骰子是否取到}j\text{这个数,显然这个东西的期望} \\ \text{是}\frac{1}{K}\left( \text{一次扔骰子只能出现一个数} \right) \\ \text{所以}a_i=\sum_{1\le k\le n}{f\left( k,i \right)} \\ \text{所以}E\left( \prod_{1\le i\le L}{a_i^F} \right) =\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F} \\ \text{考虑}\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F}\text{展开后的一项,如果这一项中存在}f\left( i,a \right) \text{和}f\left( i,b \right) ,\text{并且}a\ne b, \\ \text{那么这一整项的期望必须是}0\left( \text{否则出现了不合法情况} \right) \\ \text{对于其他的项,比如}f\left( \text{1,}a \right) ^2f\left( \text{2,}b \right) ^2….,\text{假设这一项由}p\text{个不同的变量构成} \\ \text{那么这一项的期望为}\frac{1}{K^p}\left( \text{独立变量期望可乘} \right) \\ \text{现在我们的问题在于如何计算出}\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F}\text{中由}p\text{个独立变量构成的项的个数。} \\ \text{首先考虑}\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F,\text{显然这个式子展开后不存在期望为0的项。} \\ \,\,\text{我们忽略项中是由哪些变量}\left( \text{比如}f\left( \text{1,}i \right) ,f\left( \text{2,}i \right) \text{等等} \right) \text{构成的,也忽略掉项之间的顺序关系} \\ \left( \begin{array}{c} \text{在}\left( f\left( \text{1,}i \right) +f\left( \text{2,}i \right) +f\left( \text{3,}i \right) .. \right) \times \left( f\left( \text{1,}i \right) +f\left( \text{2,}i \right) +f\left( \text{3,}i \right) .. \right) \times \left( f\left( \text{1,}i \right) +f\left( \text{2,}i \right) +f\left( \text{3,}i \right) .. \right) \text{中,}f\left( \text{1,}i \right) ^2\times f\left( \text{2,}i \right)\\ \text{会被计算6次}\\ \end{array} \right) \\ \text{则}\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F\text{展开后忽略上述限制存在}i\text{个无关变量的项的系数和为}S_2\left( F,i \right) ,\text{其中}S_2\text{为第二类斯特林数。} \\ \text{则}S_2\left( F,i \right) \times \left( \begin{array}{c} F\\ i\\ \end{array} \right) \times i!\text{为不忽略限制的符合条件的项的系数和。} \\ \text{构造生成函数}F\left( x \right) =\sum_{i\ge 0}{S_2\left( F,i \right) x^i} \\ \text{则}F\left( x \right) ^L\text{展开后}i\text{次项前的系数即为忽略上述限制后,存在}i\text{个无关变量的项的系数和。} \\ \text{为了还原至原式,将这个系数乘上}\left( \begin{array}{c} n\\ i\\ \end{array} \right) i!\text{后即为}\prod_{1\le i\le L}{\left( \sum_{1\le j\le N}{f\left( j,i \right)} \right) ^F}\text{展开后存在}i\text{个无关变量的项的系数和。} \\ \text{显然只有不超过}L\times F\text{个无关变量的项才是有意义的}. \\ \text{复杂度}O\left( LF\log LF\times \log L \right) \\ \text{注意}F=\text{1的时候需要特判,否则太慢了跑不过去。} \\ $$

    #include <assert.h>
    #include <algorithm>
    #include <cmath>
    // #include <complex>
    #include <cstring>
    #include <iostream>
    using real_t = double;
    
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    using cpx_t = struct complex;
    
    struct complex {
        real_t real, imag;
        complex(real_t real = 0, real_t imag = 0) {
            this->real = real;
            this->imag = imag;
        }
        complex operator+(const complex& rhs) const {
            return complex{real + rhs.real, imag + rhs.imag};
        }
        complex operator-(const complex& rhs) const {
            return complex{real - rhs.real, imag - rhs.imag};
        }
        complex operator*(const complex& rhs) const {
            return complex{real * rhs.real - imag * rhs.imag,
                           real * rhs.imag + imag * rhs.real};
        }
        complex& operator*=(const complex& rhs) {
            *this = (*this) * rhs;
            return *this;
        }
    };
    
    const int_t LARGE = 300000;
    const int_t mod = 2333;
    int_t S2[mod][mod];
    int_t fact[LARGE + 1], factInv[LARGE + 1], inv[LARGE + 1];
    const real_t PI = acosl(-1);
    int_t flips[LARGE + 1];
    void transform(cpx_t* A, int_t size2, int_t arg) {
        for (int_t i = 0; i < (1 << size2); i++) {
            int_t x = flips[i];
            if (x > i) std::swap(A[i], A[x]);
        }
        for (int i = 2; i <= (1 << size2); i *= 2) {
            auto mr = cpx_t(cos(2 * PI / i), arg * sin(2 * PI / i));
            for (int j = 0; j < (1 << size2); j += i) {
                auto curr = cpx_t(1, 0);
                for (int k = 0; k < i / 2; k++) {
                    auto u = A[j + k], t = A[j + k + i / 2] * curr;
                    A[j + k] = u + t;
                    A[j + k + i / 2] = u - t;
                    curr *= mr;
                }
            }
        }
    }
    void poly_mul(const int_t* Ax, const int_t* Bx, int_t size2, int_t n,
                  int_t* Cx) {
        static cpx_t A[LARGE * 2], B[LARGE * 2];
        for (int_t i = 0; i < (1 << size2); i++) {
            if (i < n) {
                A[i] = Ax[i];
                B[i] = Bx[i];
            } else {
                A[i] = B[i] = 0;
            }
        }
        transform(A, size2, 1);
        transform(B, size2, 1);
        for (int i = 0; i < (1 << size2); i++) A[i] *= B[i];
        transform(A, size2, -1);
        for (int_t i = 0; i < n; i++) {
            Cx[i] = (int_t((A[i]).real / (1 << size2) + 0.5) % mod);
            assert(Cx[i] >= 0);
        }
    }
    int_t process() {
        static int_t A[LARGE + 1], B[LARGE + 1];
    
        int_t n, k, l, f;
        cin >> n >> k >> l >> f;
        int_t size2 = 0;
        int_t nx = l * f + 1;
        while ((1 << size2) < (nx * 2)) size2++;
        for (int_t i = 0; i < (1 << size2); i++) {
            const auto flip = [&](int_t x) {
                int_t result = 0;
                for (int_t i = 0; i < size2 - 1; i++) {
                    result |= (x & 1);
                    x >>= 1;
                    result <<= 1;
                }
                return result | (x & 1);
            };
            flips[i] = flip(i);
        }
        std::fill(A, A + (1 << size2), 0);
        std::fill(B, B + (1 << size2), 0);
    
        if (k % mod == 0) return 0;
        for (int_t i = 1; i <= f; i++) {
            A[i] = S2[f][i];
        }
        B[0] = 1;
    #ifdef DEBUG
        cout << "size2 = " << size2 << " len = " << (1 << size2) << endl;
        cout << "nx  = " << nx << endl;
        cout << "A = ";
        for (int_t i = 0; i <= f; i++) cout << A[i] << " ";
        cout << endl;
    
    #endif
    
        if (f == 1) {
            B[l] = 1;
        } else {
            while (l) {
                assert(l >= 0);
                if (l & 1) poly_mul(A, B, size2, nx, B);
                poly_mul(A, A, size2, nx, A);
                l >>= 1;
            }
        }
    #ifdef DEBUG
        cout << "mul ok B = ";
        for (int_t i = 0; i < nx; i += 1) cout << B[i] << " ";
        cout << endl;
    
    #endif
        int_t result = 0;
        int_t kinv = inv[k % mod];
        int_t under = kinv;
        int_t prod = n;
        for (int_t i = 1; i < nx; i++) {
            result = (result + B[i] % mod * under * prod % mod) % mod;
    #ifdef DEBUG
            cout << "at " << i << " count on " << B[i] * under * prod % mod << endl;
    #endif
            under = under * kinv % mod;
            prod = prod * (n - i) % mod;
            if (n <= i) break;
        }
        return result;
    }
    int main() {
        S2[0][0] = 1;
        for (int_t i = 1; i < mod; i++) {
            for (int_t j = 1; j <= i; j++) {
                S2[i][j] = (S2[i - 1][j - 1] + S2[i - 1][j] * j) % mod;
            }
        }
        fact[0] = fact[1] = factInv[0] = factInv[1] = inv[1] = 1;
        for (int_t i = 2; i < mod; i++) {
            inv[i] = (mod - mod / i) * inv[mod % i] % mod;
            factInv[i] = factInv[i - 1] * inv[i] % mod;
            fact[i] = fact[i - 1] * i % mod;
        }
        int_t T;
        cin >> T;
        while (T--) cout << process() << endl;
    
        return 0;
    }

     

  • 洛谷3750 [六省联考2017]分手是祝愿

    $$ \text{首先可以证明,把原局面从到大小一个一个关所经历的步数一定最少} \\ \left( \text{每次必定选择编号最大的没有关闭的灯} \right) \text{所以先求出来这个步数。} \\ \text{如果这个步数不超过}k\text{,那么这个步数就是答案。} \\ \text{否则设}f\left( i \right) \text{表示还剩}i\text{步的时候,随机选择一个开关是转移到还剩}i-\text{1步的期望步数。} \\ f\left( i \right) =\frac{i}{n}\left( \text{选中了一个正确的开关} \right) +\left( 1-\frac{i}{n} \right) \left( f\left( i \right) +f\left( i+1 \right) +1 \right) \left( \text{选中了一个错误的开关,那么步数}+1 \right) \\ \text{显然}f\left( n \right) =\text{1,此时随便选一个开关都会让步数少}1. \\ \text{整理得}f\left( i \right) =1+\frac{n-i}{i}\left( f\left( i+1 \right) +1 \right) \\ \text{计算出}f\left( k+1 \right) \text{到}f\left( n \right) \text{的答案累加起来,加上}k\text{即可。} \\ $$

    #include <algorithm>
    #include <iostream>
    #include <vector>
    using int_t = long long int;
    using std::cin;
    using std::cout;
    using std::endl;
    
    const int_t mod = 100003;
    const int_t LARGE = 1e5;
    int_t inv[LARGE + 1];
    int_t state[LARGE + 1];
    int_t f[LARGE + 1];
    std::vector<int_t> divisors[LARGE + 1];
    int_t n, k;
    int main() {
        for (int_t i = 1; i <= LARGE; i++) {
            for (int_t j = 1; j * i <= LARGE; j++) {
                divisors[j * i].push_back(i);
            }
        }
        inv[1] = 1;
        for (int_t i = 2; i <= LARGE; i++) {
            inv[i] = (mod - mod / i) * inv[mod % i] % mod;
        }
        cin >> n >> k;
        for (int_t i = 1; i <= n; i++) {
            cin >> state[i];
        }
        int_t steps = 0;
        for (int_t i = n; i >= 1; i--) {
            if (state[i]) {
                steps++;
                for (int_t x : divisors[i]) state[x] ^= 1;
            }
        }
        int_t result = 0;
        if (steps <= k) {
            result = steps;
        } else {
            f[n] = 1;
            for (int_t i = n - 1; i >= 1; i--) {
                f[i] = (1 + (n - i) * inv[i] % mod * (f[i + 1] + 1) % mod) % mod;
            }
            for (int_t i = steps; i > k; i--) result = (result + f[i]) % mod;
            result = (result + k) % mod;
        }
        for (int_t i = 1; i <= n; i++) result = result * i % mod;
        cout << result << endl;
        return 0;
    }

     

  • SDOI2017 硬币游戏

    $$ \text{设}p\left( S \right) \text{表示字符串}S\text{出现的概率} \\ \text{设}P_i\text{为第}i\text{个模式串,}p_i\text{为第}i\text{个人赢的概率} \\ \text{对于第}i\text{个人,我们可以列方程} \\ \text{设}N_i\text{为不能使第}i\text{个人赢的字符串} \\ p\left( N_iP_i \right) =\sum_{1\le j\le n}{p_j\sum_{1\le l\le m}{\frac{\left[ Suffix\left( P_j,l \right) =Prefix\left( P_i,l \right) \right]}{2^{m-l}}}} \\ \text{即考虑所有能让第}i\text{个人赢的字符串,这个串一定是由一个不能使}i\text{赢的字符串加上一些其他东西得到的}. \\ \text{因为我们要列方程组,所以如果}P_j\text{的一个长度为}l\text{的后缀与}P_i\text{的一个前缀相等,那么使得}P_j\text{能赢的字符串加上若干字符}\left( m-l\text{个} \right) \\ \text{就可以使得}i\text{赢} \\ \text{同时,出现任意长度为}i\text{的字符串的概率是}2^{-i} \\ \text{同时}p\left( N_iP_i \right) =p\left( N_i \right) p\left( P_i \right) \\ \text{又因为所有}p\left( N_i \right) \text{相等} \\ \text{所以现在有了}n+\text{1个未知数}\left( p\left( N \right) \text{和}p_1,p_2….p_3 \right) \\ \text{算上}\sum_{1\le i\le n}{p_i}=\text{1后共计}n+\text{1个方程,高斯消元即可} $$

    #include <algorithm>
    #include <cstdio>
    #include <iomanip>
    #include <iostream>
    using int_t = long long int;
    using real_t = long double;
    using std::cin;
    using std::cout;
    using std::endl;
    const int_t LARGE = 600;
    const int_t SEED = 233;
    const int_t mod = 998244353;
    const real_t EPS = 1e-8;
    int_t prefix[LARGE + 1][LARGE + 1];
    real_t pown2[LARGE + 1] = {1};
    real_t mat[LARGE + 1][LARGE + 1];
    int_t pow_inv[LARGE + 1] = {1};
    int_t n, m;
    int_t power(int_t base, int_t index) {
        const auto phi = mod - 1;
        index = (index % phi + phi) % phi;
        base = (base % mod + mod) % mod;
        int_t result = 1;
        while (index) {
            if (index & 1) result = result * base % mod;
            index >>= 1;
            base = base * base % mod;
        }
        return result;
    }
    int main() {
        cin >> n >> m;
        for (int_t i = 1; i <= n; i++) {
            static char buff[1 << 10];
            scanf("%s", buff + 1);
            for (int_t j = 1; j <= m; j++) {
                prefix[i][j] =
                    (prefix[i][j - 1] + power(SEED, j - 1) * buff[j] % mod) % mod;
                pow_inv[j] = power(SEED, -j);
                pown2[j] = pown2[j - 1] / 2;
            }
        }
        for (int_t i = 1; i <= n; i++) {
            for (int_t j = 1; j <= n; j++) {
                for (int_t k = 1; k <= m; k++) {
                    if (prefix[i][k] ==
                        (2 * mod + prefix[j][m] - prefix[j][m - k]) % mod *
                            pow_inv[m - k] % mod) {
                        mat[i][j] += pown2[m - k];
                    }
                }
            }
            mat[i][n + 1] = -pown2[m];
            mat[n + 1][i] = 1;
        }
        mat[n + 1][n + 2] = 1;
        int_t len = n + 1;
        for (int_t i = 1; i <= len; i++) {
            for (int_t j = i + 1; j <= len; j++) {
                if (mat[i][i] == 0) mat[i][i] = EPS;
                real_t f = mat[j][i] / mat[i][i];
                for (int_t k = i; k <= len + 1; k++) {
                    mat[j][k] -= f * mat[i][k];
                }
            }
        }
        for (int_t i = len; i >= 1; i--) {
            if (mat[i][i] == 0 && mat[i][len + 1] == 0) continue;
            mat[i][len + 1] /= mat[i][i];
            mat[i][i] = 1;
            for (int_t j = i - 1; j >= 1; j--) {
                real_t f = mat[j][i];
                mat[j][i] = 0;
                mat[j][len + 1] -= f * mat[i][len + 1];
            }
        }
        std::cout.setf(std::ios::fixed);
        cout << std::setprecision(10);
        for (int_t i = 1; i <= n; i++) {
            cout << mat[i][len + 1] << endl;
        }
        return 0;
    }

     

  • CF 113D Museum

    $$ \text{设}\left( i,j \right) \text{表示第一个人在}i\text{号点,第二个人在}j\text{号点} \\ \text{首先一开始需要钦定一个两个人的汇合点,设这个汇合点为}\left( s,s \right) \\ \text{那么我们可以设计状态}P\left( i,j \right) \text{,表示从}\left( i,j \right) \text{出发,到达汇合点}\left( s,s \right) \text{的概率} \\ \text{显然,}P\left( s,s \right) =\text{1,}P\left( a,a \right) =0\left( a\ne s \right) \\ \text{对于其他的状态,可以构建方程} \\ P\left( i,j \right) =\sum_{\text{可以从}\left( i,j \right) \text{出发转移到的状态}\left( x,y \right)}{P\left( x,y \right) \times A_{\left( i,j \right) ,\left( x,y \right)}} \\ \text{其中}A_{\left( i,j \right) ,\left( x,y \right)}\text{表示从状态}\left( i,j \right) \text{转移到状态}\left( x,y \right) \text{的概率} \\ \text{设}f\left( i \right) \text{表示}i\text{号点的出边数量} \\ \text{分情况讨论一下} \\ i=x\text{且}j=y\text{时,两个人都没有动} \\ A_{\left( i,j \right) ,\left( x,y \right)}=p_i\times p_j \\ i=x\text{且}j\ne y,\text{第一个人没有走,第二个人随机走到了相邻的一个点} \\ A_{\left( i,j \right) ,\left( x,y \right)}=p_i\times \left( 1-p_j \right) \times \frac{1}{f\left( j \right)} \\ i\ne x\text{且}j=y\text{同理} \\ i\ne x\text{且}j\ne y,\text{这个时候}i\text{和}j\text{都走了} \\ A_{\left( i,j \right) ,\left( x,y \right)}=\left( 1-p_i \right) \times \frac{1}{f\left( i \right)}\times \left( 1-p_j \right) \times \frac{1}{f\left( j \right)} \\ \text{所以我们一共有}n^2-n\text{个方程组,可以高斯消元了} \\ \text{然而这种做法我们需要依次枚举汇合点在哪个点,所以时间复杂度是}O\left( n^7 \right) \text{的,会炸} \\ \text{但考虑到,每次枚举一个汇合点,我们改变的只有}P\left( i,i \right) \text{的取值。} \\ \text{所以我们可以把所有的}P\left( i,i \right) \text{视为常数项。} \\ \text{即解方程组} \\ Ax=B \\ \text{满足}A\text{是系数矩阵,}x\text{是未知数列向量,}B\text{是}n^2-n\text{个关于}P\left( i,i \right) \text{的行向量}\left( 1\le i\le n \right) \\ \text{执行完高斯消元后,我们可以得到}P\left( a,b \right) =\sum_{1\le i\le n}{A_i\times P\left( i,i \right)} \\ \text{然后依次令}P\left( i,i \right) \text{为1\ 求解即可} \\ \text{注意要特判}a=b\text{的情况} \\ \\ $$

    #include <iostream>
    #include <algorithm>
    #include <cstdio>
    #include <vector>
    #include <inttypes.h>
    #include <iomanip>
    #define debug(x) std::cout << #x << " = " << x << std::endl;
    
    typedef int int_t;
    typedef long double real_t;
    const int_t LARGE = 22;
    
    const real_t EPS = 1e-6;
    
    using std::cin;
    using std::endl;
    using std::cout;
    
    real_t mat[501][501];
    
    bool graph[23][23];
    int_t degree[LARGE + 1];
    real_t prob[LARGE + 1];
    int_t n,m,a,b;
    template<class T>
    T xabs(T x) {
    	if(x<0) return -x;
    	return x;
    }
    int_t encode(int_t a,int_t b) {
    //	cout<<"encode "<<a<<" "<<b<<" to "<< n * (a - 1) + b - ((a>b)?(a-1):a)<<endl;
    	return n * (a - 1) + b - ((a>b)?(a-1):a);
    }
    //´Ó״̬(i,j)×ªÒÆµ½(x,y)µÄ´ú¼Û
    real_t cost(int_t i,int_t j ,int_t x,int_t y) {
    	real_t result = 0;
    	if(i == x && j == y) result = prob[i] * prob[j];
    	if(i == x && j != y) {
    		result = prob[i] * (1 - prob[j]) * (1.0 / degree[j]) * graph[j][y];
    	}
    	if(i != x && j == y) result = (1 - prob[i]) * (1.0 / degree[i]) * graph[i][x] * prob[j];
    	if(i != x && j != y) result = (1 - prob[i]) * (1.0 / degree[i]) * graph[i][x] * (1 - prob[j]) * (1.0 / degree[j]) * graph[j][y];
    	return result;
    }
    int main() {
    	cin >> n >> m >> a >> b;
    	if(a == b) {
    		for(int_t i=1; i<=n; i++) {
    			if(i==a) cout<<1<<" ";
    			else cout<<0<<" ";
    		}
    		return 0;
    	}
    	for(int_t i = 1; i <= m; i++) {
    		int_t from,to;
    		cin >> from >> to;
    		graph[from][to] = graph[to][from] = 1;
    		degree[from] += 1;
    		degree[to] += 1;
    	}
    	for(int_t i = 1 ; i <= n; i++) {
    		cin >> prob[i];
    	}
    	int_t len = n * n - n;
    	//¹¹ÔìϵÊý¾ØÕó
    	for(int_t i = 1 ; i <= n; i++) {
    		for(int_t j = 1; j<= n ; j++) {
    			if(i == j) continue;
    			int_t curr = encode(i,j);
    			mat[curr][curr] = -1;
    			//ö¾Ù³ö±ß
    			for(int_t x = 1 ; x <= n; x++) {
    				for(int_t y = 1; y <= n ; y++) {
    					if(x == y) {
    						mat[curr][len + x] -= cost(i,j,x,y);
    					} else {
    						mat[curr][encode(x,y)] += cost(i,j,x,y);
    					}
    				}
    			}
    		}
    	}
    	//ö¾ÙÖ÷Ôª
    	for(int_t i = 1; i <= len; i++) {
    		if(xabs(mat[i][i]) < EPS) {
    			int_t pos = -1;
    			for(int_t j = i + 1; j <= len ; j++) {
    				if(xabs(mat[j][i]) >= EPS) {
    					pos = j;
    					break;
    				}
    			}
    			if(pos != -1) {
    				std::swap(mat[i],mat[pos]);
    			} else {
    				cout <<" exm ?? "<< endl;
    				return 0;
    			}
    		}
    		for(int_t j = i + 1; j <= len ; j++) {
    			real_t f = mat[j][i] / mat[i][i];
    			for(int_t k = i ; k <= n * n; k++) {
    				mat[j][k] -= f* mat[i][k];
    			}
    		}
    	}
    	//»Ø´ú
    	for(int_t i = len ; i >=1 ; i--) {
    		real_t f = mat[i][i];
    		for(int_t j = len + 1; j <= n * n ; j++) {
    			mat[i][j] /= f;
    		}
    		mat[i][i] /= f;
    		for(int_t j = i - 1; j >= 1; j--) {
    			real_t f2 =  mat[j][i];
    			mat[j][i] = 0;
    			for(int_t k = len + 1; k<= n * n; k++) {
    				mat[j][k] -= f2 * mat[i][k];
    			}
    		}
    	}
    //	print();
    	for(int_t i = 1; i <= n; i++) {
    		cout.setf(std::ios::fixed);
    		cout << std::setprecision(10) << mat[encode(a,b)][len + i] << " ";
    	}
    	return 0;
    }