bzoj3899: 仙人掌树的同构

Description

首先，先介绍仙人掌树。仙人掌树是一张无向图，但是每个节点最多只会在一个环里面，而且这张图的环全部都是简单环，即A->B->C->A这种。

比如下图就是一颗仙人掌树。

好的，知道了仙人掌树之后，我们现在要计算一个东西。

我们现在已经知道了一个N个节点的仙人掌树，称作为原图。接下来，我们要用1-N的一个排列A[1]-A[N]去变换这棵树，具体的，如果原图中有一条边i-j，那么变换出来的图中必须有一条A[i]-A[j]的边。同样的，如果变换出来的图中有一条A[i]-A[j]的边，那么原图中必有一条i-j的边。（简单而言就是点重新编号）

小A为了超脱宇宙的束缚，必须要知道，有多少种排列，可以使得变换出来的新图和原图是一模一样的，具体的，原图中如果存在一条i-j的边，新图也存在一条i-j的边，新图中存在一条i-j的边，原图中也存在i-j的边。

方案数目答案mod 1000000003。

Input

第一行有两个正整数，N和M，节点个数和边的个数。

接下来M行，每行有2个正整数S，T，表示一条原图的无向边。数据保证没有重边。

Output

一行一个正整数表示方案书目。

转成圆方树的自同构，找到重心变为有根树，对一个圆点，其k棵子树相同对答案贡献为k!，对一个方点，若为根直接考虑旋转和翻转的每种情况是否与原先同构，否则只考虑翻转后是否同构，贡献为2，全部用hash判定并将贡献累乘即可。方点相邻的点是有序的，在计算时要把相邻点排成环上的顺序的再处理。

#include<cstdio>
#include<algorithm>
typedef unsigned long long u64;
const int N=20007,P=1000000003;
int n,m,fac[N];
int es[N],enx[N],e0[N],e1[N],ep=2;
bool _c[N];
void ae(int*e,int a,int b){
    es[ep]=b;enx[ep]=e[a];e[a]=ep++;
    es[ep]=a;enx[ep]=e[b];e[b]=ep++;
}
void de(int*e,int a,int b){
    for(int*i=e+a;*i;i=enx+*i){
        int u=es[*i];
        if(u==b){
            *i=enx[*i];
            return;
        }
    }
}
int dfn[N],low[N],tk=0,ss[N],sp=0,n1,ri[N];
void tj(int w,int pa){
    dfn[w]=low[w]=++tk;
    ss[++sp]=w;
    for(int i=e0[w];i;i=enx[i]){
        int u=es[i];
        if(u==pa)continue;
        if(!dfn[u]){
            tj(u,w);
            if(ss[sp]==u)--sp,ae(e1,w,u);
        }else if(dfn[u]<dfn[w]){
            int rp=0;
            ae(e1,u,++n1);
            _c[n1]=1;
            while(dfn[ss[sp]]>dfn[u])ri[ss[sp]]=++rp,ae(e1,ss[sp--],n1);
        }
    }
}
int sz[N],cg[3],cgp=0,rt;
void f1(int w,int pa){
    bool is=1;
    sz[w]=1;
    for(int i=e1[w];i;i=enx[i]){
        int u=es[i];
        if(u==pa)continue;
        f1(u,w);
        sz[w]+=sz[u];
        if(sz[u]*2>n)is=0;
    }
    if(sz[w]*2<n)is=0;
    if(is)cg[cgp++]=w;
}
u64 h[N],h1[N],h2[N];
int cs[N],cp,ans=1;
bool cmp_h(int a,int b){return h[a]<h[b];}
bool cmp_ri(int a,int b){return ri[a]<ri[b];}
void f2(int w,int pa){
    for(int i=e1[w];i;i=enx[i]){
        int u=es[i];
        if(u!=pa)f2(u,w);
    }
    cp=0;
    for(int i=e1[w];i;i=enx[i]){
        int u=es[i];
        if(u!=pa)cs[cp++]=u;
    }
    if(_c[w]){
        if(w!=rt)
        for(int i=0;i<cp;++i)ri[cs[i]]=(ri[cs[i]]-ri[pa]+(cp+1))%(cp+1);
        std::sort(cs,cs+cp,cmp_ri);
        if(w==rt){
            u64 p0=29399999,pp=1;
            for(int i=1;i<=cp;++i){
                h1[i]=h1[i+cp]=h2[i]=h2[i+cp]=h[cs[i-1]];
                pp*=p0;
            }
            for(int i=1;i<=cp*2;++i)h1[i]+=h1[i-1]*p0;
            for(int i=cp*2;i;--i)h2[i]+=h2[i+1]*p0;
            int c=0;
            for(int i=cp;i<cp*2;++i)c+=(h1[i]-h1[i-cp]*pp==h1[cp]);
            for(int i=cp;i;--i)c+=(h2[i]-h2[i+cp]*pp==h1[cp]);
            ans=u64(ans)*c%P;
        }else{
            u64 h1=0,h2=0;
            for(int i=0;i<cp;++i)h1=h1*1399913+h[cs[i]];
            for(int i=cp-1;i>=0;--i)h2=h2*1399913+h[cs[i]];
            if(h1==h2)ans=ans*2%P;
            h[w]=std::min(h1,h2);
            h[w]^=h[w]>>13^h[w]*17<<31^41546541635416351llu;
        }
    }else{
        std::sort(cs,cs+cp,cmp_h);
        h[w]=0;
        for(int i=0,j=0;i<cp;i=j){
            for(;j<cp&&h[cs[i]]==h[cs[j]];++j);
            ans=u64(ans)*fac[j-i]%P;
        }
        for(int i=0;i<cp;++i)h[w]=h[w]*1844677+h[cs[i]];
        h[w]^=h[w]>>17^h[w]*11<<23^12218653252152541llu;
    }
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=fac[0]=1;i<=n;++i)fac[i]=fac[i-1]*u64(i)%P;
    n1=n;
    for(int i=0,a,b;i<m;++i){
        scanf("%d%d",&a,&b);
        ae(e0,a,b);
    }
    tj(1,0);
    f1(1,0);
    if(cgp==2){
        if(_c[cg[0]])rt=cg[0];
        else if(_c[cg[1]])rt=cg[1];
        else{
            rt=++n1;
            ae(e1,n1,cg[0]);
            ae(e1,n1,cg[1]);
            de(e1,cg[0],cg[1]);
            de(e1,cg[1],cg[0]);
        }
    }else rt=cg[0];
    f2(rt,0);
    printf("%d\n",ans);
    return 0;
}