JZOJ 3984 宝石纪念币

宝石纪念币

Description

给出一个有$n$个点的圆环，现在要用$17$种不同的颜色填充它。
若两种填色方案经过旋转后可以重合，则视这两种填色方案为同一种填色方案。
保证$17$种颜色通必须用上。
求方案数，答案保留最后$120$位，不够的用$0$补齐。

Data Constraint

$n$<=$10^9$

Solution

将这$n$种旋转方法视为$n$种不同的置换。
对于第$i$种置换（向右旋转$i$位的置换），通过简单的数学推导可以得出单个循环节的长度为$n \over {gcd(i,n)}$，故循环节数为${gcd(i,n)}$。
设仅用$m$种颜色染色的方案数为$S_m$。
由$Burnside$引理和$Polya$定理可得

$$S_m={\sum_{i=1}^n m^{gcd(i,n)}\over n}$$
考虑枚举 $gcd$( $i$,$n$)，可得
$$S_m*n=\sum_{i|n} m^{i} \sum_{j=1}^n[gcd(j,n)=i]$$
也就是
$$S_m*n=\sum_{i|n} m^{i} \sum_{j=1}^{n\over i}[gcd(j*i,n)=i]$$
$$S_m*n=\sum_{i|n} m^{i} \sum_{j=1}^{n\over i}[gcd(j,{n\over i})=1]$$
$$S_m*n=\sum_{i|n} m^{i} \varphi ({n\over i})$$
考虑到 $n$的约数不会有很多，因此可以暴力分解质因数然后暴力枚举约数，求解即可。
Ps:记得要用高精度。

接下来看到条件$17$种颜色都要用到，考虑容斥

$$Answer=\sum_{m=1}^{17}S_m*(-1)^{(i+1)}*C_{17}^m$$
接下来只需要求 $S_m$这道题就做完了。
把数转成高精度后我们只能求出 $S_m*n$，我们最后还要除以 $n$，这就比较难处理了。
于是我们可以把每一个高精度数保存为$q$* $n$+$p$的形式，其中 $p$为这个高精度数$mod$ $n$的余数，$q$为令一个高精度数，换言之，就是把最后一位的十进制改成 $n$进制，最后答案直接输出$q$部分即可，要进行压位操作卡常。

Ps:丧心病狂的出题人为了增加这道群论裸题的难度加了这么恶心的操作，细节这么多，害我比赛时调了好久，不可理喻(•́へ•́╬)

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>

#define fo(i,j,l) for(int i=j;i<=l;++i)
#define fd(i,j,l) for(int i=j;i>=l;--i)

using namespace std;
typedef long long ll;
const ll N=35,M=2e3,jw=1e8,sx=25,Q=40;

int oo,n,m,lp,j,k,l,i,o;
ll ss[N][2];
ll cj[M],lj[M],ans[M],t[M],re[M];
ll ys[M],phi[M];

inline ll gcd(ll a,ll b)
{return b==0?a:gcd(b,a%b);}

inline int min(int a,int b)
{return a<b?a:b;}

inline int max(int a,int b)
{return a>b?a:b;}

void dg(int o,ll p)
{
    if(o==oo+1){
        ys[++lp]=p;
        p=n/p;
        fo(i,1,oo)if(p%ss[i][0]==0)p=(p/ss[i][0])*(ss[i][0]-1);
        phi[lp]=p;
        return;
    }
    fo(i,0,ss[o][1])dg(o+1,p),p*=ss[o][0];
}

inline void times(ll *a,ll p)
{
    a[Q]=a[Q]*p; 
    ll cqy=a[Q]/n; a[Q]=a[Q]%n;
    fo(i,1,a[0])
    {
        a[i]=cqy+a[i]*p;
        cqy=a[i]/jw;
        a[i]=a[i]%jw;
    }
    while(cqy)a[++a[0]]=cqy%jw,cqy/=jw;
    while(a[0]>sx)a[a[0]]=0,--a[0];
}

inline void twotimes(ll *a,ll *b,ll *c)
{
    fo(i,1,sx+1)t[i]=0; 
    t[0]=1; 
    fo(i,1,min(sx,a[0]))
    {
        fo(l,1,min(b[0],sx-i+1))
        t[i+l-1]=t[i+l-1]+a[i]*b[l];
    }
    t[0]=sx; ll cqy=0;
    fo(i,1,sx)t[i]=cqy+t[i],cqy=t[i]/jw,t[i]=t[i]%jw;
    fo(i,a[0]+1,sx)a[i]=0; 
    fo(i,b[0]+1,sx)b[i]=0; 
    fo(i,1,sx)t[i]=t[i]*n+a[i]*b[Q]+b[i]*a[Q];
    t[Q]=a[Q]*b[Q]; 
    t[1]=t[1]+t[Q]/n; 
    t[Q]=t[Q]%n;
    t[0]=sx; cqy=0;
    fo(i,1,sx)t[i]=cqy+t[i],cqy=t[i]/jw,t[i]=t[i]%jw;
    while(t[t[0]]==0&&t[0]!=0)--t[0];
    fo(i,0,t[0])c[i]=t[i]; c[Q]=t[Q];
}

inline void pplus(ll *a,ll *b,ll *c)
{
    t[Q]=a[Q]+b[Q]; ll cqy=t[Q]/n; t[Q]=t[Q]%n;
    t[0]=min(sx,max(a[0],b[0]));
    fo(i,1,t[0])t[i]=a[i]+b[i]+cqy,cqy=t[i]/jw,t[i]=t[i]%jw;
    while(cqy)t[++t[0]]=cqy,cqy/=jw;
    while(t[0]>sx)t[t[0]]=0,--t[0];
    fo(i,0,t[0])c[i]=t[i]; c[Q]=t[Q];
}

inline void mminus(ll *a,ll *b,ll *c)
{
    fo(i,0,sx+1)t[i]=0;
    t[0]=min(sx,max(a[0],b[0]));
    fo(i,1,t[0])t[i]=a[i]-b[i];
    t[Q]=a[Q]-b[Q]; if(t[Q]<0)t[Q]+=n,--t[1];
    fo(i,1,t[0])if(t[i]<0)t[i]=t[i]+jw,--t[i+1];
    while(t[t[0]+1]<0&&t[0]<sx)t[0]++,t[t[0]]+=jw,--t[t[0]+1];
    fo(i,0,t[0])c[i]=t[i]; c[Q]=t[Q];
}

int main()
{
    cin>>n; ll sq=(ll)sqrt(n);
    oo=0; ll last=n;
    fo(i,2,sq)
    if(last%i==0){
        ss[++oo][0]=i;
        while(last%i==0)++ss[oo][1],last/=i;
    }
    if(last>1)ss[++oo][0]=last,ss[oo][1]=1;
    dg(1,1);
    ll u=1;
    fd(i,17,1){
        ll xs=1;
        fo(l,17-i+1,17)xs=xs*l;
        fo(l,1,i)xs=xs/l;
        re[0]=1; fo(l,1,sx+1)re[l]=0; re[Q]=0;
        fo(l,1,lp){
            lj[Q]=i; lj[1]=cj[1]=0; lj[0]=cj[0]=1; cj[Q]=1;
            int zs=ys[l];
            for(;zs;zs>>=1,twotimes(lj,lj,lj))
            if(zs&1)twotimes(lj,cj,cj);
            times(cj,phi[l]);
            pplus(re,cj,re);
        }
        times(re,xs);
        if(u==1)pplus(ans,re,ans);else mminus(ans,re,ans);
        u*=-1;
    }
    fd(i,15,1){
        int u=ans[i],k=0;
        while(u)u/=10,++k;
        fo(l,1,8-k)putchar('0');
        if(ans[i])printf("%lld",ans[i]);
    }
}