ra (数论 , 莫比乌斯反演 , 整点统计)
题意
求 $$\displaystyle \sum_{i=1}^{n} \sum_{j=1}^{n} [\mathrm{lcm} (i,j) > n] \pmod {10^9 + 7}$$ .
$ n \le 10^{10}$ .
题解
这是我们考试的一道题 ... 考试的时候以为能找出规律 , 后来发现还是一道数论题 qwq
而且部分分很不良心啊 , 只给了 \(O(n)\) 多的一点分 , 我 \(O(n \ln n)\) 根本没活路 .. 还是直接开始推吧 ~
\[\begin{align}
\sum_{i=1}^{n} \sum_{j=1}^{n} [\mathrm{lcm} (i,j) > n] &=n^2- \sum_{i=1}^{n} \sum_{j=1}^{n} [\mathrm{lcm} (i,j) \le n]\\
&= n^2 - \sum_{d=1}^{n} \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} [ijd \le n] \cdot [i \bot j] \\
&= n^2 - \sum_{d=1}^{n} \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}[ijd \le n] \sum_{x|(i,j)}\mu(x) \\
&= n^2 - \sum_{d=1}^{n} \sum_{x=1}^{\lfloor \frac{n}{d}\rfloor} \mu(x) \sum_{i=1}^{\lfloor \frac{n}{dx} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{dx} \rfloor}[ijdx^2 \le n] \\
&= n^2 - \sum_{x=1}^{n} \mu(x) \sum_{d=1}^{\lfloor \frac{n}{x} \rfloor} \sum_{i=1}^{\lfloor \frac{n}{dx} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{dx} \rfloor}[ijdx^2 \le n] \\
\end{align}
\]
到这一步不难发现由于 \([ijdx^2 \le n]\) 可以缩减很多范围了 比如 \(x \le \lfloor \sqrt n \rfloor\) ... 直接一波缩范围
\[\displaystyle = n^2 - \sum_{x=1}^{\lfloor \sqrt n \rfloor} \mu(x) \sum_{d=1}^{\lfloor \frac{n}{x^2} \rfloor} \sum_{i=1}^{\lfloor \frac{n}{dx^2} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{dx^2} \rfloor}[ij \le \lfloor \frac{n}{dx^2} \rfloor]
\]
我们可以考虑看一下后面两个 \(\sum\) 好像很有特点。令
\[\displaystyle f(x) = \sum_{i=1}^{x} \sum_{j=1}^{x} [ij \le x]
\]
那么原式就是
\[\displaystyle = n^2 - \sum_{x=1}^{\lfloor \sqrt n \rfloor} \mu(x) \sum_{d=1}^{\lfloor \frac{n}{x^2} \rfloor} f(\lfloor \frac{n}{dx^2} \rfloor)
\]
观察一下 \(f(x)\) 好像也可以进行转化
考虑枚举一个 \(i,j\) 的积 , 看有多少对 \((i,j)\) 可以 .
\[\displaystyle f(x) = \sum_{d=1}^{x} \lfloor \frac{x}{d} \rfloor
\]
这个容易在 \(O(\sqrt n)\) 直接分块解决 . 这样带入直接做就有 60 分了 \((n \le 10^8)\) , 不会积分证明复杂度QAQ ....
卡一卡 , 在本机上能跑 \(10^9\) 能拿80分 爽歪歪 qwq
后来我意识到瓶颈在 \(f(x)\) 处 , 各种问人是否有公式计算 .... 后来才发现 这个竟然是今年集训队论文 ??!!!
**《一些特殊的数论函数求和问题》 —— 安徽师范大学附属中学 朱震霆 **
考虑我最初的那个式子
\[\displaystyle f(n) = \sum_{i=1}^{n} \sum_{j=1}^{n} [ij \le n]
\]
难道不就是数 \(xy=n\) 下面的整点个数吗 !!
我认真看了论文许久，可还是看不懂，只知道大概就是用很多根切线去分割，然后去数切线下方的点.
过几天看懂了再来理解 .... 只知道复杂度是 \(O(n^{\frac{1}{3}})\) 的，十分优秀 ~
然后直接找到 whzzt 的代码 ，尝试着放进去我的程序...
竟然过了！！！跑了 \(0.5s\) 就过了.... (原来要跑 \(6s\) )
挂一波代码就跑 qwq
代码
#include <bits/stdc++.h>
#define For(i, l, r) for (register ll i = (ll)(l), i##end = (ll)(r); i <= i##end; ++ i)
#define Fordown(i, r, l) for (register ll i = (ll)(r), i##end = (ll)(l); i >= i##end; -- i)
#define Set(a, v) memset(a, v, sizeof(a))
using namespace std;
typedef long long ll;
inline bool chkmin(ll &a, ll b) { return b < a ? a = b, 1 : 0; }
inline bool chkmax(ll &a, ll b) { return b > a ? a = b, 1 : 0; }
inline ll read() {
ll x = 0, fh = 1; char ch = getchar();
for (; !isdigit(ch); ch = getchar()) if (ch == '-') fh = -1;
for (; isdigit(ch); ch = getchar()) x = (x * 10) + (ch ^ 48);
return x * fh;
}
void File() {
freopen ("ra.in", "r", stdin);
freopen ("ra.out", "w", stdout);
}
const ll Mod = 1e9 + 7;
const ll N = 1e6 + 1e3;
ll mu[N], prime[N], cnt = 0; bitset<N> is_prime;
void Init(ll maxn) {
is_prime.set(); is_prime[0] = is_prime[1] = false; mu[1]= 1;
For (i, 2, maxn) {
if (is_prime[i])
prime[++ cnt] = i, mu[i] = -1;
For (j, 1, cnt) {
ll res = prime[j] * i;
if (res > maxn) break ;
is_prime[res] = false;
if (i % prime[j]) mu[res] = - mu[i];
else { mu[res] = 0; break ; }
}
}
}
/*inline ll SumDown(ll a) {
ll res = M[a]; if (res) return res;
For (i, 1, a) {
register ll now = a / i, Nexti = a / now;
res += now * (Nexti - i + 1); i = Nexti;
}
return (M[a] = res % Mod);
}*/
typedef unsigned long long uLL;
typedef unsigned long long ull;
typedef unsigned int uint;
unordered_map<ull, uLL> M;
namespace ds {
namespace stac {
const int N = 100005;
uint qu[N][2]; int qr;
inline void pop () { qr --; }
inline void push (uint x, uint y) { qr ++; qu[qr][0] = x; qu[qr][1] = y; }
inline void top (uint &x, uint &y) { x = qu[qr][0]; y = qu[qr][1]; }
}
using stac :: push;
using stac :: pop;
using stac :: top;
inline uLL solve (ull n) {
uLL ret = M[n];
if (ret) return ret;
ull w = pow (n, 0.38), v = sqrtl (n), x, y;
uint dx, dy, ux, uy, mx, my;
while (v * v <= n) v ++; while (v * v > n) v --;
x = n / v, y = n / x + 1;
push (1, 0); push (1, 1);
auto outside = [&] (ull x, ull y) { return x * y > n; };
auto cut_off = [&] (ull x, uint dx, uint dy) { return (uLL)x * x * dy >= (uLL)n * dx; };
while (stac :: qr) {
top (dx, dy);
while (outside (x + dx, y - dy)) {
ret += x * dy + ull(dy + 1) * (dx - 1) / 2;
x += dx, y -= dy;
}
if (y <= w) break;
while (true) {
pop (), ux = dx, uy = dy, top (dx, dy);
if (outside (x + dx, y - dy)) break;
}
while (true) {
mx = ux + dx, my = uy + dy;
if (!outside (x + mx, y - my)) {
if (cut_off (x + mx, dx, dy)) break;
ux = mx, uy = my;
} else push (dx = mx, dy = my);
}
}
for (y --; y; y --) ret += n / y;
return stac :: qr = 0, (M[n] = ret * 2 - v * v);
}
}
int main() {
File();
ll n = read(), res = 0;
Init(1e6);
For (x, 1, sqrt(n)) if (mu[x]) {
register ll Lim = n / (x * x), tot = 0;;
For (d, 1, Lim) {
register ll now = Lim / d, Nextd = Lim / now;
tot += ds :: solve(now) * (Nextd - d + 1); d = Nextd;
}
(res += Mod + tot % Mod * mu[x]) %= Mod;
}
res = ((n % Mod) * (n % Mod) % Mod - res + Mod) % Mod;
cout << res << endl;
#ifdef zjp_shadow
cerr << (double) clock() / CLOCKS_PER_SEC << endl;
#endif
return 0;
}
