[2017集训队作业自选题#149]小c的岛屿

最新推荐文章于 2020-10-30 21:25:47 发布

WerKeyTom_FTD

最新推荐文章于 2020-10-30 21:25:47 发布

阅读量1.4k

点赞数 1

分类专栏：高斯消元概率与期望一般动规与递推容斥原理

本文链接：https://blog.csdn.net/WerKeyTom_FTD/article/details/78947101

版权

一般动规与递推同时被 3 个专栏收录

120 篇文章 0 订阅

订阅专栏

容斥原理

41 篇文章 0 订阅

订阅专栏

概率与期望

31 篇文章 0 订阅

订阅专栏

前言

感觉这题挺棒棒的。

题目描述

小c有n个岛屿。她为了管理岛上的人，决定让这些岛屿连通。
每个岛屿i有一个概率值pi，和一个友好列表Ai。
小c首先来到了1号岛屿，她会依次执行下面操作：
1、设现在在岛屿x，有px的概率产生一条无向边连接两个随机岛屿，这两个岛屿不会相同，也不会已经有边相连。（即在尚不存在的无向边中随机一条加入图中，不会加自环）
2、如果此时所有岛屿连通，她就会心满意足地离开。
3、否则她会在当前点的友好列表中，随机选择一个，走到那座岛屿上。她的不满意度会+1，并重复第1步。
求小c的期望不满意度。
题目中的概率都是在模10^9+7域下给出的，答案也需要在模域下求。
即，设答案是a/b，b的逆元是x，你需要输出ax mod 10^9+7的值。

n^6做法

我们设e(x,i)表示当前在x号点，已经连了i条边，还要期望经过多少步才能使图联通。
我们设q表示i条边不连通的图在加入了一条不重复的边之后联通的概率。
设dx表示点x的度数。
那么 $e(x,i)=p_x(1-q)\frac{1}{d_x}\sum_{x->y}(e(y,i+1)+1)+(1-p_x)\frac{1}{d_x}\sum_{x->y}(e(y,i)+1)$
我们分层来做，可以注意到如果q能够求出，每层内可以高斯消元。
那么复杂度是m*n^3=n^5。
现在问题是这个q怎么求。
显然 $q=\frac{所有i+1条边联通图的桥边个数和}{(n*(n-1)/2-i)*i条边不连通图的个数}$
我们来考虑怎么求得分母和分子。

任意条边联通图的桥边个数和

假设要求m条边的，现在一条桥边把图分成两个联通块，我们不妨枚举1号点所在联通块有k个点，l条边，那么贡献：
$C_{n-1}^{k-1}*k*(n-k)*(k个点l条边联通块数)*(n-k个点m-l-1条边联通块数)$
其中组合数表示除了1号点其余点的编号选取，k和n-k则表示两边各要选出一个编号作为桥边的端点。
求一个m需要nm的复杂度，要求任意条边，这里复杂度是m*nm=n^5。
现在要解决括号里的问题。

任意个点任意条边联通图个数

设f(n,m)表示n个点m条边的联通图个数。
可以枚举1号点所在联通块的点数和边数，乘个组合数即可。
复杂度是nmnm=n^6。
因此总复杂度是O(n^6)，显然不能过。

n^5做法

容易发现瓶颈，我们现在用斯特林数容斥来优化任意点任意边。
原理请参考容斥的原理及广义应用
我们设 $dp(i,j)=\sum_{G\ is\ a\ graph,v(G)=i,e(G)=j}(-1)^{c(G)-1}*(c(G)-1)!$
其中 $v(G)$ 表示图G的点数， $e(G)$ 表示图G每个联通块如果都是完全图，总边数是多少。 $c(G)$ 表示图G的联通块数。
假如算出了dp，那么f很好算。
$f(n,m)=\sum_{i>=m}dp(n,i)*C_i^m$
我们先考虑如何计算dp，可以设一个g。
$g(i,j,k)=\sum_{G\ is\ a\ graph,v(G)=i,e(G)=j,c(G)=k}(-1)^{k-1}*(k-1)!$
这个g就很好转移了，可以枚举1号点所在联通块的点数。
复杂度是n*m*n*n=n^5。
然后有了g很好得到dp。
考虑用dp来得到f，可以用n*m*m=n^5时间来得出。
至此，所有部分的复杂度均为n^5，可惜如果不是常数过分优秀，仍然无法通过。这段代码后将介绍复杂度更优秀的稳过算法。
本份代码只能对n<=40在时限内求解。

#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
typedef long long ll;
const int maxn=50+10,mo=1000000007;
int dp[maxn][maxn*maxn],f[maxn][maxn*maxn],ff[maxn][maxn*maxn],fff[maxn][maxn*maxn],g[maxn*maxn],c[maxn*maxn][maxn*maxn],a[maxn][maxn],b[maxn],d[maxn],p[maxn];
int e[maxn*maxn][maxn],num[maxn];
int h[maxn],go[maxn*maxn],nxt[maxn*maxn],fac[maxn*maxn],inv[maxn*maxn];
int i,j,k,l,r,t,n,m,ans,tot,top,q,qq,x,y,ni,mx;
bool czy;
int qsm(int x,int y){
    if (!y) return 1;
    int t=qsm(x,y/2);
    t=(ll)t*t%mo;
    if (y%2) t=(ll)t*x%mo;
    return t;
}
void add(int x,int y){
    d[x]++;
    go[++tot]=y;
    nxt[tot]=h[x];
    h[x]=tot;
}
int S(int n){
    return n*(n-1)/2;
}
int C(int n,int m){
    if (n<m||m<0) return 0;
    return (ll)fac[n]*inv[m]%mo*inv[n-m]%mo;
}
void gauss(){
    int i,j,k,r,t;
    fo(i,1,n){
        fo(j,i,n)
            if (a[j][i]) break;
        fo(k,i,n) swap(a[i][k],a[j][k]);
        if (!a[i][i]) czy=1;
        fo(j,i+1,n){
            if (!a[j][i]) continue;
            t=(ll)a[j][i]*qsm(a[i][i],mo-2)%mo;
            fo(k,i,n){
                r=(ll)a[i][k]*t%mo;
                //(a[j][k]=(ll)a[i][k]*t%mo-a[j][k])%=mo;
                a[j][k]=r-a[j][k]<0?r-a[j][k]+mo:r-a[j][k];
            }
            b[j]=((ll)b[i]*t%mo-b[j])%mo;
        }
        /*j=i+1;
        while (j+2<=n){
            t=(ll)a[j][i]*qsm(a[i][i],mo-2)%mo;
            fo(k,i,n) (a[j][k]=(ll)a[i][k]*t%mo-a[j][k])%=mo;
            b[j]=((ll)b[i]*t%mo-b[j])%mo;

            t=(ll)a[j+1][i]*qsm(a[i][i],mo-2)%mo;
            fo(k,i,n) (a[j+1][k]=(ll)a[i][k]*t%mo-a[j+1][k])%=mo;
            b[j+1]=((ll)b[i]*t%mo-b[j+1])%mo;

            t=(ll)a[j+2][i]*qsm(a[i][i],mo-2)%mo;
            fo(k,i,n) (a[j+2][k]=(ll)a[i][k]*t%mo-a[j+2][k])%=mo;
            b[j+2]=((ll)b[i]*t%mo-b[j+2])%mo;

            j+=3;
        }
        r=j;
        fo(j,r,n){
            t=(ll)a[j][i]*qsm(a[i][i],mo-2)%mo;
            fo(k,i,n) (a[j][k]=(ll)a[i][k]*t%mo-a[j][k])%=mo;
            b[j]=((ll)b[i]*t%mo-b[j])%mo;
        }*/
    }
    fd(i,n,1){
        fd(j,n,i+1) (b[i]-=(ll)a[i][j]*num[j]%mo)%=mo;
        num[i]=(ll)b[i]*qsm(a[i][i],mo-2)%mo;
        b[i]=0;
    }
}
int main(){
    freopen("island5.in","r",stdin);
    scanf("%d",&n);
    m=n*(n-1)/2;
    mx=max(n,m);
    fo(i,1,n) scanf("%d",&p[i]);
    fo(i,1,n){
        scanf("%d",&t);
        fo(j,1,t){
            scanf("%d",&k);
            add(i,k);
        }
    }
    fac[0]=1;
    fo(i,1,mx) fac[i]=(ll)fac[i-1]*i%mo;
    inv[mx]=qsm(fac[mx],mo-2);
    fd(i,mx-1,0) inv[i]=(ll)inv[i+1]*(i+1)%mo;
    fo(i,1,n) dp[i][S(i)]=1;
    fo(i,0,n)
        fo(j,0,m)
            if (dp[i][j])
                fo(k,1,n-i)
                    if (j+(l=S(k))<=m){
                        t=(ll)dp[i][j]*C(i+k-1,k)%mo;
                        //(dp[i+k][j+S(k)]-=(ll)dp[i][j]*C(i+k-1,k)%mo)%=mo;    
                        dp[i+k][j+l]=dp[i+k][j+l]-t<0?dp[i+k][j+l]-t+mo:dp[i+k][j+l]-t;
                    }
    fo(i,0,n)   
        fo(j,0,m)
            fo(k,j,m){
                t=(ll)dp[i][k]*C(k,j)%mo;
                //(f[i][j]+=(ll)dp[i][k]*C(k,j)%mo)%=mo;
                f[i][j]=f[i][j]+t>=mo?f[i][j]+t-mo:f[i][j]+t;
            }
    fo(i,0,n)
        fo(j,0,m){
            fff[i][j]=ff[i][j]=f[i][j];
            fff[i][j]=(ll)fff[i][j]*i%mo*inv[i]%mo;
            //if (i) ff[i][j]=(ll)ff[i][j]*i%mo*inv[i-1]%mo;
            ff[i][j]=(ll)fff[i][j]*i%mo;
        }
    //return 0;
    fo(j,n-1,m){
        fo(k,1,n-1){
            r=min(j-1,k+j-n);
            fo(l,k-1,r){
                t=(ll)fff[k][l]*ff[n-k][j-l-1]%mo;
                //(g[j]+=(ll)fff[k][l]*ff[n-k][j-l-1]%mo)%=mo;
                g[j]=g[j]+t>=mo?g[j]+t-mo:g[j]+t;
                //(g[j]+=(ll)C(n-1,k)*k%mo*(n-k)%mo*f[k][l]%mo*f[n-k][j-l-1]%mo)%=mo;
            }
        }
        g[j]=(ll)g[j]*fac[n-1]%mo;
    }
    fo(i,1,n) e[i][m]=0;
    //return 0;
    fd(i,m-1,0){
        fo(j,1,n)
            fo(k,j,n) a[j][k]=0;
        //fo(j,1,n) b[j]=0;
        q=(C(m,i)-f[n][i])%mo;
        q=(ll)q*(m-i)%mo;
        q=(ll)g[i+1]*qsm(q,mo-2)%mo;
        (q+=mo)%=mo;
        qq=1-q;
        (qq+=mo)%=mo;
        fo(x,1,n){
            ni=qsm(d[x],mo-2);
            /*(b[x]+=1-p[x])%=mo;
            (b[x]+=mo)%=mo;*/
            t=h[x];
            while (t){
                y=go[t];
                //r=(ll)p[x]*qq%mo*ni%mo*(e[y][i+1]+1)%mo;
                r=e[y][i+1]+1>=mo?e[y][i+1]+1-mo:e[y][i+1]+1;
                //(b[x]+=(ll)p[x]*(1-q)%mo*ni%mo*(e[y][i+1]+1)%mo)%=mo;
                b[x]=b[x]+r>=mo?b[x]+r-mo:b[x]+r;
                (a[x][y]-=(ll)(1-p[x])*ni%mo)%=mo;
                t=nxt[t];
            }
            b[x]=(ll)b[x]*p[x]%mo*qq%mo*ni%mo;
            (b[x]+=1-p[x])%=mo;
            (b[x]+=mo)%=mo;
            (a[x][x]+=1)%=mo;
        }
        fo(j,1,n){
            fo(k,1,n) (a[j][k]+=mo)%=mo;
            (b[j]+=mo)%=mo;
        }
        gauss();
        fo(x,1,n){
            e[x][i]=num[x];
            (e[x][i]+=mo)%=mo;
        }
    }
    ans=e[1][0];
    (ans+=mo)%=mo;
    printf("%d\n",ans);
}

n^4做法

我们可能需要把原来的思路推倒一部分，才能做到如此优秀的复杂度。
现在我们改dp状态，设e(x,i)表示现在在x号点，还需要i条边，期望步数是多少。
那么计算答案是 $\sum_{i=n-1}^me(1,i)*P(i条边恰好使图联通)$
考虑计算这个e，写出转移：
$e(x,i)=1+\frac{1}{d_x}\sum_{x->y}p_x*e(y,i-1)+(1-p_x)*e(y,i)$
注意到无论i是多少都是这个转移式，这启示我们只需要弄出系数表达，即
$e(x,i)=b(x)+\sum_{y=1}^na(x,y)*e(y,i-1)$
我们希望求出a和b，这样就可以递推得出e。
a和b很好求，只需要我们一次高斯消元即可，这里是n^3，然后递推需要的复杂度是n*m*n=n^4。
因此这个部分只需要n^4。
现在我们需要解决如何求任意条边恰好使图联通的方案数，可以发现 $P(i条边恰好使图联通)=P(i条边使图联通)-P(i-1条边使图联通)$ 。
因此我们要统计n个点任意边的联通图个数。

n个点任意条边的联通图个数

注意现在不是任意点了，搬上面的做法dp来得到f此时只需要n^4的复杂度。
但是处理dp仍然是一个n^5的东西，考虑优化。
$dp(i,j)=\sum_{G\ is\ a\ graph,v(G)=i,e(G)=j}(-1)^{c(G)-1}*(c(G)-1)!$
观察这条式子，先想想如何计算一个图 $(c(G)-1)!$ 次，我们可以每次枚举任意一个不包含1号点的联通块，因此这个问题迎刃而解。
dp的转移就是每次枚举任意一个不包含1号点的联通块，复杂度n*m*n=n^4。
于是所有部分均在n^4复杂度内解决。

#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fd(i,a,b) for(i=a;i>=b;i--)
using namespace std;
typedef long long ll;
const int maxn=50+10,mo=1000000007;
int dp[maxn][maxn*maxn],f[maxn][maxn*maxn],ff[maxn][maxn*maxn],fff[maxn][maxn*maxn],g[maxn*maxn],c[maxn*maxn][maxn*maxn],a[maxn][maxn*2],d[maxn],p[maxn];
int e[maxn][maxn*maxn],num[maxn][maxn];
int h[maxn],go[maxn*maxn],nxt[maxn*maxn],fac[maxn*maxn],inv[maxn*maxn];
int i,j,k,l,r,t,n,m,ans,tot,top,q,qq,x,y,ni,mx;
bool czy;
int qsm(int x,int y){
    if (!y) return 1;
    int t=qsm(x,y/2);
    t=(ll)t*t%mo;
    if (y%2) t=(ll)t*x%mo;
    return t;
}
void add(int x,int y){
    d[x]++;
    go[++tot]=y;
    nxt[tot]=h[x];
    h[x]=tot;
}
int S(int n){
    return n*(n-1)/2;
}
int C(int n,int m){
    if (n<m||m<0) return 0;
    return (ll)fac[n]*inv[m]%mo*inv[n-m]%mo;
}
void gauss(){
    int i,j,k,r,t;
    fo(i,1,n){
        fo(j,i,n)
            if (a[j][i]) break;
        fo(k,i,2*n+1) swap(a[i][k],a[j][k]);
        //if (!a[i][i]) czy=1;
        fo(j,i+1,n){
            if (!a[j][i]) continue;
            t=(ll)a[j][i]*qsm(a[i][i],mo-2)%mo;
            fo(k,i,2*n+1){
                r=(ll)a[i][k]*t%mo;
                //a[j][k]=r-a[j][k]<0?r-a[j][k]+mo:r-a[j][k];
                a[j][k]=(r-a[j][k])%mo;
            }
            //b[j]=((ll)b[i]*t%mo-b[j])%mo;
        }
    }
    fd(i,n,1){
        fd(j,n,i+1) 
            fo(k,n+1,2*n+1) (a[i][k]-=(ll)a[i][j]*num[j][k-n]%mo)%=mo;
        fo(k,1,n+1) num[i][k]=(ll)a[i][k+n]*qsm(a[i][i],mo-2)%mo;
        //b[i]=0;
    }
}
int main(){
    //freopen("island5.in","r",stdin);
    scanf("%d",&n);
    m=n*(n-1)/2;
    mx=max(n,m);
    fo(i,1,n) scanf("%d",&p[i]);
    fo(i,1,n){
        scanf("%d",&t);
        fo(j,1,t){
            scanf("%d",&k);
            add(i,k);
        }
    }
    fac[0]=1;
    fo(i,1,mx) fac[i]=(ll)fac[i-1]*i%mo;
    inv[mx]=qsm(fac[mx],mo-2);
    fd(i,mx-1,0) inv[i]=(ll)inv[i+1]*(i+1)%mo;
    fo(i,1,n) dp[i][S(i)]=1;
    fo(i,0,n)
        fo(j,0,m)
            if (dp[i][j])
                fo(k,1,n-i)
                    if (j+(l=S(k))<=m){
                        t=(ll)dp[i][j]*C(i+k-1,k)%mo;   
                        dp[i+k][j+l]=dp[i+k][j+l]-t<0?dp[i+k][j+l]-t+mo:dp[i+k][j+l]-t;
                    }
    fo(i,n,n)
        fo(j,0,m)
            fo(k,j,m){
                t=(ll)dp[i][k]*C(k,j)%mo;
                f[i][j]=f[i][j]+t>=mo?f[i][j]+t-mo:f[i][j]+t;
            }
    fo(i,1,m) g[i]=(ll)f[n][i]*qsm(C(m,i),mo-2)%mo;
    fd(i,m,1) (g[i]-=g[i-1])%=mo;
    fo(x,1,n){
        a[x][x]++;
        a[x][2*n+1]++;
        ni=qsm(d[x],mo-2);
        t=h[x];
        while (t){
            y=go[t];
            (a[x][y+n]+=(ll)ni*p[x]%mo)%=mo;
            (a[x][y]-=(ll)ni*(1-p[x])%mo)%=mo;
            t=nxt[t];
        }
    }
    fo(i,1,n)
        fo(j,1,2*n+1) (a[i][j]+=mo)%=mo;
    gauss();
    fo(i,1,n) e[i][0]=0;
    fo(i,1,m){
        fo(j,1,n){
            fo(k,1,n) (e[j][i]+=(ll)num[j][k]*e[k][i-1]%mo)%=mo;
            (e[j][i]+=num[j][n+1])%=mo;
        }
    }
    fo(i,1,m) (ans+=(ll)e[1][i]*g[i]%mo)%=mo;
    ans--;
    (ans+=mo)%=mo;
    printf("%d\n",ans);
}

WerKeyTom_FTD

关注

1
点赞
踩
2

收藏

觉得还不错? 一键收藏
2
评论
[2017集训队作业自选题#149]小c的岛屿

前言感觉这题挺棒棒的。题目描述小c有n个岛屿。她为了管理岛上的人，决定让这些岛屿连通。每个岛屿i有一个概率值pi，和一个友好列表Ai。小c首先来到了1号岛屿，她会依次执行下面操作： 1、设现在在岛屿x，有px的概率产生一条无向边连接两个随机岛屿，这两个岛屿不会相同，也不会已经有边相连。（即在尚不存在的无向边中随机一条加入图中，不会加自环） 2、如果此时所有岛屿连通，她就会心满意足地离开。
复制链接

扫一扫