题面
https://www.luogu.com.cn/problem/P1829
题解
前置知识
- 线性筛积性函数 https://www.cnblogs.com/zhoushuyu/p/8275530.html
- 莫比乌斯反演：《具体数学：计算机科学基础》4.9节
要求({sum_{i=1}^{n}}{sum_{j=1}^{m}}[i,j])。
首先转化成为({sum_{i=1}^{n}}{sum_{j=1}^{m}}{frac{ij}{(i,j)}})。
求和中化解gcd的小技巧
假设我们遇到求和({sum_{i=1}^{n}}{sum_{j=1}^{m}}f((i,j)))，其中f是某给定函数。
那么我们可以设函数(g = f { imes} {mu})，则由莫比乌斯反演，有(f=g{ imes}I)。（此处及以下“({ imes})”表示狄利克雷卷积）然后做如下转化：
[{sum_{i=1}^{n}}{sum_{j=1}^{m}}f((i,j))
]
[={sum_{i=1}^{n}}{sum_{j=1}^{m}}{sum_{d|(i,j)}}g(d)
]
由于有(d|(i,j))等价于(d|i)且(d|j)的性质，可以进一步化简为：
[={sum_{i=1}^{n}}{sum_{j=1}^{m}}{sum_{d|i,d|j}}g(d)
]
[={sum_d}g(d){sum_{i=1}^{n}}{sum_{j=1}^{m}}[d|i][d|j]
]
[={sum_d}g(d)({sum_{i=1}^{n}}[d|i])({sum_{j=1}^{m}}[d|j])
]
[={sum_d}g(d){{lfloor}{frac{n}{d}}{
floor}}{lfloor}{frac{m}{d}}{
floor}
]
即完成了gcd的化简。
如果求和是({sum_{i=1}^{n}}{sum_{j=1}^{m}}g(i)h(j)f((i,j)))，过程也如出一辙，只要
[g(d)+g(2d)+…+g({lfloor}{frac{n}{d}}{
floor}d)
]
[h(d)+h(2d)+…+h({lfloor}{frac{m}{d}}{
floor}d)
]
能够快速处理，也能够适用此方法。
回本题。此题中(f(x) = {frac{1}{x}})。
[{sum_{i=1}^{n}}{sum_{j=1}^{m}}ijf((i,j))
]
[={sum_{i=1}^{n}}{sum_{j=1}^{m}}ij{sum_{d|(i,j)}}g(d)
]
[={sum_{i=1}^{n}}{sum_{j=1}^{m}}ij{sum_{d|i,d|j}}g(d)
]
[={sum_{d}}g(d){sum_{i=1}^{n}}{sum_{j=1}^{m}}[d|i][d|j]ij
]
[={sum_{d}}g(d)({sum_{i=1}^{n}}[d|i]i)({sum_{j=1}^{n}}[d|j]j)
]
注意到两个括号内分别是
[d+2d+…+{lfloor}{frac{n}{d}}{
floor}d
]
[d+2d+…+{lfloor}{frac{m}{d}}{
floor}d
]
这两个等差数列之和，于是原式=
[{sum_d}g(d)d^2({frac{1}{2}}{lfloor}{frac{n}{d}}{
floor}({lfloor}{frac{n}{d}}{
floor}+1))({frac{1}{2}}{lfloor}{frac{m}{d}}{
floor}({lfloor}{frac{m}{d}}{
floor}+1))
]
其中
[g(d)={sum_{p|d}}f({frac{d}{p}}){mu}(p)
]
[={sum_{p|d}}{frac{p}{d}}{mu(p)}
]
为了避免出现小数或者分数，简化计算，我们令(h(d)=d*g(d))
那么有
[h(d)={sum_{p|d}}p{mu(p)}
]
这是一个积性函数，可以线性求出。同时有原式=
[{sum_d}h(d)d({frac{1}{2}}{lfloor}{frac{n}{d}}{
floor}({lfloor}{frac{n}{d}}{
floor}+1))({frac{1}{2}}{lfloor}{frac{m}{d}}{
floor}({lfloor}{frac{m}{d}}{
floor}+1))
]
可以枚举d，O(min(n,m))求解。
代码
#include<bits/stdc++.h>
using namespace std;
#define N 10000000
#define rg register
#define mod 20101009
#define F(x) (1ll*x*(x+1)/2%mod)
namespace ModCalc{
inline void Inc(int &x,int y){
x += y;if(x > mod)x -= mod;
}
inline void Dec(int &x,int y){
x -= y;if(x < 0)x += mod;
}
inline void Tms(int &x,int y){
x = 1ll * x * y % mod;
}
inline int Add(int x,int y){
Inc(x,y);return x;
}
inline int Sub(int x,int y){
Dec(x,y);return x;
}
inline int Mul(int x,int y){
Tms(x,y);return x;
}
}
using namespace ModCalc;
int pn;
int pri[1200000+5];
bool isp[N+5];
int h[N+5];
int P[N+5]; //P[i]表示i的素因数分解式中底数最小的那一项的大小
inline void Eular(){ //线性筛h[]
pn = 0;
h[1] = 1;
for(rg int i = 2;i <= N;i++)isp[i] = 1;
for(rg int i = 2;i <= N;i++){
if(isp[i]){
pri[++pn] = i;
P[i] = i;
h[i] = mod + 1 - i;
}
for(rg int j = 1;i * pri[j] <= N;j++){
isp[i*pri[j]] = 0;
if(i % pri[j]){
P[i*pri[j]] = pri[j];
h[i*pri[j]] = Mul(h[i],h[pri[j]]);
}
else{
int I = i * pri[j];
P[I] = Mul(P[i],pri[j]);
if(P[I] == I)h[I] = h[i];
else h[I] = Mul(h[I/P[I]],h[P[I]]);
break;
}
}
}
}
int n,m;
int main(){
Eular();
scanf("%d%d",&n,&m);
int ans = 0;
int M = min(n,m);
for(rg int d = 1;d <= M;d++)Inc(ans,1ll*h[d]*d%mod*F(n/d)%mod*F(m/d)%mod);
cout << ans << endl;
return 0;
}