题目链接

https://lydsy.com/JudgeOnline/problem.php?id=3202

题解

这题可以分成两部分:一个是统计珠子的个数,一个是统计项链的个数。

对于珠子的个数,记S3S3gcd=1\gcd=1的三个数的个数,S2S2gcd=1\gcd=1的两个数的个数,S1S1gcd=1\gcd=1的一个数的个数,容易发现答案就是
S1+13(3S23S1)+16(S33S2+2S1)=S3+3S2+2S16 S1+\frac{1}{3}(3S2-3S1)+\frac{1}{6}(S3-3S2+2S1)=\frac{S3+3S2+2S1}{6}
其中S1=1S1=1S2,S3S2,S3可以用莫比乌斯反演求得
S2=T=1nnT2μ(T),S3=T=1nnT3μ(T) S2=\sum_{T=1}^n \lfloor\frac{n}{T}\rfloor^2\mu(T),S3=\sum_{T=1}^n \lfloor\frac{n}{T}\rfloor^3\mu(T)
对于项链的个数,假设珠子的个数为mm,长度为ii的项链个数(不考虑旋转)为f(i)f(i),由polya原理可得答案就是
1ni=1nf(gcd(i,n))=1ndnf(d)φ(nd) \begin{aligned} & \frac{1}{n}\sum_{i=1}^n f(\gcd(i,n))\\ = & \frac{1}{n}\sum_{d|n}f(d)\varphi(\frac{n}{d}) \end{aligned}
由于限制了相邻珠子不能相同,因此可以考虑f(i)f(i)的递推式:要么从f(i1)f(i-1)转移,此时当前珠子不能与两端相同;要么从f(i2)f(i-2)转移,插入一个和第一个珠子相同的,再插入一个不同的,因此f(i)f(i)的递推式为
f(i)=(m2)f(i1)+(m1)f(i2) f(i)=(m-2)f(i-1)+(m-1)f(i-2)
用特征方程解得
f(i)=(1)i(m1)+(m1)i f(i)=(-1)^i(m-1)+(m-1)^i
由于nn有可能是109+710^9+7的倍数,此时没有逆元,因此可以先算模(109+7)2(10^9+7)^2的答案,如果nn109+710^9+7的倍数直接除掉就可以了。

代码

#include <cstdio>
#include <algorithm>

template<typename T>
T read()
{
  T x=0;
  int f=1;
  char ch=getchar();
  while((ch<'0')||(ch>'9'))
    {
      if(ch=='-')
        {
          f=-f;
        }
      ch=getchar();
    }
  while((ch>='0')&&(ch<='9'))
    {
      x=x*10+ch-'0';
      ch=getchar();
    }
  return x*f;
}

const int maxn=10000000;
const int mo=1000000007;
const long long mod=1ll*mo*mo;
const long long inv6=833333345000000041ll;

long long mul(long long x,long long y)
{
  return (x*y-(long long)(((long double)x*y+0.5)/(long double)mod)*mod+mod)%mod;
}

int p[maxn+10],prime[maxn+10],cnt;
long long mu[maxn+10];

int getprime()
{
  p[1]=mu[1]=1;
  for(int i=2; i<=maxn; ++i)
    {
      if(!p[i])
        {
          prime[++cnt]=i;
          mu[i]=mod-1;
        }
      for(int j=1; (j<=cnt)&&(i*prime[j]<=maxn); ++j)
        {
          int x=i*prime[j];
          p[x]=1;
          if(i%prime[j]==0)
            {
              mu[x]=0;
              break;
            }
          mu[x]=mod-mu[i];
        }
    }
  for(int i=2; i<=maxn; ++i)
    {
      mu[i]+=mu[i-1];
      if(mu[i]>=mod)
        {
          mu[i]-=mod;
        }
    }
  return 0;
}

long long quickpow(long long a,long long b)
{
  long long res=1;
  while(b)
    {
      if(b&1)
        {
          res=mul(res,a);
        }
      a=mul(a,a);
      b>>=1;
    }
  return res;
}

long long inv(long long x)
{
  return quickpow(x,1ll*mo*(mo-1)-1);
}

long long bead(int a)
{
  long long f=0,g=0;
  for(int l=1,r; l<=a; l=r+1)
    {
      r=a/(a/l);
      f+=mul(mul(a/l,a/l),mu[r]-mu[l-1]+mod);
      if(f>=mod)
        {
          f-=mod;
        }
      g+=mul(mul(a/l,a/l),mul(a/l,mu[r]-mu[l-1]+mod));
      if(g>=mod)
        {
          g-=mod;
        }
    }
  return mul(inv6,g+mul(3,f)+2);
}

long long phi(long long x)
{
  long long res=x;
  if(res%mo==0)
    {
      res=(res/mo)*(mo-1);
      x/=mo;
    }
  for(int i=2; 1ll*i*i<=x; ++i)
    {
      if(x%i==0)
        {
          res=mul(mul(res,i-1),inv(i));
        }
      while(x%i==0)
        {
          x/=i;
        }
    }
  if(x!=1)
    {
      res=mul(mul(res,x-1),inv(x));
    }
  return res;
}

long long F(long long n,long long m)
{
  long long k=quickpow(m-1,n)+mul(quickpow((n&1)?(mod-1):1,n),m-1);
  if(k>=mod)
    {
      k-=mod;
    }
  return k;
}

long long necklace(long long n,long long m)
{
  long long ans=0;
  for(int i=1; 1ll*i*i<=n; ++i)
    {
      if(n%i==0)
        {
          ans+=mul(F(i,m),phi(n/i));
          if(ans>=mod)
            {
              ans-=mod;
            }
          if(1ll*i*i!=n)
            {
              ans+=mul(F(n/i,m),phi(i));
              if(ans>=mod)
                {
                  ans-=mod;
                }
            }
        }
    }
  return ans;
}

int T,m;
long long n;

int main()
{
  getprime();
  T=read<int>();
  while(T--)
    {
      n=read<long long>();
      m=read<int>();
      long long ans=necklace(n,bead(m));
      printf("%lld\n",(ans%mo)?mul(ans,inv(n))%mo:mul((ans/mo),inv(n/mo))%mo);
    }
  return 0;
}