[SDOI2015]約數個數和
題意
求\(\sum\limits_i^n\ \sum\limits_j^m\ d(ij)\)其中\(d(x)\)為\(x\)的約數個數
想法
看到這個柿子想到的就是莫反了。
有這樣一個結論
\( d(xy)= \sum\limits_{i|x} \ \sum\limits_{j|y}[gcd(i,j) == 1] \)
所以變形原柿子(n > m)
\(\sum\limits_i^n\ \sum\limits_j^m\ d(ij)\)
\(\sum\limits_i^n\ \sum\limits_j^m\ \sum\limits_{x|i} \ \sum\limits_{y|j}[gcd(x,y) == 1])\)
用\(\mu(x)\) 的性質:
\([n == 1] = \sum\limits_{d|n}\mu(d)\)
\(\sum\limits_i^n\ \sum\limits_j^m\ \sum\limits_{x|i} \ \sum\limits_{y|j}\sum\limits_{d|gcd(x,y)}\mu(d))\)
改為列舉\(d\)-
\(\sum\limits_{d = 1}^n\sum\limits_{x = 1}^{n/d}\sum\limits_{i = 1}^{n/d/x}\sum\limits_{y = 1}^{m/d}\sum\limits_{j = 1}^{m/d/y}\mu(d)\)
\(\sum\limits_{d = 1}^n\sum\limits_{x = 1}^{n/d}\lfloor\frac{n}{d * x}\rfloor\sum\limits_{y = 1}^{n/d}\lfloor\frac{m}{d * y}\rfloor\mu(d)\)
\(\sum\limits_{d = 1}^n\mu(d)\sum\limits_{x = 1}^{n/d}\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{x}\rfloor\sum\limits_{y = 1}^{m/d}\lfloor\frac{\lfloor\frac{m}{d}\rfloor}{y}\rfloor\)
後面兩個\(\sum\)可以整除分塊
複雜度\(O(N\sqrt(N))\)
程式碼
#include<bits/stdc++.h>
using namespace std;
#define int long long
int read(){
char cc = getchar(); int cn = 0, flus = 1;
while(cc < '0' || cc > '9'){
if(cc == '-') flus = -flus;
cc = getchar();
}
while(cc >= '0' && cc <= '9')
cn = cn * 10 + cc - '0', cc = getchar();
return cn * flus;
}
const int N = 50000 + 5;
int u[N], sum[N], pm[N], n, f[N], top;
bool isp[N];
void init(){
u[1] = 1;
for(int i = 2; i <= n; i++){
if(!isp[i]){
pm[++top] = i;
u[i] = -1;
}
for(int j = 1; j <= top; j++){
int p = pm[j];
if(p * i > n) break;
if(i % p == 0){
isp[p * i] = 1;
break;
}
isp[p * i] = 1;
u[p * i] = u[p] * u[i];
}
}
int l, r;
for(int i = 1; i <= n; i++){
for(l = 1; l <= i; l = r + 1){
r = i/(i/l);
f[i] += (r - l + 1) * (i/l);
}
sum[i] = sum[i - 1] + u[i];
}
}
int solve(int x, int y){
if(x > y) swap(x, y);
int ans = 0, l = 1, r;
for(l = 1; l <= x; l = r + 1){
r = min((x/(x/l)), (y/(y/l)));
ans += (sum[r] - sum[l - 1]) * f[x/l] * f[y/l];
}
return ans;
}
signed main()
{
n = 50000;
init();
int T = read(), a, b;
while(T--){
a = read(); b = read();
printf("%lld\n", solve(a, b));
}
return 0;
}