设 d(x) 为 x 的约数个数,给定 N,M ,求
∑i=1N∑j=1Md(ij)。
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
这题有个很屌的结论:
∑i=1N∑j=1Md(ij)=∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋[gcd(i,j)=1]
根据 PoPoQQQ博客所说的,我们可以先证明这个式子的成立:
d(nm)=∑i∣n∑j∣m[gcd(i,j)=1]
我们可以证明一下:我们对每个质数 p 单独算贡献,设 n=n′×pk1 , m=m′×pk2 。那么,该质数 p 对答案的贡献显然为 k1+k2+1 。于是我们考虑
d(nm)=∑i∣n∑j∣m[gcd(i,j)=1]
这个式子,发现
p
对之有贡献的数对
(i,j)
仍然是
(pk1,1),(pk1−1,1)⋯(1,1)⋯(1,pk2−1),(1,pk2)
这
k1+k2+1
个,因此得证。
代入得
∑n=1N∑m=1Md(nm)=∑n=1N∑m=1M∑i∣n∑j∣m[gcd(i,j)=1]
我们转变枚举量,先枚举 i,j 就有
∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋[gcd(i,j)=1]
于是
∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋[gcd(i,j)=1]
这个式子我们可以上反演了。
反演化为
∑i=1N∑j=1M⌊Ni⌋⌊Mj⌋∑g∣ig∣jμ(g)
转而枚举 g ,于是就可得到
∑g=1Nμ(g)∑i=1⌊Ng⌋∑j=1⌊Mg⌋⌊Nig⌋⌊Mjg⌋
再化一下就可得到
∑g=1Nμ(g)∑i=1⌊Ng⌋⌊Nig⌋∑j=1⌊Mg⌋⌊Mjg⌋
又有
⌊Nab⌋=⌊⌊Na⌋b⌋
于是我们发现 ∑⌊Ng⌋i=1⌊Nig⌋ 只与 ⌊Ng⌋ 有关,我们可以 O(nn−√) 预处理
fx=∑i=1x⌊xi⌋
有了这个后再化简
∑g=1Nμ(g)f⌊Ng⌋f⌊Mg⌋
就可在
O(n−√)
分段求了。皆大欢喜。
#include<iostream>
#include<cstdio>
#include<cstdlib>
using namespace std;
typedef long long ll;
#define maxn (50010)
int f[maxn],mu[maxn],prime[maxn],n,m,tot; bool exist[maxn];
inline int calc(int x)
{
int ret = 0;
for (int i = 1,last;i <= x;i = last+1)
{
last = min(x,x/(x/i));
ret += (x/i)*(last-i+1);
}
return ret;
}
inline void ready()
{
mu[1] = 1;
for (int i = 2;i <= 50000;++i)
{
if (!exist[i]) { prime[++tot] = i; mu[i] = -1; }
for (int j = 1;j <= tot&&prime[j]*i <= 50000;++j)
{
exist[i*prime[j]] = true;
if (i % prime[j] == 0) { mu[i*prime[j]] = 0; break; }
mu[i*prime[j]] = -mu[i];
}
}
for (int i = 1;i <= 50000;++i) mu[i] += mu[i-1],f[i] = calc(i);
}
inline ll work()
{
if (n > m) swap(n,m);
ll ret = 0;
for (int i = 1,last;i <= n;i = last+1)
{
last = min(n,min(n/(n/i),m/(m/i)));
ret += (ll)(mu[last]-mu[i-1])*((ll)f[n/i]*f[m/i]);
}
return ret;
}
int main()
{
freopen("3994.in","r",stdin);
freopen("3994.out","w",stdout);
ready();
int T; scanf("%d",&T);
while (T--) scanf("%d %d",&n,&m),printf("%lld\n",work());
fclose(stdin); fclose(stdout);
return 0;
}