【JZOJ3298】【SDOI2013】项链 莫比乌斯反演+Polya计数法+欧拉函数+通项公式

任务

$对于 100\%的 数据：所有的n<=10^14 ，a<=10^7，T<=10；$

解法

显然我们可以分两步走：
1.计算珠子的种类数
2.有多少条本质不同的项链

---

Part Ⅰ

我们计算珠子的种类数，
显然是一个莫比乌斯反演的简单应用。
利用

$$Ans'=\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor^3*\mu(i)$$
我们可以求出有序三元组的个数。
但是我们要求的是无序三元组的个数。
显然我们会有算重的部分。
考虑：
$1\ 2\ 3$型会被算 $6$次，
$1\ 2\ 2$型会被算 $6$次，实际只需计算$3$次，需补 $3$次；
$1\ 1\ 1$型同样会被算 $6$次，实际只需计算$1$次，需补 $5$次；

所以我们需要把多算的加回进去。
容易算出有序二元组的个数，

$$Ans''=\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor^2*\mu(i)$$
所以 $$m=\frac{1}{6}(Ans'+Ans''*3+2)$$
（由于 $1\ 1\ 1$只有1 1 1这种情况，又被二元组补了 $3$次，所以只需补$2$次）

这样我们就求出了珠子的种类数$m$。

Part Ⅱ

接下来我们求有多少个本质不同的项链；
对于这种求本质不同的计数问题，我们通常会想到Burnside引理+Polya计数法
由Polya计数法，我们有
对于一个旋转置换$i$，我们有它的循环个数为$gcd(i,n)$；
又有，这些循环的起始位置都是相邻的；
为了满足题目的“相邻珠子必须不同”的条件。
必须求出一个$f(k)$，表示$k$个满足相邻不同的珠子的环的方案数。

---

设$g(k)$表示$k$个珠子相邻不同，但首尾必须相同的方案数。
显然

$$g(i)=f(i-2)$$
又有，

$$f(i)=f(i-1)*(m-2)+g(i-1)*(m-1)$$
则，
$$f(i)=f(i-1)*(m-2)+f(i-2)*(m-1)$$
这个递推式我们可以利用 矩阵乘法来做，
但由于题目的苛刻时限，使得我们需要进一步探寻通项公式。

---

对于这一类，$f(i)=p*f(i-1)+q*f(i-2)$的递推式，
求通项公式都有惯用的套路，以本题为例：

$$\begin{align} f(i)&=f(i-1)*(m-2)+f(i-2)*(m-1)\\ f(i)+\lambda*f(i-1)&=f(i-1)*(m-2+\lambda)+f(i-2)*(m-1)\\ 令m-1=(m-2+\lambda)*\lambda，则\lambda=1或m-1\\ 将\lambda=1代入得，\\ f(i)+f(i-1)&=(m-1)*(f(i-1)+f(i-2))\\ f(i)&=(m-1)^{i-3}*(f(1)+f(2))-f(i-1) \end{align}$$
其中 $f(1)=0,f(2)=m(m-1)$；
显然我们可以利用等比数列继续推导，我们可以得到通项公式 $(i>2)$：
令 $h=m-1$，
$$f(i)=\begin{cases} h^{\frac{h^{i-1}+h}{h+1}*f(2)}-f(2),&i是奇数\\ h^{\frac{h^{i-1}-h}{h+1}*f(2)}+f(2),&i是偶数\\ \end{cases}$$
这样，我们只需使用快速幂就能求 $f$。

---

于是，我们有

$$ans=\frac{1}{n}*\sum_{i=1}^nf(gcd(n,i))$$
很显然，这是超时的。
考虑有多少个 $j$满足，$gcd(n,j)=i$。
$gcd(n,j)=i\Rightarrow gcd(\frac{n}{i},\frac{j}{i})=1$
即 $\varphi(\frac{n}{i})$个。
考虑更换主体，则有
$$ans=\frac{1}{n}*\sum_{i|n}f(i)*\varphi(\frac{n}{i})$$
于是就完成了。
另外在， $mo|n$的情况下，我们需要特殊处理，因为费马小定理失效了。

代码

#include<iostream>
#include<algorithm>
#include<stdio.h>
#include<math.h>
#include<string.h>
#define ll long long
using namespace std;
const char* fin="ex3298.in";
const char* fout="ex3298.out";
const ll inf=0x7fffffff;
const ll maxn=10000007;
ll t,n,m,i,j,k,ans3,ans2,ans,mo,h,f2,ni;
ll mu[maxn],p[maxn],phi[maxn],c[maxn],d[maxn];
bool bz[maxn];
void pre(){
    ll i,j,k;
    mu[1]=1;
    for (i=2;i<maxn;i++){
        if (!bz[i]){
            mu[i]=-1;
            p[++p[0]]=i;
        }
        for (j=1;j<=p[0];j++){
            k=p[j]*i;
            if (k>=maxn) break;
            bz[k]=true;
            if (i%p[j]==0){
                mu[k]=0;
                break;
            }else mu[k]=-mu[i];
        }
    }
}
ll qmul(ll a,ll b){
    ll c=0,v=1;
    if (a<0 && b>0) a=-a,v=-1;
    else if (b<0 && a>0) v=-1,b=-b;
    while (b){
        if (b&1) c=(a+c)%mo;
        a=(a+a)%mo;
        b>>=1;
    }
    return c*v;
}
ll qpower(ll a,ll b){
    ll c=1;
    while (b){
        if (b&1) c=qmul(c,a);
        a=qmul(a,a);
        b>>=1;
    }
    return c;
}
ll N(ll a){
    return qpower(a,mo-2);
}
void mobius(){
    ll i,j,k;
    ans3=0;
    ans2=0;
    for (i=1;i<=m;i++){
        ans3=(ans3+qmul(qmul((m/i),(m/i)),(m/i))*mu[i])%mo;
        ans2=(ans2+qmul((m/i),(m/i))*mu[i])%mo;
    }
}
ll f(ll v){
    if (v==1) return 0;
    if (v==2) return f2;
    if (v&1) return qmul(qmul(qpower(h,v-1)+h,ni),f2)-f2;
    else return qmul(qmul(qpower(h,v-1)-h,ni),f2)+f2;
}
void dfs(ll v,ll phi,ll l,bool bz,bool yi){
    if (bz) ans=(ans+qmul(phi,f(n/v)))%mo;
    if (d[l]){
        d[l]--;
        dfs(v*c[l],(yi?phi*c[l]:phi*(c[l]-1)),l,true,true);
        d[l]++;
    }
    if (l<c[0]) dfs(v,phi,l+1,false,false);
}
int main(){
    scanf("%lld",&t);
    ll _mo=1e9+7;
    pre();
    while (t--){
        scanf("%lld%lld",&n,&m);
        if (n%_mo==0) mo=_mo*_mo;
        else mo=_mo;
        mobius();
        m=qmul(N(6),ans3+ans2*3+2);
        m=(m+mo)%mo;
        j=(ll)sqrt(n);
        h=m-1,f2=qmul(m,m-1),ni=N(h+1);
        c[0]=0;
        ll nn=n;
        for (i=2;i<=j;i++){
            if (nn%i==0){
                c[++c[0]]=i;
                d[c[0]]=0;
                while (nn%i==0){
                    d[c[0]]++;
                    nn/=i;
                }
            }
        }
        if (nn>1){
            c[++c[0]]=nn;
            d[c[0]]=1;
        }
        ans=0;
        dfs(1,1,1,true,false);
        if (mo==_mo) ans=ans*N(n)%mo;
        else{
            ans/=_mo;
            ans=ans*N(n/_mo)%_mo;
        }
        ans=(ans%_mo+_mo)%_mo;
        printf("%lld\n",ans);
    }
    return 0;
}