重拾莫比乌斯反演
起因：单纯的不想写代码，做不动题，状态非常差……于是就写点清新莫反。
记 \(f(x)\) 表示 \(x\) 的约数个数。
\[\sum_{i=1}^{n}\sum_{j=1}^{m}f(i)f(j)f(\gcd(i,j))\\
\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=d]f(d)f(i)f(j)\\
\sum_{d=1}^{n}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}[\gcd(i,j)=1]f(d)f(id)f(jd)\\
=\sum_{d=1}^{n}f(d)\sum_{i=1}^{\frac{n}{d}}f(id)\sum_{j=1}^{\frac{m}{d}}f(jd)\sum_{x|i,x|j}\mu(x)\\
=\sum_{d=1}^{n}f(d)\sum_{x=1}^{\frac{n}{d}}\mu(x)\sum_{i=1}^{\frac{n}{dx}}f(idx)\sum_{j=1}^{\frac{m}{dx}}f(jdx)\\
=\sum_{d=1}^{n}f(d)\sum_{x=1}^{\frac{n}{d}}\mu(x)g(dx)h(dx)\\
\]
其中 \(g(x)=\sum\limits_{i=1}^{n}[x|i]f(i),h(x)=\sum\limits_{i=1}^{m}[x|i]f(i)\)。
于是就可以 \(O(n\log n)\) 暴力求了，预处理 \(f,g,h\)，最后统计答案都是调和级数。不过跑的非常的慢，大力卡常之后才贴着时限跑。
Code
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N = 2000005;
int f[N], g[N], h[N], n, m, mod, ans;
int pri[N], pct, mu[N];
bool vis[N];
inline void Sieve(const int &lim = N - 1) {
for(int i = 1; i <= lim; ++i)
for(int j = 1; i * j <= lim; ++j)
++f[i * j];
for(int i = 1; i <= n; ++i) {
LL tmp = 0;
for(int j = i; j <= n; j += i)
tmp += f[j];
g[i] = tmp % mod;
}
for(int i = 1; i <= m; ++i) {
LL tmp = 0;
for(int j = i; j <= m; j += i)
tmp += f[j];
h[i] = tmp % mod;
}
mu[1] = 1;
for(int i = 2; i <= lim; ++i) {
if(!vis[i]) mu[i] = -1, pri[++pct] = i;
for(int j = 1; j <= pct && i * pri[j] <= lim; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
mu[i * pri[j]] = 0;
break;
} else mu[i * pri[j]] = -mu[i];
}
}
}
signed main() {
cin >> n >> m >> mod;
//n = 2e6, m = 2e6, mod = 998244353;//264180408
if(n > m) n ^= m ^= n ^= m;
Sieve();
rep(d, 1, n) {
int res = 0;
rep(x, 1, n / d) {
if(mu[x] == 1)
res = (res + 1ll * g[d * x] * h[d * x]) % mod;
else if(mu[x] == -1)
res = (res + mod - 1ll * g[d * x] * h[d * x] % mod) % mod;
}
ans = (ans + 1ll * f[d] * res) % mod;
}
cout << ans << '\n';
}
然后发现提交记录有人不到 100ms，就去看了下。好 nb 啊！接下去才是重点。
这题是可以 \(O(n\log\log n)\) 的，一部分一部分看。
首先约数个数和线性筛，多维护每个数最小质因子的次数，降为 \(O(n)\)。
\(g,h\) 可以狄利克雷后缀和 \(O(n\log\log n)\)。
统计答案的式子可以再化，可以看出有个狄利克雷卷积：
\[=\sum_{d=1}^{n}f(d)\sum_{x=1}^{\frac{n}{d}}\mu(x)g(dx)h(dx)\\
=\sum_{i=1}^{n}g(i)h(i)\sum_{j|i}f(j)\mu(\dfrac{i}{j})
\]
看一看 \(d*\mu\) 是什么：\(d=I*I,d*\mu=I*I*\mu=I*\epsilon=I\)。
于是后面那个恒为 \(1\)，答案就是 \(\sum\limits_{i=1}^{n}g(i)h(i)\)。
甚至不用筛 \(\mu\)……
Code
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N = 2000005;
int f[N], g[N], h[N], n, m, mod, ans, sf[N];
int pri[N], pct, low[N];
bool vis[N];
inline int add(int x, int y) { return (x += y) >= mod ? x - mod : x; }
inline void Sieve(const int &lim = N - 1) {
f[1] = 1;
for(int i = 2; i <= lim; ++i) {
if(!vis[i]) low[i] = 1, f[i] = 2, pri[++pct] = i;
for(int j = 1; j <= pct && i * pri[j] <= lim; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
f[i * pri[j]] = f[i] / (low[i] + 1) * (low[i] + 2),
low[i * pri[j]] = low[i] + 1;
break;
} else {
f[i * pri[j]] = f[i] * 2,
low[i * pri[j]] = 1;
}
}
}
rep(i, 1, n) g[i] = f[i];
rep(i, 1, m) h[i] = f[i];
for(int i = 1; i <= lim; ++i) sf[i] = add(sf[i], f[i - 1]);
for(int i = 1; i <= pct; ++i)
for(int j = n / pri[i]; j >= 1; --j)
g[j] = add(g[j], g[j * pri[i]]);
for(int i =1; i <= pct; ++i)
for(int j = m / pri[i]; j >= 1; --j)
h[j] = add(h[j], h[j * pri[i]]);
}
signed main() {
cin >> n >> m >> mod;
// n = 2e6, m = 2e6, mod = 998244353;//264180408
if(n > m) n ^= m ^= n ^= m;
Sieve();
rep(i, 1, n) ans = add(ans, 1ll * g[i] * h[i] % mod);
cout << ans << '\n';
}
\[\sum_{i_1=L}^{H}\sum_{i_2=L}^{H}\cdots\sum_{i_+n=L}^{H}[\gcd(i_1,i_2,\cdots,i_n)=K]\\
\]
记 \(A=\lfloor\dfrac{L-1}{K}\rfloor +1,B=\lfloor\dfrac{R}{K}\rfloor\)，上式可以化简：
\[\sum_{i_1=A}^{B}\sum_{i_2=A}^{B}\cdots\sum_{i_n=A}^{B}[\gcd(i_1,i_2,\cdots,i_n)=1]\\=\sum_{i_1=A}^{B}\sum_{i_2=A}^{B}\cdots\sum_{i_n=A}^{B}\sum_{x|i_1,x|i_2,\cdots,x|i_n}\mu(x)\\=\sum_{x=1}^{B}\mu(x)\sum_{i_1=A}^{B}[x|i_1]\sum_{i_2=A}^{B}[x|i_2]\cdots\sum_{i_n=A}^{B}[x|i_n]\\=\sum_{x=1}^{B}\mu(x)(\lfloor\dfrac{B}{x}\rfloor-\lfloor\dfrac{A-1}{x}\rfloor)^n
\]
后面整除分块前面杜教筛，不懂 \(H-L\le 10^5\) 有什么用……
Code
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
#define mod 1000000007
const int N = 1000005;
int n, K, L, H, ans;
int pri[N], pct, smu[N];
bool vis[N];
inline void Sieve(const int&n = N - 1) {
smu[1] = 1;
for(int i = 2; i <= n; ++i) {
if(!vis[i]) pri[++pct] = i, smu[i] = mod - 1;
for(int j = 1; i * pri[j] <= n; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
smu[i * pri[j]] = 0;
break;
}
smu[i * pri[j]] = mod - smu[i];
}
}
for(int i = 1; i <= n; ++i) smu[i] = (smu[i] + smu[i - 1]) % mod;
}
unordered_map<int, int> _mu;
int sum_mu(int n) {
if(n < N) return smu[n];
if(_mu.find(n) != _mu.end()) return _mu[n];
int res = 1;
for(int l = 2, r; l <= n; l = r + 1)
r = n / (n / l), res = (res + mod - 1ll * (r - l + 1) * sum_mu(n / l) % mod) % mod;
return _mu[n] = res;
}
inline int qpow(int n, int k) {
int res = 1;
for(; k; k >>= 1, n = 1ll * n * n % mod)
if(k & 1) res = 1ll * n * res % mod;
return res;
}
signed main() {
Sieve();
cin >> n >> K >> L >> H, --L;
int A = L / K, B = H / K;
for(int l = 1, r; l <= B; l = r + 1) {
r = B / (B / l);
if(A >= l) ckmin(r, A / (A / l));
ans = (ans + 1ll * qpow(B / l - A / l, n) * (sum_mu(r) - sum_mu(l - 1) + mod)) % mod;
}
cout << ans << '\n';
}
挺不错的一道题，然而我没做出来/kk，体现了我只会做板子题的事实。
\[\dfrac{1}{a}+\dfrac{1}{b}=\dfrac{1}{c}\\
bc+ac=ab\\
ab-bc-ac+c^2=c^2\\
(a-c)(b-c)=c^2
\]
因为 \(\gcd(a,b,c)=1\)，所以 \(\gcd(a-c,b-c,c)=1\)。
然后就是非常神奇的一步，我就卡在这里了：
\[\gcd(a-c,b-c,c)=1\Rightarrow \gcd(a-c,b-c)=1
\]
原因：如果存在一个数 \(x\) 满足 \(x|a-c,x|b-c\)，那么根据 \((a-c)(b-c)=c^2\)，左边必然有因子 \(x^2\)，那么就有 \(x|c\)。
由于 \(\gcd(a-c,b-c)=1\)，那么 \(c^2\) 的每一种质因子只能全给 \(a-c\) 或 \(b-c\)。所以 \(a-c,b-c\) 都是完全平方数。
设 \(x^2=a-c,y^2=b-c\)，那么 \(c=xy\)。
因为 \(\gcd(a-c,b-c)=1\)，所以 \(\gcd(x,y)=1\)。
所以我们可以枚举 \(x\)，这样枚举的范围只有 \(\sqrt{n}\)。记 \(m=\sqrt{n}\)
式子是：
\[\sum_{x=1}^{m}\sum_{y\ge 1}[\gcd(x,y)=1][1\le a\le n][1\le b\le n][1\le c\le n]\\=\sum_{x=1}^{m}\sum_{y\ge 1}[\gcd(x,y)=1][1\le x^2+xy\le n][1\le y^2+xy\le n][1\le xy\le n]
\]
对 \(c\) 的限制显然弱于对 \(a,b\) 的限制，直接去掉：
\[=\sum_{x=1}^{m}\sum_{y\ge 1}[\gcd(x,y)=1][1\le x^2+xy\le n][1\le y^2+xy\le n]
\]
发现后两个限制只是限制了 \(y\) 的范围，不妨直接解出来。
\[x^2+xy\le n\Rightarrow y\le\dfrac{n-x^2}{x}\\y^2+xy\le n\Leftrightarrow y^2+xy-n\le 0\Leftrightarrow \\\dfrac{-x-\sqrt{x^2+4n}}{2}\le y\le \dfrac{-x+\sqrt{x^2+4n}}{2}
\]
令 \(a_x=\min(\dfrac{n-x^2}{x},\dfrac{-x+\sqrt{x^2+4n}}{2})\)，那么答案就是
\[\sum_{i=1}^{m}\sum_{j=1}^{a_i}[\gcd(i,j)=1]
\]
一个人口普查形式的莫反。不过注意不能根号。算了还是写到底吧。
\[\sum_{i=1}^{m}\sum_{j=1}^{a_i}[\gcd(i,j)=1]\\=\sum_{i=1}^{m}\sum_{x|i}\mu(x)\sum_{j=1}^{a_i}[x|j]\\=\sum_{x=1}^{m}\mu(x)\sum_{x|i}\dfrac{a_i}{x}
\]
可以调和级数计算。
Code
#include<bits/stdc++.h>
using namespace std;
#define fi first
#define se second
#define mkp make_pair
#define pb push_back
#define sz(v) (int)(v).size()
typedef long long LL;
typedef double db;
template<class T>bool ckmax(T&x,T y){return x<y?x=y,1:0;}
template<class T>bool ckmin(T&x,T y){return x>y?x=y,1:0;}
#define rep(i,x,y) for(int i=x,i##end=y;i<=i##end;++i)
#define per(i,x,y) for(int i=x,i##end=y;i>=i##end;--i)
inline int read(){
int x=0,f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=0;ch=getchar();}
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
return f?x:-x;
}
const int N = 1000005;
int mu[N], pct, pri[N];
bool vis[N];
LL n, ans, a[N];
int m;
inline void Sieve(const int &n = N - 1) {
mu[1] = 1;
for(int i = 2; i <= n; ++i) {
if(!vis[i]) pri[++pct] = i, mu[i] = -1;
for(int j = 1; j <= pct && i * pri[j] <= n; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
mu[i * pri[j]] = 0;
break;
} else mu[i * pri[j]] = -mu[i];
}
}
}
signed main() {
cin >> n, m = sqrt(n);
Sieve(m);
rep(i, 1, m) a[i] = min((n - 1ll * i * i) / i, (LL)(sqrt(1ll * i * i + 4ll * n) - i) / 2);
rep(i, 1, m) {
LL tmp = 0;
for(int j = 1; i * j <= m; ++j) tmp += 1ll * a[i * j] / i;
ans += tmp * mu[i];
}
cout << ans << '\n';
}
好像是曾经的月赛题，场上没做出来，现在还是不会做……
考虑怎么化简这个 \(p(a,b)\)。
首先注意到，\(\gcd(a,b)>1\) 就不行。又注意到 \(a\bmod 2 = b\bmod 2\) 也不行。剩下的凭直觉就感觉都可以啊，打个表感觉很对啊，往上一交就 20 了啊。
也就是这个东西：
inline int p(int x, int y) {
return ((x & 1) ^ (y & 1)) && (__gcd(x, y) == 1);
}
所以答案就是
\[\sum_{i=1}^{n}\sum_{j=1}^{n}[\gcd(i,j)=1][i\bmod 2 \not= j\bmod 2]
\]
不妨钦定 \(i\) 为偶数 \(j\) 为奇数，问题转化成 \(\sum_{i=1}^{\frac{n}{2}}\sum_{j=1}^{n}[\gcd(i,j)=1]\)，后面那个东西可以根号，人口普查形式的莫反不再多说，复杂度 \(O(n\sqrt{n})\) 得到 35。
感觉这个思路没得救，换个思路。
考虑计算 \(f(x)=\sum_{i=1}^{x}[i\bmod 2 \not= x \bmod 2][\gcd(x,i)=1]\)。
对于 \(x\bmod 2 = 0\)，直接 \(f(x)=\varphi(x)\) 即可，因为所有偶数都不会被统计在欧拉函数内。
对于 \(x\bmod 2 = 1\)，\(f(x)=\dfrac{\varphi(x)}{2}\)，因为欧拉函数每统计到一个与 \(x\) 互质的奇数数 \(y\)，都会再统计到一个与 \(x\) 互质的偶数 \(x-y\)，所以除以二。
特判 \(x=1\)，最后答案乘 \(2\) 即可，线性筛 \(\varphi\)，可以得到 50 分。
看看现在的式子（把 \(2\) 先乘进去）：
\[\sum_{i=1}^{n}\left([i\bmod 2 = 0]2\varphi(i)+[i\bmod 2 = 1]\varphi(i)\right)-1
\]
设 \(f(n)=\sum_{i=1}^{n}[i\bmod 2 = 0]2\varphi(i)+[i\bmod 2 = 1]\varphi(i)\)
注意到可以变形成 \(f(n)=\sum_{i=1}^{n}\varphi(i)+\sum_{i=1}^{n}[i\bmod 2 = 0]\varphi(i)\)。
而 \(\sum_{i=1}^{n}[i\bmod 2 = 0]\varphi(i)\)，可以发现它就等于 \(f(\lfloor\dfrac{n}{2}\rfloor)\)！
具体分析就是，如果 \(i\bmod 2 = 0\)，那么 \(\varphi(2i)=2\varphi(i)\)，否则 \(\varphi(2i)=\varphi(i)\)。分讨一下 \(\le \lfloor\dfrac{n}{2}\rfloor\) 的奇偶性，发现刚好就是。
所以 \(f(n)=\sum_{i=1}^{n}\varphi(i)+f(\lfloor\dfrac{n}{2}\rfloor)\)，直接递归，加个杜教筛，就好了！
无限 orz 出题人 QuantAsk！！！
Code
#include <cstdio>
#include <tr1/unordered_map>
using namespace std::tr1;
const int N = 10000005;
typedef unsigned long long ull;
int pri[N], pct;
ull phi[N], g[N];
bool vis[N];
inline void Sieve(const int &n = N - 1) {
phi[1] = 1;
for(int i = 2; i <= n; ++i) {
if(!vis[i]) pri[++pct] = i, phi[i] = i - 1;
for(int j = 1; j <= pct && i * pri[j] <= n; ++j) {
vis[i * pri[j]] = 1;
if(i % pri[j] == 0) {
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
phi[i * pri[j]] = phi[i] * (pri[j] - 1);
}
}
for(int i = 1; i <= n; ++i)
g[i] = g[i - 1] + (i & 1 ? phi[i] : phi[i] + phi[i]), phi[i] += phi[i - 1];
}
unordered_map<ull, ull> _phi;
ull sum_phi(ull n) {
if(n < N) return phi[n];
if(_phi.find(n) != _phi.end()) return _phi[n];
ull res = n & 1 ? (n + 1) / 2 * n : n / 2 * (n + 1);
for(ull l = 2, r; l <= n; l = r + 1)
r = n / (n / l), res -= (r - l + 1) * sum_phi(n / l);
return _phi[n] = res;
}
ull f(ull n) {
if(n < N) return g[n];
return f(n / 2) + sum_phi(n);
}
signed main() {
Sieve();
int T; ull n;
for(scanf("%d", &T); T--;)
scanf("%lld", &n), printf("%llu\n", f(n) - 1);
}
路漫漫其修远兮,吾将上下而求索
