题解 P3307 【[SDOI2013]项链】

2019-01-14 11:02:22


丧心病狂的数论二合一题。。

可以分成两部分来考虑,珠子的种类数和组成环的方案数。


珠子的种类数

设 $s_3$为最多三元环的种类数,即:

$$s_3=\sum_{i=1}^a\sum_{j=1}^a\sum_{k=1}^a[\gcd(i,j,k)=1]$$

$s_2,s_1$依次类推,即: $$s_2=\sum_{i=1}^a\sum_{j=1}^a[\gcd(i,j)=1],s_1=1$$ 可以容斥得到:

总种类数 $=1+(3s_2-3)/3+(s_3-(3s_2-3)-1)/6=(s_3+3s_2+2)/6.$

然后莫比乌斯反演求一下 $s$就行了,式子: $$s_3=\sum_{d=1}^a\mu(d)\lfloor\frac{a}{d}\rfloor^3,s_2=\sum_{d=1}^a\mu(d)\lfloor\frac{a}{d}\rfloor^2$$


环的方案数

设先前求出来的方案数为 $m$。

看到环旋转前后属于同一种方案,可以想到 $polya$定理: $$ans=\sum_{i=1}^{n}f(\gcd(n,i))=\sum_{d|n}\varphi(\frac{n}{d})f(d)=\sum_{d|n}\varphi(d)f(\frac{n}{d})$$ 由于有限制:相邻的颜色不同,设 $f(n)$为 $n$个点的环满足条件的涂色方案,则有: $$f(n)=(m-2)f(n-1)+(m-1)f(n-2)$$ 这个式子可以理解为:枚举一个合法方案的断点,两边颜色一定不同,这时候要么当前方案是由 $(n-1)$个点组成的,只有 $(m-2)$种颜色可选,要么由 $(n-2)$个点组成的,加一个与一边相同的点,然后只有 $(m-1)$中颜色选。

然后发现这是一个二阶齐次线性方程,可以搬出特征方程的套路,设这是一个等比数列,可得特征方程: $$x^2=(m-2)x+m-1$$ 解得: $$x_1=m-1,x_2=-1$$ 设当前数列通项公式为: $$f(n)=\alpha x_1^n+\beta x_2^n$$ 由于: $$f(1)=0,f(2)=m^2-m$$ 列出方程:

$$\alpha(m-1)-\beta=0$$ $$\alpha(m-1)^2+\beta=m^2-m$$

解得 $\alpha=1,\beta=m-1$,所以 $f$的通项公式为: $$f(n)=(m-1)^n+(-1)^n(m-1)$$ 答案是: $$\sum_{d|n}\varphi(\frac{n}{d})f(d)$$所以直接一路暴力就好了。

记得 $long\,\,long$会乘爆,用爆 $long\,\,long$小技巧就好了,具体看代码 $mul$函数。

由于 $n$有可能是 $1e9+7$的倍数,所以特判下,如果是就模 $(1e9+7)^2$,最后除掉 $1e9+7$即可。

然后细节还有挺多的,,注意下写法,不然调试很恶心。

#include<bits/stdc++.h>
using namespace std;

#define ll long long 
#define lf long double
#define ull unsigned long long

void read(ll &x) {
    x=0;ll f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

#define sqr(x) mul(x,x)
#define cul(x) mul(sqr(x),x)

const int maxn = 1e7+10;
const int MOD = 1e9+7;
const int MAXN = 1e4+10;

ll mod;

int pri[maxn],vis[maxn],mu[maxn],tot;

void sieve(int N) {
    mu[1]=1;
    for(int i=2;i<=N;i++) {
        if(!vis[i]) pri[++tot]=i,mu[i]=-1;
        for(int j=1;j<=tot&&i*pri[j]<=N;j++) {
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=N;i++) mu[i]+=mu[i-1];
}

ll mul(ll a,ll b) {
    lf res=(lf)a*b/mod;
    ll t=a*b,t2=(t-(ull)res*mod+mod)%mod;return t2;
}

ll qpow(ll a,ll x) {
    ll res=1;
    for(;x;x>>=1,a=mul(a,a)) if(x&1) res=mul(res,a);
    return res;
}

ll a[MAXN],p[MAXN],cnt;

ll m,ans,inv6,nown;

void mobius(int n) {
    m=0;int T=1;
    while(T<=n) {
        int pre=T;T=n/(n/T);
        m=(m+1ll*mul((mu[T]-mu[pre-1]),cul(n/T)))%mod;
        m=(m+1ll*mul((mu[T]-mu[pre-1]),sqr(n/T))*3ll%mod)%mod;T++;
    }m=(m+2)%mod;m=mul(m,inv6)%mod;
}

void prepare(ll n) {
    cnt=0;
    for(int i=2;1ll*i*i<=n;i++) {
        if(n%i) continue;
        a[++cnt]=i;p[cnt]=0;
        while(n%i==0) n/=i,p[cnt]++;
    }
    if(n!=1) a[++cnt]=n,p[cnt]=1;
}

ll f(ll n) {
    ll Ans=qpow(m-1,n);
    if(n&1) Ans=(Ans+mod-m+1)%mod;
    else Ans=(Ans+m-1)%mod;
    return Ans;
}

void polya(int now,ll d,ll phi) {
    if(now==cnt+1) return ans=(ans+mul(phi,f(nown/d)))%mod,void();
    polya(now+1,d,phi);
    d*=a[now],phi*=a[now]-1;
    polya(now+1,d,phi);
    for(int i=2;i<=p[now];i++)
        d*=a[now],phi*=a[now],polya(now+1,d,phi);
}

ll inn[20],ina[20],mx,inv[10];

void pre_inv() {
    inv[1]=1;
    for(int i=2;i<=6;i++) inv[i]=mul(mod-mod/i,inv[mod%i]);
    inv6=inv[6];
}

int main() {
    int t;scanf("%d",&t);
    for(int i=1;i<=t;i++) read(inn[i]),read(ina[i]),mx=max(mx,ina[i]);
    sieve(mx);
    for(int i=1;i<=t;i++) {
        nown=inn[i];
        if(inn[i]%MOD) mod=MOD;
        else mod=1ll*MOD*MOD;
        pre_inv();
        mobius(ina[i]);
        prepare(inn[i]);ans=0;
        polya(1,1,1);
        if(inn[i]%MOD) ans=mul(ans,qpow(inn[i],mod-2));
        else ans=mul(ans/MOD,qpow(inn[i]/MOD,MOD-2))%MOD;
        write(ans);
    }
    return 0;
}