题解-SDOI2013 项链

题面

SDOI2013 项链

(T) 组测试数据,每次给 (n,a)。有一个长度为 (n) 的珍珠项链,每个珍珠上有 (3) 个数构成一个环。旋转、翻转同构的珍珠相同。相邻的珍珠必须不同,旋转同构的项链相同。求本质不同的项链个数。答案对 (10^9+7) 取模。

数据范围:(1le nle 10^{14})(1le ale 10^7)


题解

话说最近我写的题怎么都这么大啊,动不动就 (4k+),而且好多是紫题,这道短一点却是黑题。

分成两步计算:计算有几种不同的珍珠,有几种不同的项链。


计算有几种不同的珍珠

(g(i)) 表示 (3) 个数 (gcd)(i) 的倍数的方案数。

(f(i))(3) 个数的 (gcd) 恰好为 (i) 的方案数。

[g(i)=sum_{i|d} f(d)Longleftrightarrow f(i)=sum_{i|d} g(d)mu(frac{d}{i}) ]

[ herefore f(1)=sum_{d=1}^{a} g(d)mu(d) ]

计算 (g(d)),相当于 (3) 个数的选取范围是 (lfloorfrac ad floor)

利用 Burnside 定理,考虑 (G) 的元素,设 (m=lfloorfrac ad floor)

[g(d)=frac{1}{6}left(2m+3m^2+m^3 ight) ]

然后就可以计算出 (f(1)) 了,设 (c=f(1)),需要线性筛和整除分块。


计算有几种不同的项链

同样利用 Burnside 定理,枚举循环节。

[ans=frac{1}{n}sum_{d|n} p(d)varphi(frac nd) ]

先想想这个 (varphi(d)) 怎么做?不能线性筛,直接求也会 TLE

可以利用 (n) 的每个质因数,统一计算 (varphi(d)=prod_{p^c|d} left(p^c-p^{c-1} ight))

这东西可以把 (n) 的质因数和幂次预处理出来,用一个 dfs 实现。

这个 (p(d)) 是长度为 (d) 的鸽子序列,每个鸽子可以染 (c) 种颜色,相邻、首尾两只鸽子不能同色的方案数。

我曾经犯过的错误(点击展开)。

[p(d)=c(c-1)^{d-2}(c-2) ]

分类讨论:第 (1) 个鸽子和第 (d-1) 个鸽子颜色是否一样。

  • 如果一样,相当于 (d-2) 个鸽子,最后一只切成两只,然后往中间放第 (d) 个鸽子:(p(d)leftarrow (c-1)p(d-2))

  • 如果不一样,相当于 (d-1) 只鸽子,在最后两只之间放第 (d) 个鸽子:(p(d)leftarrow (c-2)p(d-1))

边界情况:(p(1)=0)(p(2)=c(c-1))

矩阵快速幂会 TLE,但是这个东西可以通过特征方程递推:

[p(d)=(c-2)p(d-1)+(c-1)p(d-2) o x^2=(c-2)x+c-1 ]

[(x-(c-1))(x+1)=0 o x=c-1 { m or} -1 ]

[egin{cases} p(d)-(c-1)p(d-1)=-(p(d-1)-(c-1)p(d-2))\ p(d)+p(d-1)=(c-1)p(d-1)+(c-1)p(d-2)\ end{cases}]

[egin{cases} p(d)-(c-1)p(d-1)=(-1)^{d-2}c(c-1)\ p(d)+p(d-1)=(c-1)^{d-2}c(c-1)\ end{cases}]

[p(d)=(c-1)((-1)^{d-2}+(c-1)^{d-1}) ]


最后不可忽略的细节!

注意了,还有 ((10^9+7)|n) 的情况,要把 (mod) 变成 ((10^9+7)^2),这样 (mod mid n)

所有的乘法变成龟速乘。最后先把答案除以 (10^9+7)(因为答案如果不取模一定是 (n) 的倍数),再除以 (frac{n}{10^9+7}) 对于 (10^9+7) 的逆元。

时间复杂度 (Theta(T(sqrt alog mod+sqrt n+{ m maxd}(n)log^2 mod)))


代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
#define x first
#define y second
#define bg begin()
#define ed end()
#define pb push_back
#define mp make_pair
#define sz(a) int((a).size())
#define R(i,n) for(int i(0);i<(n);++i)
#define L(i,n) for(int i((n)-1);i>=0;--i)
const int iinf=0x3f3f3f3f;
const ll linf=0x3f3f3f3f3f3f3f3f;

//Math
const ll mod1=1e9+7;
const ll mod2=mod1*mod1;
ll mod,i6;
ll mymul(ll a,ll x,ll mod=::mod,ll res=0){
    for(;x;x>>=1,(a+=a)%=mod)
        if(x&1) (res+=a)%=mod;
    return res;
}
ll mypow(ll a,ll x,ll mod=::mod,ll res=1){
    for(;x;x>>=1,a=mymul(a,a,mod))
        if(x&1) res=mymul(res,a,mod);
    return res;
}
const ll i61=mypow(6,mod1-2,mod1);
const ll i62=mypow(6,(mod1-1)*mod1-1,mod2);

//Math//Number
namespace nb{
    const int N=1e7+1;
    bool np[N];
    int mu[N],sm[N|1];
    vector<int> prime;
    void sieve(int n=N){
        np[1]=true,mu[1]=1;
        for(int i=2;i<n;++i){
            if(!np[i]) prime.pb(i),mu[i]=-1;
            for(int p:prime){
                if(i*p>=n) break;
                np[i*p]=true;
                if(i%p==0){mu[i*p]=0; break;}
                mu[i*p]=-mu[i];
            }
        }
        R(i,n) sm[i+1]=sm[i]+mu[i];
    }
    vector<pair<ll,int>> pfac(ll n){
        vector<pair<ll,int>> res;
        for(ll d=2;d*d<=n;++d)if(n%d==0)
            for(res.pb({d,0});n%d==0;++res.back().y,n/=d);
        if(n>1) res.pb({n,1});
        return res;
    }
    void dfs(vector<pair<ll,int>> &p,
        vector<pair<ll,ll>> &np,int dep,pair<ll,ll> now){
        if(dep==sz(p)) return np.pb(now);
        pair<ll,ll> de=now;
        R(i,p[dep].y+1){
            dfs(p,np,dep+1,de);
            if(i) de.x*=p[dep].x,de.y*=p[dep].x;
            else de.x*=p[dep].x,de.y*=p[dep].x-1;
        }
    }
    vector<pair<ll,ll>> facphi(ll n){
        vector<pair<ll,ll>> res;
        vector<pair<ll,int>> pf=pfac(n);
        dfs(pf,res,0,mp(1,1));
        return res;
    }
}

//Matrix//Function
ll mycalc(ll c,ll n){
    if(n==1) return 0;
    if(n==2) return mymul(c,c-1+mod);
    c=(c-1+mod)%mod;
    ll res=mypow(c,n);
    if(n&1) (res+=mod-c)%=mod;
    else (res+=c)%=mod;
    return res;
}

//Main//Data
unordered_map<ll,ll> phi;

//Main
void mymain(){
    ll n,c=0,res=0; int a; cin>>n>>a;
    if(n%mod1==0) mod=mod2,i6=i62;
    else mod=mod1,i6=i61;
    // cout<<"mod="<<mod<<'
';
    for(int l=1,r=-1;l<=a;l=r){
        int m=a/l; ll g=2*m; r=a/m+1;
        g=(mymul(3,mymul(m,m))+g)%mod;
        g=(mymul(mymul(m,m),m)+g)%mod;
        g=mymul(mymul(g,i6),nb::sm[r]-nb::sm[l]+mod);
        (c+=g)%=mod;
    }
    // cout<<"c="<<c<<'
';
    vector<pair<ll,ll>> fp=nb::facphi(n);
    for(auto &u:fp) phi[u.x]=u.y;
    for(auto &u:fp) res=(mymul(mycalc(c,u.x),phi[n/u.x])+res)%mod;
    // cout<<"res="<<res<<'
';
    if(n%mod1==0) res=(res/mod1*mypow(n/mod1,mod1-2,mod1))%mod1;
    else (res*=mypow(n,mod1-2,mod1))%=mod1;
    cout<<res<<'
';
}
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0),cout.tie(0);
    nb::sieve();
    int t; for(cin>>t;t--;mymain());
    return 0;
}

祝大家学习愉快!

原文地址:https://www.cnblogs.com/George1123/p/14277797.html