ra (数论 , 莫比乌斯反演 , 整点统计)

Posted 2020-11-10 zjp_shadow

• 题意 :

求 \[\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)\) 根本没活路 ..

还是直接开始推吧 ~

\[\displaystyle \sum_{i=1}^{n} \sum_{j=1}^{n} [\mathrm{lcm} (i,j) > n]\]

\[=n^2-\displaystyle \sum_{i=1}^{n} \sum_{j=1}^{n} [\mathrm{lcm} (i,j) \le n]\]

\[\displaystyle = 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] \]

\[\displaystyle = 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) \]

\[\displaystyle = n^2 - \sum_{d=1}^{n} \sum_{x=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{i=1}^{\lfloor \frac{n}{dx} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{dx} \rfloor}[ijdx^2 \le n] \]

\[\displaystyle = n^2 - \sum_{x=1}^{n} \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] \]

到这一步不难发现可以缩减很多范围了 ... 直接一波缩范围

\[\displaystyle = n^2 - \sum_{x=1}^{\lfloor \sqrt n \rfloor} \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} \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}^{x} \sum_{j=1}^{x} [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;
}