洛谷 P3307 [SDOI2013]项链 burnside引理+polya定理+莫比乌斯反演

题目描述
项链是人体的装饰品之一，是最早出现的首饰。项链除了具有装饰功能之外，有些项 链还具有特殊显示作用，如天主教徒的十字架链和佛教徒的念珠。

从古至今人们为了美化人体本身，也美 化环境，制造了各种不同风格，不同特点、不同式样的项链，满足了不同肤色、不同民族、不同审美观的人的审美需要。就材料而论，首饰市场上的项链有黄金、白银、珠宝等几种。

珍珠项链为珍珠制成的饰品，即将珍珠 钻孔后用线串在一起，佩戴于项间。天然珍珠项链具有一定的护养作用。 最近，铭铭迷恋上了一种项链。与其他珍珠项链基本上相同，不过这种项链的珠子却 与众不同，是正三菱柱的泰山石雕刻而成的。

三菱柱的侧面是正方形构成的，上面刻有数字。 能够让铭铭满意的项链必须满足下面的条件：

1：这串项链由n颗珠子构成的。

2：每一个珠子上面的数字x,必须满足0<x<=a，且珠子上面的数字的最大公约数要恰 好为1。两个珠子被认为是相同的，当且仅当他们经过旋转，或者翻转后能够变成一样的。

3：相邻的两个珠子必须不同。

4：两串项链如果能够经过旋转变成一样的，那么这两串项链就是相同的！ 铭铭很好奇如果给定n和a，能够找到多少不同串项链。由于答案可能很大，所以对输 出的答案mod 1000000007。

输入输出格式

输入格式：
数据由多组数据构成： 第一行给定一个T<=10，代表由T组数据。 接下来T行，每行两个数n和a。

输出格式：
对于每组数据输出有多少不同的串。

输入输出样例

输入样例#1：
1
2 2
输出样例#1：
3
说明

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

样例解释：由三种珠子：[1,1,1],[1,1,2],[1,2,2].组成的串有：[1,2],[1,3],[2,3]。

分析：
我们先求出不同的珠子有多少个。相当于求有序三元组$(x,y,z)$使得$gcd(x,y,z)=1$。
我们考虑化成无序三元组。对于任意一个$(x,y,z)$无序时会被算6次。而当有$x=y$时，少算三次。当$x=y=z$时，此时$x=y=z=1$，这个少算5次，但是有3次已经在$x=y$算过，所欲只要加2次。
即，设$s_1$为无序三元组答案，$s_2$为无序二元组方案，那么不同珠子数$m=\frac{1}{6}\ast (s_1+3\ast s_2+2)$。

考虑求$s_1$和$s_2$，很套路的莫比乌斯反演一下就得到
$$s_1=\sum _{i=1}^a\lfloor \frac{a}{i}{\rfloor }^3\ast \mu (i)$$
$$s_2=\sum _{i=1}^a\lfloor \frac{a}{i}{\rfloor }^2\ast \mu (i)$$

考虑能拼出多少个项链。
可以考虑用burnside引理+polya定理。显然置换有$n$个。分别是$(1,2,3,....n)$，$(2,3,4,...,n,1),$,…
考虑循环，假设第$i$个置换是$1$对$i$的那个，显然循环有$gcd(n,i)$个，又因为这些循环是相连的，而循环中的取值相同。
我们设$f(x)$为一个大小为$x$的环，相邻两个颜色都不同的方案。那么第$i$个置换的答案就是$f(gcd(n,i))$。
也就是$ans={\sum }_{i=1}^{n}f(gcd(n,i))$。
我们设$d=gcd(n,i)$，则$ans={\sum }_{d|n}f(d)\ast {\sum }_{i=1}^n[gcd(n,i)=d]$
然后也就是$ans={\sum }_{d|n}f(d)\ast \varphi (n/d)$。
而$f(x)=f(x-1)\ast (m-2)+f(x-2)\ast (m-1)$
可以理解为在$x-1$个环的任意一个位置，由于两边颜色不同，可以取$m-2$种；或者在$x-2$个环中放一个颜色相同的，再在两个相同颜色中放一个颜色不同的，有$m-1$种情况。可以矩阵乘搞。

还有就是$n$有可能是$mod$的倍数，所以上述运算要在$mo{d}^2$意义下跑。如果$n$不是倍数，那么直接把算出来的答案模$mod$再乘逆元；否则把结果先除n，然后乘$\frac{n}{p}$在$mod$意义下的值。

代码：

// luogu-judger-enable-o2
#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long

const LL p=1e9+7;
const LL mod=p*p;
const int maxn=1e7+7;
const LL inv6=833333345000000041;

using namespace std;

int T;
LL n,m,cnt,k,ans,num;
LL prime[maxn],mu[maxn],q[10007],a[10007];
bool not_prime[maxn];

LL mul(LL x,LL y)
{
    LL tmp=x*y-(LL)((long double)x*y/mod+0.1)*mod;
    if (tmp<0) tmp+=mod;
    return tmp;
}

void getmu(LL n)
{
    mu[1]=1;
    for (LL i=2;i<=n;i++)
    {
        if (!not_prime[i])
        {
            prime[++cnt]=i;
            mu[i]=-1;
        }
        for (LL j=1;j<=cnt;j++)
        {
            if (i*prime[j]>n) break;
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for (LL i=1;i<=n;i++) mu[i]+=mu[i-1];
}

LL power(LL x,LL y)
{
    if (!y) return 1;
    LL c=power(x,y/2);
    c=mul(c,c);
    if (y%2) c=mul(c,x);
    return c;
}

LL power1(LL x,LL y)
{
    if (!y) return 1;
    LL c=power1(x,y/2);
    c=(c*c)%p;
    if (y%2) c=(c*x)%p;
    return c;
}

void calc()
{
    LL s1=0,s2=0;
    for (int i=1,last;i<=m;i=last+1)
    {
        last=m/(m/i);
        s1=(s1+mul(power(m/i,3),(mu[last]+mod-mu[i-1])%mod))%mod;
        s2=(s2+mul(power(m/i,2),(mu[last]+mod-mu[i-1])%mod))%mod;
    }
    k=mul(inv6,((s1+mul(s2,3)+2)%mod));
}

LL F(LL x)
{
    LL tmp=power(k-1,x);
    if (x&1) tmp=(tmp+1+mod-k)%mod;
        else tmp=(tmp+k-1)%mod;
    return tmp;
}

void dfs(LL x,LL d,LL phi)
{
    if (x>num)
    {
        ans=(ans+mul(phi,F(n/d)))%mod;
        return;
    }
    dfs(x+1,d,phi);
    int s=1;
    for (int i=1;i<=q[x];i++)
    {
        d*=a[x];
        phi*=a[x]-s;
        s=0;
        dfs(x+1,d,phi);
    }
}

int main()
{
    scanf("%d",&T);
    getmu(1e7);		
    while (T--)
    {
        scanf("%lld%lld",&n,&m);
        calc();
        LL x=n;
        num=0;		
        for (int i=1;i<=cnt&&(prime[i]*prime[i]<=n);i++)
        {
            if (x%prime[i]==0)
            {
                a[++num]=prime[i];
                q[num]=0;
                while (x%prime[i]==0) q[num]++,x/=prime[i];
            }
        }
        if (x) a[++num]=x,q[num]=1;
        ans=0;
        dfs(1,1,1);
        if (n%p==0) ans=(ans/p)*power1(n/p,p-2)%p;
               else ans=(ans%p)*power1(n%p,p-2)%p;
        printf("%lld\n",ans);
    }
}