BZOJ 3202 [Sdoi2013]项链

题目链接

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

题解

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

对于珠子的个数，记$S3$为$\gcd=1$的三个数的个数，$S2$为$\gcd=1$的两个数的个数，$S1$为$\gcd=1$的一个数的个数，容易发现答案就是
$S1+\frac{1}{3}(3S2-3S1)+\frac{1}{6}(S3-3S2+2S1)=\frac{S3+3S2+2S1}{6}$
其中$S1=1$，$S2,S3$可以用莫比乌斯反演求得
$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)$
对于项链的个数，假设珠子的个数为$m$，长度为$i$的项链个数（不考虑旋转）为$f(i)$，由polya原理可得答案就是
$\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-1)$转移，此时当前珠子不能与两端相同；要么从$f(i-2)$转移，插入一个和第一个珠子相同的，再插入一个不同的，因此$f(i)$的递推式为
$f(i)=(m-2)f(i-1)+(m-1)f(i-2)$
用特征方程解得
$f(i)=(-1)^i(m-1)+(m-1)^i$
由于$n$有可能是$10^9+7$的倍数，此时没有逆元，因此可以先算模$(10^9+7)^2$的答案，如果$n$是$10^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;
}