In mathematics you don't understand things, you just get used to them.

标签 probablities 下的文章

Desc.

Link.

有 $n$ 个二元组, 用 $(a_i, t_i)$ 描述. 等概率在 $[1, k] \bigcap \mathbb{N}^*$ 中选取一个 $x$, 执行以下操作 $x$ 次:

  • 取一空数组 $\{b_n\}$, 令 $\displaystyle b_i = \sum_{t_j = i} a_j$, 然后再用 $b_i$ 替换每个 $a_i$;

问最终 $a_i$ 的期望值对 $\mathbf{998244353}$ 取模的结果;

$1\leqslant n \leqslant 2 \times 10^5$, $1 \leqslant k \leqslant 10^{18}$, $k$ is not a multiple of $\mathbf{998244353}$.

Sol.

令 $\mathbf A = \left( \begin{matrix} a_1, \dots, a_n \end{matrix}\right)$, $\mathbf T = \left (\begin{matrix} 0/1 & \cdots & 0/1 \\ \vdots & \ddots & \vdots \\ 0/1 & \cdots & 0/1 \end{matrix}\right)$, 其中 $\mathbf T_{i, t_i} = 1$, 其余为 $0$. 则答案为 (以下均省去 $\frac 1k$ 不写):

$$ f(\mathbf A, \mathbf T, k) = \sum_{i=1}^k \mathbf A \times \mathbf T^i $$

进行分治处理 (以下忽略 $\frac 1k$ 的系数).

$$ \begin{aligned} &f(\mathbf A, \mathbf T, 2m)\\ ={} & \mathbf A\times \mathbf T + (\mathbf A \times \mathbf T^2) + \dots + (\mathbf A \times \mathbf T ^{2m}) \\ ={} & f(\mathbf A\mathbf T+\mathbf A\mathbf T^2, \mathbf T^2, m) \end{aligned} $$

对于 $k$ 为奇数的情况:

$$ \begin{aligned} &f(\mathbf A, \mathbf T, k)\\ ={}& \mathbf A+f(\mathbf A \mathbf T, \mathbf T, k-1) \end{aligned} $$

于是我们用 $\mathcal O(n\log k)$ 的时间解决了问题. 当然亦可以建出内向基环树.

int main()
{
    ios::sync_with_stdio(0);
    cin.tie(nullptr);
    int n; ll k; cin >> n >> k;
    const int MOD = 998244353;
    using mint = modint<MOD>;
    using vm = vector<mint>;
    vi tmp(n), T(n); rds(T), rds(tmp);
    vm A(n);
    transform(allu(T), T.begin(), [&](int x) { return x-1; });
    transform(allu(tmp), A.begin(), [&](ll x) { return x%MOD; });
    const auto add = [&](const vm& lhs, const vm& rhs) {
        vm res(n);
        for (int i=0;i<n;++i) res[i] = lhs[i]+rhs[i];
        return res;
    };
    const auto mul = [&](const vm& lhs, const vi& to) {
        vm res(n);
        for (int i=0;i<n;++i) res[to[i]] += lhs[i];
        return res;
    };
    const auto trans = [&](const vi& v) {
        vi res(n);
        for (int i=0;i<n;++i) res[i] = v[v[i]];
        return res;
    };
    vm res(n);
    ll tmp_k = k%MOD;
    A = mul(A, T);
    while (k) {
        if (k&1) res = add(res, A), A = mul(A, T);
        A = add(A, mul(A, T));
        T = trans(T);
        k /= 2;
    }
    mint inv = 1/mint(tmp_k);
    for (int i=0;i<n;++i) cout << res[i]*inv << " \n"[i == n-1];
}

/ 相处的时间 你我已命运相连 /

/ 音乐频率凋谢 默契超越一切 /

/ 已记不清有多少个深夜你独自伏案桌前 /

/ 彻夜无眠 只为让我更耀眼 /

—— Kide - 星愿 ft. 星尘


date: '2023-10-12'
title: 'Solution -「ABC 313F」Flip Machines'


Desc.

Link.

有 $N$ 张卡片,第 $i$ 张卡片正面印着一个数 $A_i$,反面印着一个数 $B_i$。一开始所有数正面朝上。

有 $M$ 种操作,第 $i$ 种操作表示为:

  • 有 $50\%$ 的概率将卡片 $X_i$ 翻转,否则将 $Y_i$ 翻转。

要求你求一个集合 $S\subseteq \mathbb{N} \bigcap [1,m]$,使得进行了集合中所有的编号的操作之后正面朝上的所有数的和的期望最大。输出这个最大值。

Sol.

容易发现有概率被任意正数台机器翻转的卡片对答案的贡献是 $\frac{a_i+b_i}{2}$, 其他卡片的贡献是 $a_i$. 将卡片分为两个集合 $A$ 和 $B$, 其中 $\forall i \in A$ 有 $a_i > b_i$, $\forall i \in B$ 有 $a_i \leqslant b_i$. 把一台机器的操作 $(x_i, y_i)$ 分为三个集合 $P, Q$ 和 $X$, 其中 $P$ 中的边只涉及 $A$ 内的卡片, $Q$ 中的边只涉及 $B$ 中的卡片, $X$ 中的边同时涉及两边的卡片. $P$ 中的边一定不会连, $Q$ 中的边一定会连, 留给我们决策的是 $X$ 中的边.

分类讨论 $|A|$ 和 $|B|$ 的相对大小:

  • $|A| > |B|$:

此时 $|B|<\frac n2$, 采取以下的动态规划, 设 $f_{i, S}$ 表示 $A$ 中的前 $i$ 个点与 $B$ 的点集 $S$ 有连边的最大期望答案. 有转移

$$ f_{i, S}\rightarrow f_{i+1, S\cup\{j\}} $$

  • $|A| \leqslant |B|$:

朴素的应用搜索即可.

代码有点脑溢血, 抄了官方题解的实现.

int main()
{
    ios::sync_with_stdio(0);
    cin.tie(nullptr);
    int n, m; cin >> n >> m;
    vi a(n), b(n), u(m), v(m);
    for (int i=0;i<n;++i) cin >> a[i] >> b[i];
    for (int i=0;i<m;++i) cin >> u[i] >> v[i];
    for (int i=0;i<m;++i) u[i]--, v[i]--;
    for (int i=0;i<m;++i) {
        if (u[i] == v[i] && a[u[i]] < b[u[i]]) swap(a[u[i]], b[u[i]]);
    }
    int ans = 0;
    vi A, B, id(n);
    for (int i=0;i<n;++i) {
        ans += a[i]*2;
        if (a[i] >= b[i]) id[i] = A.size(), A.pb(a[i]-b[i]);
        else id[i] = B.size(), B.pb(b[i]-a[i]);
    }
    int nA = A.size(), nB = B.size();
    vector<ll> ls(nA);
    const int INF = 1e9;
    for (int i=0;i<m;++i) {
        bool xp = a[u[i]] >= b[u[i]],
             yp = a[v[i]] >= b[v[i]];
        if (xp == yp) {
            if (!xp) {
                ans += B[id[u[i]]]+B[id[v[i]]];
                B[id[u[i]]] = B[id[v[i]]] = 0;
            }
               
        } else {
            if (!xp) swap(u[i], v[i]);
            ls[id[u[i]]] |= 1ll<<id[v[i]];
        }
    }
    if (nA <= nB) {
        int opt = 0;
        for (int i=0;i<1<<nA;++i) {
            int cur = 0; ll st = 0;
            for (int j=0;j<nA;++j) if (i&(1<<j)) cur -= A[j], st |= ls[j];
            for (int j=0;j<nB;++j) if (st&(1ll<<j)) cur += B[j];
            chkmax(opt, cur);
        }
        ans += opt;
    } else {
        vector dp(nA+1, vi(1<<nB, -INF));
        dp[0][0] = 0;
        for (int i=0;i<nA;++i) {
            for (int j=0;j<1<<nB;++j) {
                chkmax(dp[i+1][j], dp[i][j]);
                chkmax(dp[i+1][j|ls[i]], dp[i][j]-A[i]);
            }
        }
        int mx = 0;
        for (int j=0;j<1<<nB;++j) {
            int cur = dp[nA][j];
            for (int i=0;i<nB;++i) if (j&(1<<i)) cur += B[i];
            chkmax(mx, cur);
        }
        ans += mx;
    }
    cout << fixed << ans / 2. << "\n";
}

/ 当最后一缕余温消亡 /

/ 空躯壳何处能够安放 /

/ 睁眼醒于温暖的慌 /

/ 一笑而过继续坚强 /

—— 孤寂code - 逆光 ft. 洛天依

link.

对于 pass 1, 你把他考虑成 $\frac{\sum x}{i}$ 的形式, 于是每次操作的贡献就是 $\frac{2}{i}$, 那么答案就是 $\sum_{i=2}^n\frac{2}{i}$.

对于 pass 2, 首先由全概率公式 $\textbf{E}(n)=\sum_{i=0}^{\infty}\textbf{P}(n\geqslant i)$. 那么我们考虑固定 $k=i$, 求出 $\textbf P(n\geqslant k)$. 其实直接求这个有点麻烦, 强化一下条件, 我们求 $\textbf P(n=k)$, 最后后缀和即可取得 $\textbf P(n\geqslant k)$.

设 $f(i,j)$ 表示共 $i$ 个叶结点, 树深度为 $j$ 时的概率. 对于二叉树我们有一种经典的考虑子结构的方式, 即考虑左子树的信息. 对于此题, 我们考虑左子树 (相对于根结点来说) 的叶结点数量 $k$. 那么有 $f(i,\max\{k,l\}+1)\gets f(i,\max\{k,l\}+1)+\frac{f(j,k)\times f(i-j,l)}{i-1}$. 其中 $i$ 为全树的叶子个数, $k$, $l$ 分别为左, 右子树的深度, $j$ 为左子树的叶结点个数. $\frac{1}{i-1}$ 是对于一个有 $i$ 个叶子的二叉树, 有 $k$ 个结点在左子树的概率, which 与 $k$ 无关. 至于为什么是 $\frac{1}{i-1}$, 你考虑操作次数为 $i-1$ 即可.

#include <bits/stdc++.h>
using namespace std;
typedef double db;
int Q, N;
db f[200][200];
signed main() {
  ios::sync_with_stdio(0), cin.tie(0);
  cin >> Q >> N;
  cout << fixed;
  cout << setprecision(6);
  if (Q == 1) {
    db ans = 0;
    for (int i = 2; i <= N; ++i) {
      ans += 2.0 / i;
    }
    cout << ans << "\n";
  } else if (Q == 2) {
    f[1][0] = f[2][1] = f[3][2] = 1;
    for (int i = 4; i <= N; ++i) {
      for (int j = 1; j < i; ++j) {
        for (int k = 0; k < j; ++k) {
          for (int l = 0; l < i - j; ++l) {
            f[i][max(k, l) + 1] += f[j][k] * f[i - j][l] / (i - 1);
          }
        }
      }
    }
    db ans = 0;
    for (int i = N; i >= 1; --i) {
      f[N][i] += f[N][i + 1];
    }
    for (int i = 1; i < N; ++i) {
      ans += f[N][i];
    }
    cout << ans << "\n";
  }
  return 0;
}

Desciption

Link.

给定一个值域在 $[1,x]$ 的长度为 $n$ 的序列(由随机数构成),求给定一组区间中的最小值的最大值的期望。

Solution

记:

$$ w=\max\{\min\{a_{l_{j}},a_{l_{j}+1},\cdots,a_{r_{j}}\}|j\in[1,q]\} $$

因为我们最后取的是 $\max$,不能直接用全概率公式,转化一下:

$$ E(w)=\sum_{i=0}^{\infty}P(w\ge i)=\sum_{i=0}^{\infty}1-P(w<i) $$

这意味着每一个被询问区间中的最小值都需 $<i$。也就是说,每一个区间至少需要一个 $<i$ 的数。

这对于每一个区间来说概率为 $\frac{i-1}{x}$。又因为区间可能出现相交,所以我们考虑用点去被包含于区间。

当然,一个区间包含另一个区间,这个区间肯定是没有用的。然后把区间按左右端点分别为第一、第二关键字排序。

枚举 $w$,设 $f_{i}$ 表示区间右端点在 $i$ 之前的所有区间满足条件的概率。

$$ f_{i}=\frac{w-1}{x}\times\sum_{j=0}^{i}f_{j}\times(1-\frac{w-1}{x})^{i-j-1} $$

#include <cstdio>

using i64 = long long;

const int MOD = 666623333;
const int MAXN = 2e3 + 5;

int n, x, q, ar[MAXN];
i64 f[MAXN][2], ff[MAXN][2];

void imax ( int& a, const int b ) { a = a < b ? b : a; }
int add ( const int a, const int b, const int p = MOD ) { return a + b < p ? a + b : ( a + b ) % p; }
int sub ( const int a, const int b, const int p = MOD ) { return a - b < 0 ? a - b + p : a - b; }
int mul ( const i64 a, const i64 b, const int p = MOD ) { return a * b % p; }
int cpow ( int bas, int idx = MOD - 2 ) {
    int res = 1;
    while ( idx ) {
        if ( idx & 1 )    res = mul ( res, bas );
        bas = mul ( bas, bas ), idx >>= 1;
    }
    return res % MOD;
}

int main () {
    scanf ( "%d%d%d", &n, &x, &q );
    for ( int i = 1, tmpl, tmpr; i <= q; ++ i )    scanf ( "%d%d", &tmpl, &tmpr ), imax ( ar[tmpr + 1], tmpl );
    for ( int i = 1; i <= n + 1; ++ i )    imax ( ar[i], ar[i - 1] );
    i64 ix = cpow ( x ), ans = 0;
    for ( int i = 1; i <= x; ++ i ) {
        i64 p = mul ( i - 1, ix ) % MOD, ip = cpow ( 1 - p ), s;
        ff[0][0] = ff[0][1] = 1;
        for ( int j = 1; j <= n; ++ j )    ff[j][0] = mul ( ff[j - 1][0], 1 - p ) % MOD, ff[j][1] = mul ( ff[j - 1][1], ip ) % MOD;
        f[0][0] = 0, f[0][1] = 1;
        for ( int j = 1; j <= n; ++ j ) {
            f[j][0] = mul ( mul ( p, sub ( f[j - 1][1], ar[j] ? f[ar[j] - 1][1] : 0 ) ) % MOD, ff[j - 1][0] ) % MOD;
            f[j][1] = add ( mul ( f[j][0], ff[j][1] ) % MOD, f[j - 1][1] ) % MOD;
        }
        s = 0;
        for ( int j = ar[n + 1]; j <= n; ++ j )    s = add ( s, mul ( f[j][0], ff[n - j][0] ) % MOD ) % MOD;
        ans = sub ( add ( ans, 1 ) % MOD, s );
    }
    printf ( "%lld\n", ans % MOD );
    return 0;
}

Description

Link.

这道题 的加强版。

Solution

题解里面大多数都是概率 DP,或者是期望 DP 然后是逆推。甚至不给 DP 的转移式。机房 yyds Reanap 发了一篇逆推的题解,那我就来补一篇正推的期望 DP 的填表法做法。

首先这道题看上去好像可以状压的样子,我们可以设 $f_{i,S}$ 表示当前打了 $i$ 次,敌方的情况是 $S$ 的期望。

不过仔细想一下发现我们只需要知道各种血量的奴隶主有多少即可。

于是我们重新设计 DP 的状态:$f_{s,i,j,k}$ 表示目前打了 $s$ 次,敌方分别有 $i$、$j$、$k$ 个 1hp、2hp、3hp 的奴隶主。

首先我们令 $T=[i+j+k<K]$

那么我们的方程就是:

$$ f_{s,i,j,k}=\begin{cases} f_{s-1,i-1,j,k}+\frac{i}{i+j+k+1},M=1\land i\neq0 \\ f_{s-1,i-1,j,k}+\frac{i}{i+j+k+1},M=2\land i\neq0 \\ f_{s-1,i+1,j-1+T,k}+\frac{j}{i+j+k+1},M=2\land j\neq0 \\ f_{s-1,i-1,j,k}+\frac{i}{i+j+k+1},M=3\land i\neq0 \\ f_{s-1,i+1,j-1,k+T}+\frac{j}{i+j+k+1},M=3\land j\neq0 \\ f_{s-1,i,j+1,k-1+T}+\frac{k}{i+j+k+1},M=3\land k\neq0 \end{cases} $$

这个方程挺好理解的,基本就等于照题意模拟。

然后我们发现转移式中的系数部分和 $f$ 数组没有关系,所以我们可以用矩阵来加速这个东西。

数一数状态数,直接加速直接 T 飞。

有一个矩阵加速常用的 trick,预处理矩阵 2 的幂。

然后取模卡卡常即可。

(代码不保证稳定能过)

#include <cstdio>
#define mod ( 998244353 )

using namespace std;
typedef long long LL;

char buf[1 << 21], *p1 = buf, *p2 = buf;
#define getchar( ) ( p1 == p2 && ( p2 = ( p1 = buf ) + fread( buf, 1, 1 << 21, stdin ), p1 == p2 ) ? EOF : *p1 ++ )

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

template<typename _T, typename... Args>
void read( _T &t, Args&... args ) { read( t ), read( args... ); }

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

template<typename _T>
void Add( _T &x, _T y )
{
    if( y >= mod )  y %= mod;
    x += y;
    if( x >= mod )  x -= mod;
}

template<typename _T>
_T square( _T x ) { return x * x; }

const int Maxn = 10 + 5, Maxk = 170 + 5;
int T, M, K, S, Unite[Maxn][Maxn][Maxn];
LL tmp[Maxk], Ans[Maxk], Inv[Maxn];
struct Matrix
{
    LL mat[Maxk][Maxk];

    friend Matrix operator * ( const Matrix &one, const Matrix &another )
    {
        Matrix res;
        for( int i = 0; i <= S + 1; ++ i )
        {
            for( int j = 0; j <= S + 1; ++ j )
            {
                res.mat[i][j] = 0;
                for( int k = 0; k <= S + 1; ++ k )    Add( res.mat[i][j], one.mat[i][k] * another.mat[k][j] );
            }
        }
        return res;
    }
} dp[Maxk];

template<typename _T>
_T qkpow( _T base, _T times )
{
    _T res = 1;
    while( times )
    {
        if( times & 1 )    res = ( LL )res * base % mod;
        base = ( LL )base * base % mod;
        times >>= 1;
    }
    return res;
}

void progressInversions( ) { for( int i = 0; i <= 10; ++ i )    Inv[i] = qkpow( i, mod - 2 ); }
signed main( )
{
    progressInversions( );
    read( T, M, K );
    for( int i = 0; i <= K; ++ i )
    {
        int UpI;
        if( M > 1 )    UpI = K - i;
        else    UpI = 0;
        for( int j = 0; j <= UpI; ++ j )
        {
            int UpJ;
            if( M > 2 )    UpJ = K - i - j;
            else    UpJ = 0;
            for( int k = 0; k <= UpJ; ++ k )    Unite[i][j][k] = ++ S;
        }
    }
    for( int i = 0; i <= K; ++ i )
    {
        int UpI;
        if( M > 1 )    UpI = K - i;
        else    UpI = 0;
        for( int j = 0; j <= UpI; ++ j )
        {
            int UpJ;
            if( M > 2 )    UpJ = K - i - j;
            else    UpJ = 0;
            for( int k = 0; k <= UpJ; ++ k )
            {
                int Add;
                if( i + j + k < K )    Add = 1;
                else    Add = 0;
                if( M == 1 && i )    dp[0].mat[Unite[i][j][k]][Unite[i - 1][j][k]] = ( LL )i * Inv[i + j + k + 1] % mod;
                else if( M == 2 )
                {
                    if( i )    dp[0].mat[Unite[i][j][k]][Unite[i - 1][j][k]] = ( LL )i * Inv[i + j + k + 1] % mod;
                    if( j )    dp[0].mat[Unite[i][j][k]][Unite[i + 1][j - 1 + Add][k]] = ( LL )j * Inv[i + j + k + 1] % mod;
                }
                else if( M == 3 )
                {
                    if( i )    dp[0].mat[Unite[i][j][k]][Unite[i - 1][j][k]] = ( LL )i * Inv[i + j + k + 1] % mod;
                    if( j )    dp[0].mat[Unite[i][j][k]][Unite[i + 1][j - 1][k + Add]] = ( LL )j * Inv[i + j + k + 1] % mod;
                    if( k )    dp[0].mat[Unite[i][j][k]][Unite[i][j + 1][k - 1 + Add]] = ( LL )k * Inv[i + j + k + 1] % mod;
                }
                dp[0].mat[Unite[i][j][k]][Unite[i][j][k]] = dp[0].mat[Unite[i][j][k]][S + 1] = Inv[i + j + k + 1];
            }
        }
    }
    dp[0].mat[S + 1][S + 1] = 1;
    for( int i = 1; i <= 60; ++ i )    dp[i] = square( dp[i - 1] );
    while( T -- > 0 )
    {
        LL N;
        read( N );
        for( int i = 0; i <= S + 1; ++ i )  Ans[i] = 0;
        if( M == 1 )    Ans[Unite[1][0][0]] = 1;
        else if( M == 2 )    Ans[Unite[0][1][0]] = 1;
        else    Ans[Unite[0][0][1]] = 1;
        for( int i = 0; i <= 60; ++ i )
        {
            if( ( N >> i ) & 1 )
            {
                for( int j = 0; j <= S + 1; ++ j )
                {
                    tmp[j] = 0;
                    for( int k = 0; k <= S + 1; ++ k )    Add( tmp[j], Ans[k] * dp[i].mat[k][j] );
                }
                for( int j = 0; j <= S + 1; ++ j )    Ans[j] = tmp[j];
            }
        }
        write( Ans[S + 1] ), putchar( '\n' );
    }
    return 0;
}