【LuoguP3600】随机数生成器-概率DP+双指针

测试地址：随机数生成器
做法：本题需要用到概率DP+双指针。
考虑离散概率情况下的期望公式：
E[ans]=∑xs=1s⋅P(ans=s)$E[ans]=\sum _{s=1}^{x}s\cdot P(ans=s)$
也就相当于：
E[ans]=∑xs=1P(ans≥s)$E[ans]=\sum _{s=1}^{x}P(ans\ge s)$
考虑到P(ans≥s)=1−P(ans≤s−1)$P(ans\ge s)=1-P(ans\le s-1)$，我们只需计算出所有P(ans≤s)$P(ans\le s)$即可。
所有最小值的最大值≤s$\le s$，就相当于在每个区间中至少存在一个≤s$\le s$的数。那么问题就变成，每个位置有sx$\frac{s}{x}$的概率被选，求每个区间都至少有一个数被选中的概率。
但是原本所给的区间比较杂乱，我们需要对这些限制条件做一些简化。显然，如果一个区间严格包含另一个区间，那么这个区间就没有用了，因为被包含那个区间中的数被选中也就意味着该区间的数被选中。所以我们把这样的区间去掉之后，剩下的区间按左端点从小到大排序，右端点也是不降的。
你有可能在此时想到了直接概率DP的办法，但我这里要介绍另一种办法。显然因为区间的单调性，每个位置被选中所能影响到的区间的区间（没有打错）也是单调的。这样问题就变成，每个区间有P$P$的概率被选，求选中的区间覆盖整个序列的概率。
那么我们就可以DP了。令f(i)$f(i)$为前i$i$个区间中，选择第i$i$个并且能覆盖从头到第i$i$个区间右端点这一段序列的概率。有状态转移方程：
f(i)=P⋅(∑rj≥li−1f(j)⋅(1−P)i−j−1+[li=1]⋅(1−P)i−1)$f(i)=P\cdot (\sum _{r_j\ge l_i-1}f(j)\cdot (1-P{)}^{i-j-1}+[l_i=1]\cdot (1-P{)}^{i-1})$
实际上就是枚举上一个选择的区间j$j$，显然必须要满足rj≥li−1$r_j\ge l_i-1$才可能覆盖整段。而既然要满足j$j$是“上一个”，那么中间的i−j−1$i-j-1$个就不能选，因此要乘上(1−P)i−j−1$(1-P{)}^{i-j-1}$。而后面要加上的那个东西，是因为当li=1$l_i=1$时，是有可能没有“上一个”选择的区间的，所以要多算上这种情况的概率。在这种情况下，最后的答案显然为∑ri=nf(i)⋅(1−P)t−i$\sum _{r_i=n}f(i)\cdot (1-P{)}^{t-i}$，其中n$n$表示处理后询问的总数，t$t$表示有用的位置总数（不能影响任何区间的位置就是没用的位置，对答案没有影响）。
这个方程是O(n2)$O(n^2)$的，不能从方程的角度优化了，因此我们考虑从转移的角度优化。注意到可以把(1−P)i−j−1$(1-P{)}^{i-j-1}$拆成(1−P)−j⋅(1−P)i−1$(1-P{)}^{-j}\cdot (1-P{)}^{i-1}$，那么我们就是要求出∑rj≥li−1f(j)⋅(1−P)−j$\sum _{r_j\ge l_i-1}f(j)\cdot (1-P{)}^{-j}$。而因为区间的单调性，显然可以用双指针的方法维护这个东西，那么转移就是O(n)$O(n)$的了，而对于1−P$1-P$的一些幂可以O(n)$O(n)$预处理。于是我们就以O(n2)$O(n^2)$的总时间复杂度完成了这一题。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=666623333;
int n,q;
ll x;
struct interval
{
    int l,r;
}s[2010];
bool vis[2010]={0};
int L[2010],R[2010];
ll f[2010],antiP[2010][2];

ll power(ll a,ll b)
{
    ll s=1,ss=a;
    while(b)
    {
        if (b&1) s=s*ss%mod;
        b>>=1;ss=ss*ss%mod;
    }
    return s;
}

bool cmp(interval a,interval b)
{
    return a.l<b.l;
}

int main()
{
    scanf("%d%lld%d",&n,&x,&q);
    for(int i=1;i<=q;i++)
        scanf("%d%d",&s[i].l,&s[i].r);
    for(int i=1;i<=q;i++)
        for(int j=1;j<=q;j++)
            if (i!=j&&s[i].l<=s[j].l&&s[i].r>=s[j].r)
                if (s[i].l<s[j].l||s[i].r>s[j].r)
                    vis[i]=1;
    int tot=0;
    for(int i=1;i<=q;i++)
        if (!vis[i])
        {
            s[++tot].l=s[i].l;
            s[tot].r=s[i].r;
        }
    q=tot;
    sort(s+1,s+q+1,cmp);

    int l=1,r=0;
    tot=0;
    for(int i=1;i<=n;i++)
    {
        while(r<q&&s[r+1].l<=i) r++;
        while(l<=q&&s[l].r<i) l++;
        if (l<=r) L[++tot]=l,R[tot]=r;
    }

    ll ans=0,lastansP=0;
    for(ll p=1;p<=x-1;p++)
    {
        ll ansP=0;
        ll P=p*power(x,mod-2ll)%mod;
        antiP[0][0]=antiP[0][1]=1;
        antiP[1][0]=(1-P+mod)%mod;
        antiP[1][1]=power((1-P+mod)%mod,mod-2);
        for(int i=2;i<=tot;i++)
        {
            antiP[i][0]=(antiP[i-1][0]*antiP[1][0])%mod;
            antiP[i][1]=(antiP[i-1][1]*antiP[1][1])%mod;
        }

        l=0;
        ll nowsum=0;
        for(int i=1;i<=tot;i++)
        {
            if (i>1) nowsum=(nowsum+f[i-1]*antiP[i-1][1])%mod;
            while((l==0&&i>1)||R[l]<L[i]-1)
            {
                l++;
                if (l>1) nowsum=((nowsum-f[l-1]*antiP[l-1][1])%mod+mod)%mod;
            }
            f[i]=nowsum*antiP[i-1][0]%mod;
            if (L[i]==1) f[i]=(f[i]+antiP[i-1][0])%mod;
            f[i]=f[i]*P%mod;
            if (R[i]==q) ansP=(ansP+f[i]*antiP[tot-i][0])%mod;
        }

        ans=((ans+p*(ansP-lastansP))%mod+mod)%mod;
        lastansP=ansP;
    }
    printf("%lld",((ans+x*(1ll-lastansP))%mod+mod)%mod);

    return 0;
}