bzoj3899: 仙人掌树的同构

#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;
}