BZOJ3994: [SDOI2015]约数个数和

3994: [SDOI2015]约数个数和

Time Limit: 20 Sec  Memory Limit: 128 MB
Submit: 1072  Solved: 736
[Submit][Status][Discuss]

Description

 设d(x)为x的约数个数，给定N、M，求  

 

 

Input

输入文件包含多组测试数据。

第一行，一个整数T，表示测试数据的组数。

接下来的T行，每行两个整数N、M。

 

Output

 T行，每行一个整数，表示你所求的答案。

 

Sample Input

2
7 4
5 6

Sample Output

110
121

HINT

 

 1<=N, M<=50000

1<=T<=50000

 

思路{

　　首先有D(N)=π(ki+1)------①其中ki为质因子pi的出现次数.一个质因子对答案的贡献就是ki+1

　　∑(i|n)∑(j|m)[gcd(i,j)=1]------②对于pi这个因子只能在i或j中出现,故合理的方案有(0,k1), (0,k1-1)...(0,0),(2,0),(3,0),,,(k2,0)有k1+k2+1对和①算出的贡献相等.即①等价于②

　　所以Ans= (n)∑(i=1) (m)∑(j=1) ∑(k|i) ∑(l|j) [gcd(k,l)=1]

　　那么直接把强上莫比乌斯反演:

　　Ans=(n)∑(i=1) (m)∑(j=1) (min(n,m))∑(d=1) (μ(d)* ∑(k|i) ∑(l|j) [gcd(k,l)|d])

　　考虑把莫比乌斯函数和后面的一大坨提到前面.

　　Ans=(min(n,m)) ∑(d=1) (μ(d) (n)∑(k=1) (m)∑(l=1) ( [gcd(k,l)|d]*( (n/d)*(m/d) ) )

　　再去掉gcd,

　　Ans=(min(n,m)) ∑(d=1) (μ(d) (n/d)∑(k=1) (m/d)∑(l=1) ( (n/ (d*k) ) *( (m / (d*l) ) ) )

　　里面可以变成

　　Ans=(min(n,m)) ∑(d=1) μ(d) [ (n/d)∑(k=1) (n/ (d*k) ) ] [ (m/d)∑(l=1) *(m / (d*l) ) ]

　　设F(x)=(x)∑(i=1) (x/i)

　　那么Ans=(min(n,m)) ∑(d=1) μ(d) F(n/d) F(m/d)

　　预处理μ的前缀和和F函数,再数论分块就可以了.

}

#include<bits/stdc++.h>
#define LL long long
#define RG register
#define il inline
#define N 50010
#define LL long long
using namespace std;
LL mu[N],f[N],p[N],n,m;bool vis[N];
void pre(){
  mu[1]=1;
  for(int i=2;i<N;++i){
    if(!vis[i])p[++p[0]]=i,mu[i]=-1;
    for(int j=1;j<=p[0]&&p[j]*i<N;++j){
      vis[i*p[j]]=true;
      if(i%p[j])mu[i*p[j]]=-mu[i];
      else {
    mu[i*p[j]]=0;
    break;
      }
    }
  }
  for(int i=2;i<N;++i)mu[i]+=mu[i-1];
  for(int i=1;i<N;++i){
    for(int l=1,r;l<=i;l=r+1){
      r=i/(i/l);
      f[i]+=(i/l)*(r-l+1);
    }
  }
}
int main(){
  int T;scanf("%d",&T);
  pre();
  while(T--){
    scanf("%lld%lld",&n,&m);
    if(n>m)swap(n,m);
    LL ans(0);
    for(int l=1,r;l<=n;l=r+1){
      r=min(n/(n/l),m/(m/l));
      ans+=(mu[r]-mu[l-1])*(f[n/l]*f[m/l]);
    }printf("%lld\n",ans);
  }return 0;
}