Description
设d(x)为x的约数个数,给定N、M,求(sumlimits_{i=1}^nsumlimits_{j=1}^md(ij))
Input
输入文件包含多组测试数据。
第一行,一个整数T,表示测试数据的组数。
接下来的T行,每行两个整数N、M。
Output
T行,每行一个整数,表示你所求的答案。
Sample Input
2
7 4
5 6
Sample Output
110
121
HINT
(1leqslant N, Mleqslant50000)
(1leqslant Tleqslant50000)
(d)太丑了,我之后都写成(sigma_0)了
首先你需要知道(sigma_0(nm)=sumlimits_{i|n}sumlimits_{j|m}[gcd(i,j)=1])
(我就是不知道这个推了1h无果)
然后我们证明一下这个东西(我只能证明等式左右两边相等)
令(n=prodlimits_{i=1}^kp_i^{a_i},m=prodlimits_{i=1}^kp_i^{b_i}),(a_i,b_i)可以为0
然后有(sigma_0(nm)=prodlimits_{i=1}^k(a_i+b_i+1)),我们考虑某个质因子(p_i)对答案的贡献, (p_i)要么在(n)中有(a_i)中取法,要么在(m)中有(b_i)种取法,要么不取有一种方法,所以共计(a_i+b_i+1)中取法,再乘法原理累乘一下即为(sigma_0(n,m))
然后我们来推式子
我们调整一下枚举顺序
(u,v)太丑啦,我们改一下变量名,并且对最后用面的式子反演一下
我们把枚举(x)提到前面,并且进行分组,得到
如果我们能够求出(sumlimits_{i=1}^nlfloordfrac{n}{i} floor),我们就可以对原式进行数论分块了
然后我们参考一下[AHOI2005]约数研究这题,可以发现(sumlimits_{i=1}^nlfloordfrac{n}{i} floor=sumlimits_{i=1}^nsigma_0(i)),因为(lfloordfrac{n}{i} floor)相当于统计(1sim n)中有多少个数的因数有(i)
(sigma_0(n))显然是个积性函数,因此我们可以线筛求出前缀和,然后对原式进行数论分块即可
/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline char gc(){
static char buf[1000000],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
inline int frd(){
int x=0,f=1; char ch=gc();
for (;ch<'0'||ch>'9';ch=gc()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=gc()) x=(x<<3)+(x<<1)+ch-'0';
return x*f;
}
inline int read(){
int x=0,f=1; char ch=getchar();
for (;ch<'0'||ch>'9';ch=getchar()) if (ch=='-') f=-1;
for (;ch>='0'&&ch<='9';ch=getchar()) x=(x<<3)+(x<<1)+ch-'0';
return x*f;
}
inline void print(int x){
if (x<0) putchar('-'),x=-x;
if (x>9) print(x/10);
putchar(x%10+'0');
}
const int N=5e4;
int prime[N+10],mu[N+10],smu[N+10],f[N+10],Frm[N+10];
ll sf[N+10];
bool inprime[N+10];
void prepare(){
int tot=0;
mu[1]=f[1]=smu[1]=sf[1]=1;
for (int i=2;i<=N;i++){
if (!inprime[i]) prime[++tot]=i,mu[i]=-1,f[i]=Frm[i]=2;
for (int j=1;j<=tot&&i*prime[j]<=N;j++){
inprime[i*prime[j]]=1;
if (i%prime[j]==0){
mu[i*prime[j]]=0;
Frm[i*prime[j]]=Frm[i]+1;
f[i*prime[j]]=f[i]/Frm[i]*Frm[i*prime[j]];
break;
}
mu[i*prime[j]]=-mu[i];
Frm[i*prime[j]]=2;
f[i*prime[j]]=f[i]<<1;
}
sf[i]=sf[i-1]+f[i];
smu[i]=smu[i-1]+mu[i];
}
}
int main(){
prepare();
for (int T=read();T;T--){
int n=read(),m=read(); ll Ans=0;
if (n>m) swap(n,m);
for (int i=1,pos;i<=n;i=pos+1){
pos=min(n/(n/i),m/(m/i));
Ans+=(smu[pos]-smu[i-1])*sf[n/i]*sf[m/i];
}
printf("%lld
",Ans);
}
return 0;
}