洛谷P3327 [SDOI2015]約數個數和 【莫比烏斯反演】
題目
設d(x)為x的約數個數,給定N、M,求\(\sum_{i = 1}^{N} \sum_{j = 1}^{M} d(ij)\)
輸入格式
輸入文件包含多組測試數據。第一行,一個整數T,表示測試數據的組數。接下來的T行,每行兩個整數N、M。
輸出格式
T行,每行一個整數,表示你所求的答案。
輸入樣例
2
7 4
5 6
輸出樣例
110
121
提示
1<=N, M<=50000
1<=T<=50000
題解
好神的題【是我太弱吧】
首先上來就傷結論。。
題目所求
\(ans = \sum_{i = 1}^{N} \sum_{j = 1}^{M} d(ij)\)
有一個這樣的結論:
\(d(ij) = \sum_{x|i}\sum_{y|j} [gcd(x,y) == 1]\)
那麽就轉化為了:
\(ans =\sum_{i = 1}^{N} \sum_{j = 1}^{M} \sum_{x|i}\sum_{y|j} [gcd(x,y) == 1]\)
我們考慮對於每一對互質的x、y,x會被枚舉\(\lfloor \frac{N}{x} \rfloor\)次,y會被枚舉\(\lfloor \frac{M}{y} \rfloor\)次
所以有
\(ans =\sum_{i = 1}^{N} \sum_{j = 1}^{M} \lfloor \frac{N}{i} \rfloor \lfloor \frac{M}{j} \rfloor [gcd(i,j) == 1]\)
那麽可以進行莫比烏斯反演了
令
\(f(n) = \sum_{i = 1}^{N} \sum_{j = 1}^{M} \lfloor \frac{N}{i} \rfloor \lfloor \frac{M}{j} \rfloor [gcd(i,j) == n]\)
\(F(n) = \sum_{i = 1}^{N} \sum_{j = 1}^{M} \lfloor \frac{N}{i} \rfloor \lfloor \frac{M}{j} \rfloor [n | gcd(i,j)]\)
那麽有
\(F(n) = (\sum_{i = 1}^{N}\lfloor \frac{N}{i} \rfloor) * (\sum_{j = 1}^{M} \lfloor \frac{M}{j} \rfloor) [n | gcd(i,j)]\)
其中\(\sum_{i = 1}^{N}\lfloor \frac{N}{i} \rfloor\)可以\(O(n\sqrt{n})\)預處理出,我們記為\(sum(n)\)
\(ans = f(1) = \sum_{d = 1}^{N} \mu(d) * F(d) = \sum_{d} \mu(d) sum(\lfloor \frac{N}{d} \rfloor) * sum(\lfloor \frac{M}{d} \rfloor)\)
分塊計算
復雜度\(O(T\sqrt{N} + N\sqrt(N))\)
#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define LL long long int
#define REP(i,n) for (int i = 1; i <= (n); i++)
#define Redge(u) for (int k = h[u],to; k; k = ed[k].nxt)
#define BUG(s,n) for (int i = 1; i <= (n); i++) cout<<s[i]<<‘ ‘; puts("");
using namespace std;
const int maxn = 50005,maxm = 100005,INF = 1000000000;
inline int read(){
int out = 0,flag = 1; char c = getchar();
while (c < 48 || c > 57) {if (c == ‘-‘) flag = -1; c = getchar();}
while (c >= 48 && c <= 57) {out = (out << 3) + (out << 1) + c - ‘0‘; c = getchar();}
return out * flag;
}
int prime[maxn],primei,mu[maxn],isn[maxn];
LL sum[maxn];
void init(){
mu[1] = 1;
for (int i = 2; i < maxn; i++){
if (!isn[i]) prime[++primei] = i,mu[i] = -1;
for (int j = 1; j <= primei && i * prime[j] < maxn; j++){
isn[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 < maxn; i++) mu[i] += mu[i - 1];
for (int n = 1,nxt; n <= 50000; n++){
for (int i = 1; i <= n; i = nxt + 1){
nxt = n / (n / i);
sum[n] += (LL)(nxt - i + 1) * (n / i);
}
}
}
int main(){
init();
int T = read(),n,m;
while (T--){
n = read(); m = read();
if (n > m) swap(n,m);
LL ans = 0; int nxt;
for (int i = 1; i <= n; i = nxt + 1){
nxt = min(n / (n / i),m / (m / i));
ans += sum[n / i] * sum[m / i] * (mu[nxt] - mu[i - 1]);
}
printf("%lld\n",ans);
}
return 0;
}
洛谷P3327 [SDOI2015]約數個數和 【莫比烏斯反演】