题目链接

https://www.luogu.org/problemnew/show/P4916

题解

将项链用序列表示,11代表黑色,00代表白色,对于一个合法序列,它必定能表示成dd个循环节,每个循环节nd\frac{n}{d}个珠子,其中md\frac{m}{d}个黑色珠子。假设循环节长度nn的合法序列方案数为f(n)f(n)(不考虑旋转后相同的情况),容易发现答案就是
dn1df(d) \sum_{d|n}\frac{1}{d}f(d)
定义
g(n)=dnf(d) g(n)=\sum_{d|n}f(d)
容易发现
f(n)=dnμ(nd)g(d) f(n)=\sum_{d|n}\mu(\frac{n}{d})g(d)
考虑如何求g(d)g(d)。显然可以对每个循环节分开考虑,对于长度为aa,有bb颗黑珠子的循环节,考虑被aba-b颗白珠子划分成的黑珠子每一段的数量,g(a)g(a)就是下面方程的整数解的方案数
i=0abxi=b(i[0,ab],0xik,x0+xabk) \sum_{i=0}^{a-b}x_i=b(\forall i\in [0,a-b],0\leq x_i\leq k,x_0+x_{a-b}\leq k)
容易发现g(a)g(a)就是下面式子在xbx^b项的系数
(i=0kxi)ab1((i=0kxi)2 mod xk+1) (\sum_{i=0}^{k}x^i)^{a-b-1}((\sum_{i=0}^{k}x^i)^2\bmod x^{k+1})
展开右边
(i=0kxi)ab1(i=0k(i+1)xi) (\sum_{i=0}^k x^i)^{a-b-1}(\sum_{i=0}^{k}(i+1)x^i)

(1xk+11x)ab1(1(k+2)xk+1+(k+1)xk+2(1x)2) (\frac{1-x^{k+1}}{1-x})^{a-b-1}(\frac{1-(k+2)x^{k+1}+(k+1)x^{k+2}}{(1-x)^2})
去括号
(1xk+1)ab1(1x)(ab+1)(1(k+2)xk+1+(k+1)xk+2) (1-x^{k+1})^{a-b-1}(1-x)^{-(a-b+1)}(1-(k+2)x^{k+1}+(k+1)x^{k+2})
利用二项式定理
(i=0(ab1i)(1)ix(k+1)i)(i=0(ab+ii)xi)(1(k+2)xk+1+(k+1)xk+2) (\sum_{i=0}^{\infin}\binom{a-b-1}{i}(-1)^ix^{(k+1)i})(\sum_{i=0}^{\infin}\binom{a-b+i}{i}x^i)(1-(k+2)x^{k+1}+(k+1)x^{k+2})
定义
S(n)=(k+1)i+j=n(ab1i)(1)i(ab+ii) S(n)=\sum_{(k+1)i+j=n}\binom{a-b-1}{i}(-1)^i\binom{a-b+i}{i}
则有
g(a)=S(b)(k+2)S(bk1)+(k+1)S(bk2) g(a)=S(b)-(k+2)S(b-k-1)+(k+1)S(b-k-2)
S(n)S(n)可以枚举ii,时间复杂度就是O(σ1(n)k+1)O(\frac{\sigma_1(n)}{k+1}),由于σ1(n)\sigma_1(n)不会很大,可以通过此题。

代码

#include <cstdio>

int read()
{
  int x=0,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=100000;
const int mod=998244353;

int p[maxn+10],prime[maxn+10],cnt,mu[maxn+10],fac[maxn+10],inv[maxn+10],ifac[maxn+10];

int getprime()
{
  p[1]=mu[1]=1;
  for(int i=2; i<=maxn; ++i)
    {
      if(!p[i])
        {
          prime[++cnt]=i;
          mu[i]=-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]=-mu[i];
        }
    }
  fac[0]=1;
  for(int i=1; i<=maxn; ++i)
    {
      fac[i]=1ll*fac[i-1]*i%mod;
    }
  inv[0]=inv[1]=1;
  for(int i=2; i<=maxn; ++i)
    {
      inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    }
  ifac[0]=1;
  for(int i=1; i<=maxn; ++i)
    {
      ifac[i]=1ll*ifac[i-1]*inv[i]%mod;
    }
  return 0;
}

inline int C(int a,int b)
{
  if((a<b)||(b<0)||(a<0))
    {
      return 0;
    }
  return 1ll*fac[a]*ifac[b]%mod*ifac[a-b]%mod;
}

inline int getsum(int a,int b,int k,int mx)
{
  int ans=0;
  for(int i=0; mx-(k+1)*i>=0; ++i)
    {
      int j=mx-(k+1)*i;
      ans=(ans+((i&1)?(-1ll):(1ll))*C(a-b-1,i)*C(a-b+j,j))%mod;
      if(ans<0)
        {
          ans+=mod;
        }
    }
  return ans;
}

inline int count(int a,int b,int k)
{
  int ans=(getsum(a,b,k,b)-1ll*(k+2)*getsum(a,b,k,b-k-1)+1ll*(k+1)*getsum(a,b,k,b-k-2))%mod;
  if(ans<0)
    {
      ans+=mod;
    }
  return ans;
}

int n,m,k,f[maxn+10],g[maxn+10];

int getG(int d)
{
  g[n/d]=count(n/d,m/d,k);
  return 0;
}

int main()
{
  getprime();
  n=read();
  m=read();
  k=read();
  if(m==0)
    {
      puts("1");
      return 0;
    }
  for(int i=1; i*i<=m; ++i)
    {
      if(m%i==0)
        {
          if(n%i==0)
            {
              getG(i);
            }
          int j=m/i;
          if((i*i!=m)&&(n%j==0))
            {
              getG(j);
            }
        }
    }
  for(int i=1; i<=n; ++i)
    {
      for(int j=i; j<=n; j+=i)
        {
          f[j]+=mu[j/i]*g[i];
          if(f[j]>=mod)
            {
              f[j]-=mod;
            }
          if(f[j]<0)
            {
              f[j]+=mod;
            }
        }
    }
  int ans=0;
  for(int i=1; i<=n; ++i)
    {
      if(n%i==0)
        {
          ans=(ans+1ll*inv[i]*f[i])%mod;
        }
    }
  printf("%d\n",ans);
  return 0;
}