Luogu3307 项链 | Arahc's Home

题面

简化题意

定义一个珠子为无序三元组 $(x,y,z)$ ， $0<x,y,z\leq a,\gcd(x,y,z)=1$ 。

定义一个项链为一个 $n$ 个珠子组成的环，相邻的珠子不相同。

求有多少本质不同（经过旋转能完全相同的定义为本质相同）的项链。多测， $n\leq 10^{14},a\leq 10^7,T\leq 10$ 。

题解

因为“项链”的要求只需要相邻珠子不同，不要求珠子上的数字有什么别的关系，所以这个题是完全割裂的两个问题：

怎么求珠子的数量？
已知珠子的数量，怎么求合法的项链的数量？

珠子数量

第一问

定义一个珠子为一个无序三元组 $(x,y,z)$ ，满足 $0<x,y,z\leqslant a$ ，且 $\gcd(x,y,z)=1$ 。

给定 $a$ ，求有多少种珠子。

互质这个条件就显得非常数论，然而一般的数论枚举是有序的，因此转化为有序问题：

若 $x\neq y\neq z$ ，设有序三元组数量为 $a$ ，则无序三元组数量为 $\frac{a}{3!}=\frac{a}{6}$ 。
若有两个元素相同，设有序三元组数量为 $a$ ，则无序三元组数量为 $\frac{a}{\operatorname{A}_3^2}=\frac{a}{3}$ 。
若 $x=y=z$ ，因为三个数要互质，显然只有 $(1,1,1)$ 符合条件。

先考虑考虑第二个怎么算：

cnt_2 = -3 + 3\sum_{i=1}^a\sum_{j=1}^a [\gcd(i,j)=1]

莫反人 DNA 已经动了：

-3 + 3\sum_{i=1}^a \sum_{j=1}^a \sum_{d\mid \gcd(i,j)} \mu(d)

然后换 $\sum$ ：

-3 + 3\sum_{d=1}^a \mu(d) \left\lfloor \frac{a}{d} \right\rfloor^2

第一个？同理即可，不多赘述：

cnt_1 = -1 -cnt_2 + \sum_{d=1}^a \mu(d)\left\lfloor \frac{a}{d} \right\rfloor^3

于是这一部分答案为 $\frac{cnt_1}{6} + \frac{cnt_2}{3} + 1$ 。

项链数量

第二问

有 $m$ 种互不相同的珠子，每种无限个。求可以串成多少种本质不同的大小为 $n$ 的项链，要求相邻珠子不同种。

两个项链本质相同，当且仅当其旋转后相同。

“本质不同的项链”非常的群论，我们考虑枚举循环节，设 $f_i$ 表示循环节个数为 $i$ 时，不动点的数量：

ans = \sum_{i=1}^n f(\gcd(i,n)) = \frac{1}{n} \sum_{d\mid n} \varphi(\frac{n}{d})f(d)

考虑 $f_i$ 的现实意义：大小为 $i$ 的合法项链的染色的方案数。考虑如何递推：对于一个大小为 $n-1$ 的合法项链，枚举一个断点，加入一个和两边颜色都不同的珠子；或者对于一个大小为 $n-2$ 的合法项链，先加入一个和一边颜色相同的珠子，然后在两个珠子之间加一个不同的颜色劈开。（因为项链是合法的，因此这两种统计不重）：

f_n = (m-2)f_{n-1} + (m-1)f_{n-2}

边界 $f_1 = 0,f_2 = m(m-1)$ 。看 $n$ 的范围 $10^{14}$ ，线性递推就算了吧。因为本人太菜不会常系数齐次线性递推，只能带大家用生成函数乱搞一下。怎么用生成函数求通项公式？写出 OGF、代入递推式（注意边界系数）、求出封闭形式、换另一种方法展开。

我们写出 $f_i$ 的 OGF：

F(x) = \sum_{i=1}^\infty f_ix^i

根据递推式得到：

F(x) = x(m-2)F(x) + x^2(m-1)F(x) + x^2f_2 + xf_1

整理，将 $f_1 = 0,f_2 = m(m-1)$ 代入：

F(x) = \frac{x^2m(m-1)}{1-(m-2)x-(m-1)x^2}

如何展开？不妨待定系数法，设待定系数 $A,B,a,b$ 满足：

\frac{Ax}{1-ax}+\frac{Bx}{1-bx} = \frac{x^2m(m-1)}{1-(m-2)x-(m-1)x^2}

整理得到：

\frac{(A+B)x - (Ab+aB)x^2}{1-(a+b)x+abx^2} = \frac{x^2m(m-1)}{1-(m-2)x-(m-1)x^2}

于是得到：

\begin{cases} A+B=0 \\ Ab+aB=-m(m-1) \\ a+b = m-2 \\ ab=-(m-1) \end{cases}

可以先解出 $a,b$ ，由第三个式子， $b=m-a-2$ ，代入第四个式子解一元二次方程：

-a^2 - (2-m)a + m-1 = 0

解二次方程都会吧，因为 $m\geqslant 0$ ，容易解得 $a_1=m-1,a_2 = -1$ 。不妨令 $a=m-1,b=-1$ 。

\begin{cases} A + B = 0 \\ -A + (m-1) B = m-m^2 \end{cases}

解二元一次方程都会吧，容易解得 $A=m-1,B=1-m$ 。

F(x) = x\left(\frac{A}{1-ax}+\frac{B}{1-bx} \right) = x\left(\frac{m-1}{1-(m-1)x} + \frac{1-m}{1-(-1)x}\right)

比值为 $p$ 的等比数列的生成函数为 $G(x)=\frac{1}{1-px}$ ，因此拆开化简一下：

\sum_{i\geqslant 1} x^i \left[ (m-1)^i + (-1)^i(m-1) \right]

因此这个序列的通项公式为：

f_n = (m-1)^n + (-1)^n(m-1)

又因为：

ans = \frac{1}{n} \sum_{d\mid n} \varphi(\frac{n}{d})f(d)

于是这个部分也完成了。

细节

细节很多，调了一个上午，看题解才豁然开朗。

写这种题就是要像生产队的驴一样——多磨（模）。
$n$ 可能是模数的倍数，不能直接乘以 $n$ 。可以考虑先在模 $mod^2$ 意义下求解，最后在转化为模 $10^9+7$ 的答案。如果你 WA on #3 读到 $0$ ，就是这个原因。
用上述方法时要么龟速乘要么 __in128，像我这种懒人当然选择 __128 了。

代码

#include<bits/stdc++.h>
#define int __int128
using namespace std;
const int max_n=10000001,MOD=1000000007;int mod;
inline int read(){
    int x=0;bool w=0;char c=getchar();
    while(c<'0' || c>'9') w|=c=='-',c=getchar();
    while(c>='0' && c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return w?-x:x;
}
inline void write(int x){
    if(x<0) putchar('-'),x=-x;
    if(x>9) write(x/10);
    putchar(x%10^48);
}
inline int mi(int a,int p=mod-2){
    int res=1;
    while(p){
        if(p&1) res=res*a%mod;
        a=a*a%mod,p>>=1;
    }
    return res;
}
struct par{
    int x,y;
};
inline par makep(int a,int b){
    par res;
    res.x=a,res.y=b;
    return res;
}

int n,m,iv6,iv62=833333345000000041LL,iv;

signed pr[max_n/10],mu[max_n],cntp;
bool isp[max_n];

inline void Getprime(int n){
    memset(isp,1,sizeof(isp)),isp[1]=0;
    mu[1]=1;
    for(register int i=2;i<=n;++i){
        if(isp[i]){
            pr[++cntp]=i,
            mu[i]=-1;
        }
        for(register int j=1;j<=cntp && i*pr[j]<=n;++j){
            isp[i*pr[j]]=0;
            if(i%pr[j]==0){
                mu[i*pr[j]]=0;
                break;
            }
            mu[i*pr[j]]=-mu[i];
        }
    }
    for(register int i=1;i<=n;++i)
        mu[i]=mu[i-1]+mu[i];
}

inline int F(int n){
    return (mi(m-1,n)+(n&1?mod-1:1)*(m-1)%mod)%mod;
}

int cnt1,cnt2;

inline int sub1(int a){
    cnt1=cnt2=0;
    for(register int l=1,r,k;l<=a;l=r+1){
        r=(a/(a/l)),k=a/l;
        cnt2=(cnt2+(int)(mu[r]-mu[l-1])*k*k)%mod,
        cnt1=(cnt1+(int)(mu[r]-mu[l-1])*k*k%mod*k)%mod;
    }
    cnt1=(cnt1%mod+mod)%mod,
    cnt2=(cnt2%mod+mod)%mod;

    cnt2=(mod-3+3*cnt2)%mod,
    cnt1=(mod-1-cnt2+cnt1)%mod;
    return (iv6*cnt1+iv6*2*cnt2+1)%mod;
}

int lim,ans;
vector<par> pfac; 

inline void fact(int n){
    pfac.resize(0);
    for(register int i=1;i<=cntp && pr[i]*pr[i]<=n;++i)
        if(n%pr[i]==0){
            int k=0;
            while(n%pr[i]==0)
                n/=pr[i],++k;
            pfac.push_back(makep(pr[i],k));
        }
    if(n>1)
        pfac.push_back(makep(n,1));
    lim=pfac.size();
}

inline void dfs(int l,int a,int phi){
    if(l==lim){
        ans=(ans+phi*F(n/a))%mod;
        return;
    }
    dfs(l+1,a,phi);
    for(register int i=1;i<=pfac[l].y;++i){
        a=a*pfac[l].x;
        phi=phi*(i==1?pfac[l].x-1:pfac[l].x)%mod;
        dfs(l+1,a,phi);
    }
}

inline int sub2(){
    ans=0,fact(n);
    dfs(0,1,1);
    return ans%mod;
}

signed main(){
    mod=MOD,Getprime(10000000),iv=iv6=mi(6);
    for(register int T=read();T>0;--T){
        n=read();bool ok=0;
        if(n%MOD==0) mod=MOD*MOD,iv6=iv62;
        else mod=MOD,ok=1,iv6=iv;
        m=sub1(read());
        int ans=sub2();
        if(!ok) write(ans/MOD*mi(n/MOD)%MOD);
        else write(ans*mi(n%MOD)%MOD);
        putchar('\n');
    }
    return 0;
}