【THUSC2017】【LOJ2981】若是奇蹟有顏色 DP BM 打表 線性遞推

題目大意

　　有一個 \(n\) 個點的環，你要用 \(m\) 中顏色染這 \(n\) 個點。

　　要求連續 \(m\) 個點的顏色不能是 $1 \sim m $ 的排列。

　　兩種環相同當且僅當這兩個環能夠在旋轉以後變得如出一轍。

　　求方案數對 \({10}^9+7\) 取模的結果。

　　\(n\leq {10}^9,m\leq 7\)

題解

　　考慮 polya 定理，記 \(f(n)\) 爲 \(n\) 個點的答案，\(g(n)\) 爲 \(n\) 個點不考慮旋轉的答案。那麼就有
\[ \begin{align} f(n)&=\frac{1}{n}\sum_{i=1}^ng(\gcd(n,i))\\ &=\frac{1}{n}\sum_{i\mid n}\varphi(\frac{n}{i})g(i) \end{align} \]
　　\(g(i)\) 能夠 DP 計算。

　　記 \(h_{i,j,k}\) 爲長度爲 \(i\) 的鏈，前 \(m\) 個點的狀態（顏色）爲 \(j\)，最後 \(m\) 個點的狀態爲 \(k\) 的方案數。

　　還能夠記錄前 \(m\) 個顏色的前多少個是互不相同的，還有最後 \(m\) 個點的顏色。就前 \(j\) 個的顏色是互不相同的，第 \(j+1\) 個點顏色和前面某個點顏色相同。

　　顯然 \(g(i)\) 有一個長度不超過 \(m^m\) 的遞推式。

　　暴力打出前面的項而後 BM 求出遞推式便可。

　　開 O2 只用了 48s 就打出來了。

　　\(m=7\) 時遞推式長度爲 \(409\)。

　　表在這：https://github.com/ywwywwyww/THUSC2017/tree/master/yww/farben

代碼

const int N=1010;
const ll p=1000000007;
ll fp(ll a,ll b)
{
    ll s=1;
    for(;b;b>>=1,a=a*a%p)
        if(b&1)
            s=s*a%p;
    return s;
}

ll b[10][2010]=表;
ll a[10][2010]=表;

int c[N],d[N],t;
ll pw[N][N];
ll ans;
int n,m;
int s[N];
int len;
void add(int &a,ll b)
{
    a=(a+b)%p;
}

void mul()
{
    static int f[N];
    for(int i=0;i<=2*len;i++)
        f[i]=0;
    for(int i=0;i<=len;i++)
        if(s[i])
            for(int j=0;j<=len;j++)
                add(f[i+j],(ll)s[i]*s[j]);
    for(int i=0;i<=2*len;i++)
        s[i]=f[i];
}
void module()
{
    for(int i=2*len;i>=len;i--)
        if(s[i])
        {
            for(int j=1;j<=len;j++)
                add(s[i-j],(ll)s[i]*b[m][j]);
            s[i]=0;
        }
}

void fp(int x)
{
    if(!x)
    {
        for(int i=0;i<=2*len;i++)
            s[i]=0;
        s[0]=1;
        return;
    }
    fp(x>>1);
    mul();
    if(x&1)
    {
        for(int i=2*len;i>=1;i--)
            s[i]=s[i-1];
        s[0]=0;
    }
    module();
}

ll calc(int x)
{
    if(x<=500)
        return a[m][x];
    fp(x-1);
    ll res=0;
    for(int i=1;i<=len;i++)
        res=(res+(ll)a[m][i]*s[i-1])%p;
    return res;
}

void dfs(int x,int y,ll phi)
{
    if(x>t)
    {
        ans=(ans+calc(y)*phi)%p;
        return;
    }
    for(int i=0;i<d[x];i++)
        dfs(x+1,y*pw[x][i],phi*(c[x]-1)%p*pw[x][d[x]-i-1]%p);
    dfs(x+1,y*pw[x][d[x]],phi);
}
int main()
{
    open("farben");
    scanf("%d%d",&n,&m);
    int _n=n;
    for(int i=2;i*i<=_n;i++)
        if(_n%i==0)
        {
            c[++t]=i;
            while(_n%i==0)
            {
                d[t]++;
                _n/=i;
            }
        }
    if(_n>1)
    {
        c[++t]=_n;
        d[t]=1;
    }
    for(int i=1;i<=t;i++)
    {
        pw[i][0]=1;
        for(int j=1;j<=d[i];j++)
            pw[i][j]=pw[i][j-1]*c[i]%p;
    }
    len=b[m][0];
    dfs(1,1,1);
    ans=ans*fp(n,p-2)%p;
    ans=(ans%p+p)%p;
    printf("%lld\n",ans);
    return 0;
}

打表程序

const ll p=1000000007;
const int N=1010;
const int n=1000;
//const int m=5;
int m;

ll fp(ll a,ll b)
{
    ll s=1;
    for(;b;b>>=1,a=a*a%p)
        if(b&1)
            s=s*a%p;
    return s;
}


const int MA=2100000;

int e[MA];
int ban[MA];
vector<int> g[N];
int pw[N];
int ma;

namespace gao1
{
    int a[N],b[N],c[N];
    void dfs(int x)
    {
        if(x>m)
        {
            for(int i=1;i<=m;i++)
                b[i]=0;
            int tot=0;
            int s=0;
            int s2=0;
            int first=0;
            for(int i=1;i<=m;i++)
            {
                if(!b[a[i]])
                    b[a[i]]=++tot;
                c[i]=b[a[i]];
                s+=c[i]*pw[i-1];
                s2+=a[i]*pw[i-1];
                if(c[i]<=c[i-1]&&!first)
                    first=i-1;
            }
            if(tot==m)
                ban[s2]=1;
            else
                g[first].push_back(s);
            return;
        }
        for(int i=1;i<=m;i++)
        {
            a[x]=i;
            dfs(x+1);
        }
    }
    void gao()
    {
        dfs(1);
    }
}

namespace gao2
{
    int a[N],b[N],c[N];
    void dfs(int x)
    {
        if(x>m)
        {
            for(int i=1;i<=m;i++)
                b[i]=0;
            for(int i=m;i>=1;i--)
                if(!a[i]||!b[a[i]])
                {
                    c[i]=a[i];
                    b[a[i]]=1;
                }
                else
                    c[i]=0;
            int s=0,s1=0;
            for(int i=1;i<=m;i++)
                s+=c[i]*pw[i-1];
            for(int i=1;i<=m;i++)
                s1+=a[i]*pw[i-1];
            e[s1]=s;
            return;
        }
        for(int i=0;i<=m;i++)
        {
            a[x]=i;
            dfs(x+1);
        }
    }
    void gao()
    {
        dfs(1);
    }
}

int f[N];
int h[2][MA];
int d[MA];

void add(int &a,int b)
{
    a=(a+b)%p;
}

int append(int a,int b)
{
    return a/(m+1)+b*pw[m-1];
}

namespace gao3
{
    void gao(int x)
    {
        memset(d,0,sizeof d);
        
        for(int i=0;i<=ma;i++)
        {
            int flag=1;
            for(int y=i,j=1;j<=x;j++)
            {
                y=append(y,j);
                if(ban[y])
                {
                    flag=0;
                    break;
                }
            }
            d[i]=flag;
        }
        
        memset(h,0,sizeof h);
        int cur=0;
        for(auto v:g[x])
            h[cur][e[v]]++;
        for(int i=m;i<=n;i++)
        {
            fprintf(stderr,"%d %d %d\n",m,x,i);
            memset(h[cur^1],0,sizeof h[cur^1]);
            for(int j=0;j<=ma;j++)
                if(h[cur][j]&&!ban[j])
                {
                    add(f[i],d[j]*h[cur][j]);
                    for(int k=1;k<=m;k++)
                        add(h[cur^1][e[append(j,k)]],h[cur][j]);
                }
            cur^=1;
        }
    }
}

int main(int argc,char **argv)
{
//  freopen("farben2.txt","w",stdout);
    sscanf(argv[1],"%d",&m);
    ma=fp(m+1,m);
    pw[0]=1;
    for(int i=1;i<=m;i++)
        pw[i]=pw[i-1]*(m+1);
    gao1::gao();
    gao2::gao();
    
    for(int i=1;i<m;i++)
        f[i]=fp(m,i);
    
    for(int i=1;i<m;i++)
        gao3::gao(i);
    
    for(int i=1;i<=n;i++)
        printf("%d\n",f[i]);
    return 0;
}