题目:http://www.lydsy.com/JudgeOnline/problem.php?id=4018
题意:
T
组询问,每组询问给定
T≤1000,N,M≤2∗106
。
题解:
考虑将
∑Ni=1∑Mj=1|i−j|gcd(i,j)
化简。
首先引入一个记号,令
[x]
表示:若
x
为真,则
∑i=1N∑j=1M|i−j|gcd(i,j)=∑i=1N∑j=1M∑d|i−j|gcd(i,j)[d=gcd(i,j)]=∑d∑i=1⌊Nd⌋∑j=1⌊Md⌋|i−j|[gcd(i,j)=1]=∑d∑i=1⌊Nd⌋∑j=1⌊Md⌋|i−j|∑d′|gcd(i,j)μ(d′)=∑d∑d′d′μ(d′)∑i=1⌊Ndd′⌋∑j=1⌊Mdd′⌋|i−j|=∑d∑d′|dd′μ(d′)∑i=1⌊Nd⌋∑j=1⌊Md⌋|i−j|
(最后一步是将 dd′ 用 d 代换了,所以才会有
后面那一坨只和 Nd,Md 有关系的式子的值是:设 A=min(Nd,Md),B=max(Nd,Md) ,则 ∑Ai=1∑Bj=1|i−j|=(A−1)A(A+1)3+AB(B−A)2 ,这里懒得推导,存在更优美的式子。
而前面那一坨,令 f(n)=∑d|ndμ(d) ,由积性函数的性质可知,它也是积性函数,我们可以考虑筛出它的值。 f(1)=1 。
若 n 为质数,则
若 n 的最小质因数
若 n 的最小质因数
对于询问,容易看出 (⌊Nd⌋,⌊Md⌋) 的取值最多只有 2⌊N−−√⌋+2⌊M−−√⌋ 种,且取值对应的 d 是一段连续的区间,可以将
代码:
#include <cstdio>
#include <algorithm>
using namespace std;
const int maxn = 2000001, mod1 = 1000000007, mod2 = 1000000009;
int t, n, m, p, q, tot, prime[maxn], sf[maxn][2], ans[2];
bool vis[maxn];
void inc(int &x, int y, int mod)
{
x += y;
if(x >= mod)
x -= mod;
}
int f(int A, int B, int mod)
{
if(!A || !B)
return 0;
if(A > B)
swap(A, B);
int cnt1 = (long long)(A - 1) * A * (A + 1) / 3 % mod;
int cnt2 = (long long)A * B * (B - A) / 2 % mod;
inc(cnt1, cnt2, mod);
return cnt1;
}
int main()
{
sf[1][0] = sf[1][1] = 1;
for(int i = 2; i < maxn; ++i)
{
if(!vis[i])
{
prime[tot++] = i;
sf[i][0] = mod1 - i + 1;
sf[i][1] = mod2 - i + 1;
}
for(int j = 0; j < tot && (long long)i * prime[j] < maxn; ++j)
{
vis[i * prime[j]] = 1;
if(i % prime[j] == 0)
{
sf[i * prime[j]][0] = sf[i][0];
sf[i * prime[j]][1] = sf[i][1];
break;
}
else
{
sf[i * prime[j]][0] = (long long)sf[i][0] * sf[prime[j]][0] % mod1;
sf[i * prime[j]][1] = (long long)sf[i][1] * sf[prime[j]][1] % mod2;
}
}
}
for(int i = 2; i < maxn; ++i)
{
inc(sf[i][0], sf[i - 1][0], mod1);
inc(sf[i][1], sf[i - 1][1], mod2);
}
scanf("%d", &t);
while(t--)
{
scanf("%d%d", &n, &m);
if(n > m)
swap(n, m);
ans[0] = ans[1] = 0;
for(int i = 1, j; i <= n; i = j + 1)
{
j = min(n / (n / i), m / (m / i));
inc(ans[0], (long long)(sf[j][0] - sf[i - 1][0] + mod1) * f(n / i, m / i, mod1) % mod1, mod1);
inc(ans[1], (long long)(sf[j][1] - sf[i - 1][1] + mod2) * f(n / i, m / i, mod2) % mod2, mod2);
}
printf("%d %d\n", ans[0], ans[1]);
}
return 0;
}