BZOJ 4815 [Cqoi2017]小Q的表格

题目链接

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

题解

观察发现$bf(a,a+b)=(a+b)f(a,b)$很像更相减损术的式子，稍加推导可得
$$f(a,b)=\frac{a}{\gcd (a,b)}\frac{b}{\gcd (a,b)}f(\gcd (a,b),\gcd (a,b))$$
容易发现答案就是
$$\sum _{g=1}^{k}f(g,g)\sum _{i=1}^{\lfloor k/g\rfloor }\sum _{j=1}^{\lfloor k/g\rfloor }ij[\gcd (i,j)=1]$$
反演可得答案为
$$\sum _{g=1}^{k}f(g,g)\sum _{i=1}^{\lfloor k/g\rfloor }{i}^2\phi (i)$$
后面一部分可以预处理，那么可以在$O(\sqrt{k})$时间内处理询问。

注意修改操作是$O(1)$的，因此需要设计一种$O(\sqrt{n})$修改，$O(1)$查询的数据结构维护$f(g,g)$。

代码

#include <cmath>
#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=4000000;
const int mod=1000000007;
 
int p[maxn+10],prime[maxn+10],cnt,phi[maxn+10],inv[maxn+10];
 
int getprime()
{
  p[1]=phi[1]=1;
  for(int i=2; i<=maxn; ++i)
    {
      if(!p[i])
        {
          prime[++cnt]=i;
          phi[i]=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)
            {
              phi[x]=phi[i]*prime[j];
              break;
            }
          phi[x]=phi[i]*(prime[j]-1);
        }
    }
  for(int i=1; i<=maxn; ++i)
    {
      phi[i]=(phi[i-1]+1ll*phi[i]*i%mod*i)%mod;
    }
  inv[0]=inv[1]=1;
  for(int i=2; i<=maxn; ++i)
    {
      inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
    }
  return 0;
}
 
int block[maxn+10],sum[maxn+10],bs[maxn+10],val[maxn+10],totb;
 
int modify(int pos,int v)
{
  int k=v-val[pos];
  if(k<0)
    {
      k+=mod;
    }
  val[pos]+=k;
  if(val[pos]>=mod)
    {
      val[pos]-=mod;
    }
  for(int i=pos; block[i]==block[pos]; ++i)
    {
      sum[i]+=k;
      if(sum[i]>=mod)
        {
          sum[i]-=mod;
        }
    }
  for(int i=block[pos]; i<=totb; ++i)
    {
      bs[i]+=k;
      if(bs[i]>=mod)
        {
          bs[i]-=mod;
        }
    }
  return 0;
}
 
int getsum(int l,int r)
{
  (l==0)&&(++l);
  int v=bs[block[r]-1]-bs[block[l]-1]+sum[r];
  if(v<0)
    {
      v+=mod;
    }
  if(v>=mod)
    {
      v-=mod;
    }
  if(block[l-1]==block[l])
    {
      v-=sum[l-1];
      if(v<0)
        {
          v+=mod;
        }
    }
  return v;
}
 
int gcd(int x,int y)
{
  return y?gcd(y,x%y):x;
}
 
int quickpow(int a,int b)
{
  int res=1;
  while(b)
    {
      if(b&1)
        {
          res=1ll*res*a%mod;
        }
      a=1ll*a*a%mod;
      b>>=1;
    }
  return res;
}
 
int m,n,a,b,k;
long long x;
 
int main()
{
  getprime();
  m=read<int>();
  n=read<int>();
  totb=sqrt(n);
  for(int i=1; i<=n; ++i)
    {
      block[i]=(i-1)/totb+1;
    }
  totb=(n-1)/totb+1;
  for(int i=1; i<=n; ++i)
    {
      val[i]=1ll*i*i%mod;
      bs[block[i]]+=val[i];
      if(bs[block[i]]>=mod)
        {
          bs[block[i]]-=mod;
        }
    }
  for(int i=1; i<=totb; ++i)
    {
      bs[i]+=bs[i-1];
      if(bs[i]>=mod)
        {
          bs[i]-=mod;
        }
    }
  for(int i=1; i<=n; ++i)
    {
      sum[i]=((block[i]==block[i-1])?sum[i-1]:0)+val[i];
      if(sum[i]>=mod)
        {
          sum[i]-=mod;
        }
    }
  for(int i=1; i<=m; ++i)
    {
      a=read<int>();
      b=read<int>();
      x=read<long long>();
      k=read<int>();
      x%=mod;
      int g=gcd(a,b);
      modify(g,x*quickpow((1ll*(a/g)*(b/g))%mod,mod-2)%mod);
      int ans=0;
      for(int l=1,r; l<=k; l=r+1)
        {
          r=k/(k/l);
          ans=(ans+1ll*getsum(l,r)*phi[k/l])%mod;
        }
      printf("%d\n",ans);
    }
  return 0;
}