题解 P5572 【[CmdOI2019]简单的数论题】

command_block · 2019-10-02 18:27:16 · 题解

官方题解 : 简单的数论题

打出 \varphi 的表，然后按照题意大力计算即可。

• \Large 40'

打个答案的表就好了, \varphi 线性筛,复杂度 O(n^2\log n+T) 或者 O(n^2+T) 。

#include<algorithm>
#include<cstdio>
#define MaxN 4000500
#define MaxNum 2010
#define mod 23333
using namespace std;
int T,N,M;
bool e[MaxN+500];
int tn,p[MaxN+500],phi[MaxN+500];
int gcd[MaxNum+50][MaxNum+50];
int ans[MaxNum+50][MaxNum+50];
int gcdpre(int a,int b)
{
  if (!b)return a;
  if (gcd[a][b])return gcd[a][b];
  return gcd[a][b]=gcdpre(b,a%b); 
}
inline void getsth()
{
  e[1]=1;phi[1]=1;
  for (int i=2,t;i<=MaxN;i++){
    if (!e[i]){
      p[++tn]=i;
      phi[i]=i-1;
    }
    for (int j=1;j<=tn&&(t=p[j]*i)<=MaxN;j++){
      phi[t]=phi[i]*(i%p[j]==0 ? p[j] : p[j]-1);
      e[t]=1;
      if (i%p[j]==0)break;
    }
  }
  for (int i=1;i<=MaxNum;i++)
    for (int j=1;j<=MaxNum;j++)
      gcd[i][j]=gcdpre(i,j);
  for (int i=1;i<=MaxNum;i++)
    for (int j=1;j<=MaxNum;j++)
      ans[i][j]=(phi[i*j/gcd[i][j]/gcd[i][j]]
        +ans[i-1][j]+ans[i][j-1]-ans[i-1][j-1])%mod;
}
int main()
{
  scanf("%d",&T);
  getsth();
  while(T--){
    scanf("%d%d",&N,&M);
    printf("%d\n",(ans[N][M]+mod)%mod);
  }return 0;
}

• \Large 60'

送分结束,开始推导式子:

求 \sum\limits_{i=1}^n\sum\limits_{j=1}^m\varphi\left(\dfrac{lcm(i,j)}{gcd(i,j)}\right)

• 设 G(n,k)=\sum\limits_{i=1}^n\varphi(ik)

.=\sum\limits_{d=1}^n\sum\limits_{t=1}^{n/d}\mu(t)G(\lfloor n/dt\rfloor,t)G(\lfloor m/dt\rfloor,t)

设 L=\max(n,m)。

不难发现，我们调用的 G(n,k) 函数都满足 n,k\leq L，所以总状态量是 O\Big(\sum\limits_{x,y}[xy\leq n]\Big)=O(L\log L) 的。

不过，由于 t 的存在，上面那个式子仍然不能直接整除分块。

• 设 S(n,m,k)=\mu(k)G(n,k)G(m,k)

则 {\rm Ans}=\sum\limits_{d=1}^n\sum\limits_{t=1}^{n/d}S(\lfloor n/dt\rfloor,\lfloor m/dt\rfloor,t)

对 S 在 t 维度上做前缀和，即可整除分块。

我们预处理出 n,m\leq B 的 S(n,m,k)，由 m*k,n*k<=L 可以得到总状态量为 :

令 B≈\sqrt{L} ,预处理 S 的时空复杂度为 O(L\sqrt{L}) 。

设 f(n,m)=\sum\limits_{t=1}^{n}S(\lfloor n/t\rfloor,\lfloor m/t\rfloor,t)

当 n/t\geq B\Rightarrow t\leq n/B 时，没有预处理 S ，需要暴力求和,复杂度为 O(n/B)\leq O(\sqrt{n})。

当 t>n/B 时，可以整除分块计算，这部分复杂度为 O(\sqrt{n})。

直接分块套分块计算，复杂度为O(∫_1^{\sqrt{n}}\sqrt{n/x})=O(n^{3/4}) ,常数较大。

总复杂度 O(TL^{3/4}+L\sqrt{L})。

#include<algorithm>
#include<cstdio>
#define ll long long
#define MaxN 30500
using namespace std;
const int mod=23333,B=200,lim=30000;
int tn,p[MaxN>>2],mu[MaxN],phi[MaxN]
   ,*G[MaxN],_G[MaxN*16],*gp=_G
   ,*S[B+5][B+5],_S[MaxN*(B+B/4)],*sp=_S;
void getsth()
{
  phi[1]=mu[1]=1;
  for (int i=2,t;i<=lim;i++){
    if (!phi[i]){
      mu[p[++tn]=i]=-1;
      phi[i]=i-1;
    }
    for (int j=1;j<=tn&&(t=p[j]*i)<=lim;j++){
      if (i%p[j]==0){
        phi[t]=phi[i]*p[j];
        break;
      }mu[t]=-mu[i];
      phi[t]=phi[i]*(p[j]-1);
    }
  }G[0]=gp;gp+=lim+1;
  for (int x=1;x<=lim;x++){
    G[x]=gp+1;
    gp+=lim/x+1;
    for (int y=1;x*y<=lim;y++)
      G[x][y]=(G[x-1][y]+phi[x*y])%mod;
  }
  for (int j=1;j<=B;j++)
    for (int i=j;i<=B;i++){
      S[i][j]=sp;
      sp+=lim/i+2;
      for (int k=1;k*i<=lim;k++)
        S[i][j][k]=(S[i][j][k-1]+mu[k]*G[i][k]%mod*G[j][k])%mod;
    }
}
int f(int n,int m)
{
  ll ans=0;
  for (int i=1;i*B<=n;i++)
    ans+=mu[i]*G[n/i][i]*G[m/i][i];
  ans%=mod;
  register int l=n/B+1,r;int tn,tm;
  for (;l<=m;l=r+1){
    r=min(n/(tn=n/l),m/(tm=m/l));
    ans+=S[tn][tm][r]-S[tn][tm][l-1];
  }return ans%mod;
}
int T,N,M;
int main()
{
  scanf("%d",&T);
  getsth();
  while(T--){
    scanf("%d%d",&N,&M);
    if (N<M)swap(N,M);
    ll ans=0;
    int l=1,r;
    for (;l<=M;l=r+1){
      r=min(N/(N/l),M/(M/l));
      ans+=f(N/l,M/l)*(r-l+1)%mod;
    }printf("%d\n",(ans%mod+mod)%mod);
  }return 0;
}

• \Large \oplus\ 20'

发现预处理所有 f(n,n) 的值，就只需一次外层整除分块。

(此外, S 也不用费心)

• \Large 100'

承接上面的式子:

设 T=dk ,交换和式得到:

仍然不能整除分块，暴力的话需要 O(n\log n) ,不知T到哪里去了。

我们设 R(n,m,t)=\sum\limits_{k|T}\mu(k)G(n,k)G(m,k)

我们预处理 n\leq\sqrt{L} 的 R(n,m,t) ，状态量仍然是 O(L\sqrt{L})。

预处理时，对于 t 这一维做约数求和，容易做到 O(L\sqrt{L}\log L) (调和级数)。

实际上还可以做到 O(L\sqrt{L}\log\log L) ，感兴趣的同学可见 狄利克雷相关。

对于 n/T\leq B 的部分整除分块，复杂度为 O(\sqrt{n}).

对于 n/T>B 即 T\leq n/B 的部分我们是没有预处理的,复杂度为 O(n/B\log n)=O(\sqrt{n}\log n)。

综上,空间 O(n\sqrt{n}) ,时间 O\big((L+T)\sqrt{L}\log L\big)。

常数比较小,所以比60pts做法快到不知哪里去了

#include<algorithm>
#include<cstdio>
#define ll long long
#define MaxN 50500
using namespace std;
const int mod=23333,B=150,lim=50000;
int tn,p[MaxN>>2],mu[MaxN],phi[MaxN]
   ,*G[MaxN],_G[MaxN*16],*gp=_G
   ,*R[B+5][B+5],_R[MaxN*(B+B/4)],*rp=_R;
void getsth()
{
  phi[1]=mu[1]=1;
  for (int i=2,t;i<=lim;i++){
    if (!phi[i]){
      mu[p[++tn]=i]=-1;
      phi[i]=i-1;
    }
    for (int j=1;j<=tn&&(t=p[j]*i)<=lim;j++){
      if (i%p[j]==0){
        phi[t]=phi[i]*p[j];
        break;
      }mu[t]=-mu[i];
      phi[t]=phi[i]*(p[j]-1);
    }
  }G[0]=gp;gp+=lim+1;
  for (int x=1;x<=lim;x++){
    G[x]=gp+1;
    gp+=lim/x+1;
    for (int y=1;x*y<=lim;y++)
      G[x][y]=(G[x-1][y]+phi[x*y])%mod;
  }
  for (int j=1;j<=B;j++)
   for (int i=j;i<=B;i++){
     R[i][j]=rp;
     rp+=lim/i+2;
     int n0=lim/i;
     for (int k=1;k<=n0;k++){
       int sav=mu[k]*G[i][k]%mod*G[j][k];
       for (int p=k;p<=n0;p+=k)
         R[i][j][p]+=sav;
     }
     for (int k=1;k*i<=lim;k++)
       R[i][j][k]=(R[i][j][k-1]+R[i][j][k])%mod;
  }
}
int T,N,M;
int main()
{
  scanf("%d",&T);
  getsth();
  while(T--){
    scanf("%d%d",&N,&M);
    long long ans=0;
    for (int i=1;i*B<=N;i++)
      for (int j=i;j*B<=N;j+=i)
        ans+=mu[i]*G[N/j][i]*G[M/j][i];
    ans%=mod;
    register int l=N/B+1,r;int tn,tm;
    for (;l<=M;l=r+1){
      r=min(N/(tn=N/l),M/(tm=M/l));
      ans+=R[tn][tm][r]-R[tn][tm][l-1];
    }printf("%d\n",(ans%mod+mod)%mod);
  }return 0;
}

另解

来自 @流风之回雪 大佬的赛时解法。

发现此题并没有强制在线，可以使用基于莫队的的做法。

静待大佬发题解……

（然而一年多过去了还是没有其他题解，那我就自己写写吧）

设一行的权值 S(n,k)=\sum\limits_{i=1}^n\varphi\Big(\dfrac{{\rm lcm}(i,k)}{\gcd(i,k)}\Big)

莫队可以把问题转化成 O(n\sqrt{T}) 个 S(n,k) 的求解。

到这里可以尝试 O\big((I*d)\big(k)) 枚举约数计算。

平均复杂度大概是 O(\log^2n) 的，这也正是大佬的赛时做法，现在似乎已经无法通过。