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

版权声明:本文是蒟蒻写出来的,神犇转载也要说一声哦! https://blog.csdn.net/WerKeyTom_FTD/article/details/78947101

前言

感觉这题挺棒棒的。

题目描述

小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)=px(1q)1dxx>y(e(y,i+1)+1)+(1px)1dxx>y(e(y,i)+1)
我们分层来做,可以注意到如果q能够求出,每层内可以高斯消元。
那么复杂度是m*n^3=n^5。
现在问题是这个q怎么求。
显然q=i+1(n(n1)/2i)i
我们来考虑怎么求得分母和分子。

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

假设要求m条边的,现在一条桥边把图分成两个联通块,我们不妨枚举1号点所在联通块有k个点,l条边,那么贡献:
Ck1n1k(nk)(kl)(nkml1)
其中组合数表示除了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)=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)=i>=mdp(n,i)Cmi
我们先考虑如何计算dp,可以设一个g。
g(i,j,k)=G is a graph,v(G)=i,e(G)=j,c(G)=k(1)k1(k1)!
这个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条边,期望步数是多少。
那么计算答案是mi=n1e(1,i)P(i使)
考虑计算这个e,写出转移:
e(x,i)=1+1dxx>ypxe(y,i1)+(1px)e(y,i)
注意到无论i是多少都是这个转移式,这启示我们只需要弄出系数表达,即
e(x,i)=b(x)+ny=1a(x,y)e(y,i1)
我们希望求出a和b,这样就可以递推得出e。
a和b很好求,只需要我们一次高斯消元即可,这里是n^3,然后递推需要的复杂度是n*m*n=n^4。
因此这个部分只需要n^4。
现在我们需要解决如何求任意条边恰好使图联通的方案数,可以发现P(i使)=P(i使)P(i1使)
因此我们要统计n个点任意边的联通图个数。

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

注意现在不是任意点了,搬上面的做法dp来得到f此时只需要n^4的复杂度。
但是处理dp仍然是一个n^5的东西,考虑优化。
dp(i,j)=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);
}
阅读更多
想对作者说点什么?

博主推荐

换一批

没有更多推荐了,返回首页