【SDOI2015】BZOJ3994 约数个数和题解(Mobius反演+除法分块)

题目:BZOJ3994.
题目大意:设 d ( n ) = ∑ d ∣ n 1 d(n)=\sum_{d|n}1 d(n)=dn1,给定 T T T组询问 n , m n,m n,m,每次询问 ∑ i = 1 n ∑ j = 1 m d ( i j ) \sum_{i=1}^{n}\sum_{j=1}^{m}d(ij) i=1nj=1md(ij).
1 ≤ n , m , T ≤ 5 ∗ 1 0 4 1\leq n,m,T\leq 5*10^4 1n,m,T5104.

我们先考虑,如何在知道 i i i j j j的情况下如何用 O ( i j ) O(ij) O(ij)的时间复杂度求 d ( i j ) d(ij) d(ij).

考虑知道 i j ij ij后,由于 d d d其实就是约数个数函数,而我们知道,对于任意一个 n = ∏ i = 1 k p i c i n=\prod_{i=1}^{k}p_i^{c_i} n=i=1kpici,那么 d ( n ) = ∏ i = 1 k ( c i + 1 ) d(n)=\prod_{i=1}^{k}(c_i+1) d(n)=i=1k(ci+1).

对于两个数 n n n m m m,考虑如何计算 d ( n m ) d(nm) d(nm).首先我们把两个数分别写成 n = ∏ i = 1 k p i c i , m = ∏ i = 1 k p i q i n=\prod_{i=1}^{k}p_i^{c_i},m=\prod_{i=1}^{k}p_i^{q_i} n=i=1kpici,m=i=1kpiqi,那么:
d ( n m ) = d ( ∏ i = 1 k p i c i + q i ) = ∏ i = 1 k ( c i + q i + 1 ) d(nm)=d(\prod_{i=1}^{k}p_i^{c_i+q_i})=\prod_{i=1}^{k}(c_i+q_i+1) d(nm)=d(i=1kpici+qi)=i=1k(ci+qi+1)

对于每一个 ( c i + q i + 1 ) (c_i+q_i+1) (ci+qi+1),我们构造一个 ( 0 , 0 ) (0,0) (0,0) ( c i , q i ) (c_i,q_i) (ci,qi)的矩阵,那么我们会发现 ( c i + q i + 1 ) (c_i+q_i+1) (ci+qi+1)等于矩阵中的一个"L"形的区域,即所有坐标中含有 0 0 0的总数.我们发现罪域一个坐标 ( x , y ) (x,y) (x,y)满足 x = 0 x=0 x=0 y = 0 y=0 y=0时, g c d ( p i x , p i y ) = 1 gcd(p_i^x,p_i^y)=1 gcd(pix,piy)=1,所以可以得到:
d ( n m ) = ∏ i = 1 k ∑ j = 0 c i ∑ t = 0 q i [ g c d ( p i j , p i t ) = 1 ] = ∑ i ∣ n ∑ j ∣ m [ g c d ( i , j ) = 1 ] d(nm)=\prod_{i=1}^{k}\sum_{j=0}^{c_i}\sum_{t=0}^{q_i}[gcd(p_i^j,p_i^t)=1]=\sum_{i|n}\sum_{j|m}[gcd(i,j)=1] d(nm)=i=1kj=0cit=0qi[gcd(pij,pit)=1]=injm[gcd(i,j)=1]

得到了这个式子后,我们就可以开始推导原式:
∑ i = 1 n ∑ j = 1 m d ( i j ) = ∑ i = 1 n ∑ j = 1 m ∑ x ∣ i ∑ y ∣ j [ g c d ( x , y ) = 1 ] = ∑ x = 1 n ∑ y = 1 m ⌊ n x ⌋ ⌊ m y ⌋ ∑ t ∣ x ∧ t ∣ y μ ( t ) = ∑ t = 1 m i n ( n , m ) μ ( t ) ∑ x = 1 ⌊ n t ⌋ ⌊ n t x ⌋ ∑ y = 1 ⌊ m t ⌋ ⌊ m t y ⌋ \sum_{i=1}^{n}\sum_{j=1}^{m}d(ij)\\ =\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x|i}\sum_{y|j}[gcd(x,y)=1]\\ =\sum_{x=1}^{n}\sum_{y=1}^{m}\left\lfloor\frac{n}{x}\right\rfloor\left\lfloor\frac{m}{y}\right\rfloor\sum_{t|x\wedge t|y}\mu(t)\\ =\sum_{t=1}^{min(n,m)}\mu(t)\sum_{x=1}^{\left\lfloor\frac{n}{t}\right\rfloor}\left\lfloor\frac{n}{tx}\right\rfloor\sum_{y=1}^{\left\lfloor\frac{m}{t}\right\rfloor}\left\lfloor\frac{m}{ty}\right\rfloor i=1nj=1md(ij)=i=1nj=1mxiyj[gcd(x,y)=1]=x=1ny=1mxnymtxtyμ(t)=t=1min(n,m)μ(t)x=1tntxny=1tmtym

我们设 g ( n ) = ∑ i = 1 n ⌊ n i ⌋ g(n)=\sum_{i=1}^{n}\left\lfloor\frac{n}{i}\right\rfloor g(n)=i=1nin,然后可以得到:
= ∑ t = 1 m i n ( n , m ) μ ( t ) ∑ x = 1 ⌊ n t ⌋ ⌊ n t x ⌋ ∑ y = 1 ⌊ m t ⌋ ⌊ m t y ⌋ = ∑ t = 1 m i n ( n , m ) μ ( t ) g ( ⌊ n t ⌋ ) g ( ⌊ m t ⌋ ) =\sum_{t=1}^{min(n,m)}\mu(t)\sum_{x=1}^{\left\lfloor\frac{n}{t}\right\rfloor}\left\lfloor\frac{n}{tx}\right\rfloor\sum_{y=1}^{\left\lfloor\frac{m}{t}\right\rfloor}\left\lfloor\frac{m}{ty}\right\rfloor\\ =\sum_{t=1}^{min(n,m)}\mu(t)g(\left\lfloor\frac{n}{t}\right\rfloor)g(\left\lfloor\frac{m}{t}\right\rfloor) =t=1min(n,m)μ(t)x=1tntxny=1tmtym=t=1min(n,m)μ(t)g(tn)g(tm)

然后我们看看 g g g函数的本质是什么.复读机!
g ( n ) = ∑ i = 1 n ⌊ n i ⌋ = ∑ i = 1 n ∑ i ∣ j ∧ j ∈ [ 1 , n ] 1 = ∑ j = 1 n ∑ i ∣ j 1 = ∑ j = 1 n d ( j ) g(n)=\sum_{i=1}^{n}\left\lfloor\frac{n}{i}\right\rfloor\\ =\sum_{i=1}^{n}\sum_{i|j\wedge j\in[1,n]}1\\ =\sum_{j=1}^{n}\sum_{i|j}1\\ =\sum_{j=1}^{n}d(j) g(n)=i=1nin=i=1nijj[1,n]1=j=1nij1=j=1nd(j)

发现 g g g函数其实就是 d d d函数的前缀和,而我们知道约数个数函数是可以直接线性筛的.所以预处理出 g g g函数,然后除法分块一下即可,时间复杂度 O ( n + T n ) O(n+T\sqrt{n}) O(n+Tn ).

代码如下:

#include<bits/stdc++.h>
  using namespace std;

#define Abigail inline void
typedef long long LL;

const int N=50000;

int pr[N+9],tp,b[N+9];
LL d[N+9],num[N+9],mu[N+9];

void sieve(int n){
  for (int i=2;i<=n;++i) b[i]=1;
  d[1]=mu[1]=num[1]=1;int v;
  for (int i=2;i<=n;++i){
  	if (b[i]) pr[++tp]=i,d[i]=2,num[i]=1,mu[i]=-1;
  	for (int j=1;j<=tp&&i*pr[j]<=n;++j){
  	  b[v=i*pr[j]]=0;
  	  if (i%pr[j]==0){
  	    num[v]=num[i]+1;d[v]=d[i]/num[v]*(num[v]+1);
  	    mu[v]=0;
  	    break;
  	  }
  	  num[v]=1;d[v]=d[i]<<1;
  	  mu[v]=-mu[i];
  	}
  }
}

LL get(int n,int m){
  LL ans=0;
  if (n>m) swap(n,m);
  for (int l=1,r;l<=n;l=r+1){
  	r=min(n/(n/l),m/(m/l));
  	ans+=d[n/l]*d[m/l]*(mu[r]-mu[l-1]);
  }
  return ans;
} 

int n,m;

Abigail start(){
  sieve(N);
  for (int i=1;i<=N;++i)
    d[i]+=d[i-1],mu[i]+=mu[i-1];
}

Abigail into(){
  scanf("%d%d",&n,&m);
}

Abigail outo(){
  printf("%lld\n",get(n,m));
}

int main(){
  start();
  int T;
  scanf("%d",&T);
  while (T--){
  	into();
  	outo();
  }
  return 0;
}
  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值