Du‘s Sieve

由于是没有质量的垃圾博客，就直接放 csdn 了。

对于一个积性函数 $f(n)$，我们的目标是求出 $S(n)=\sum f(i)$，我们选择构造积性函数 $g(n)$，并以 $g(1)S(n)$ 的形式间接求出 $S(n)$（内心 p.s. 这么 constructive 的东西怎么发明的）。

我们把 Dirichlet Convolution 的形式略改一下，$\sum _{i=1}^n(f\ast g)(i)=\sum _{i=1}^n\sum _{d\mid i}f(d)\times g(\frac{i}{d})=\sum _{d=1}^{n}f(d)\sum _{i=1}^{n/d}g(i)=$，带入 $g$ 和 $S$ 得 $\sum _{i=1}^n\sum _{d\mid i}g(d)f\left(\frac{i}{d}\right)=\sum _{i=1}^{n}g(i)S(n/i)$。

于是 $g(1)S(n)=\left(\sum _{i=1}^n(f\ast g)(i)\right)-\left(\sum _{i=2}^{n}g(i)S(n/i)\right)$。

直接看题捏。

51nod 1227 平均最小公倍数

给定 $a,b$，求 $\sum _{i=a}^b\sum _{j=1}^i\frac{j}{(i,j)}$
$1⩽a,b⩽10^9$

考虑前缀答案，容易化成 $\sum _{d=1}^n\sum _{i=1}^{n/d}\sum _{j=1}^{i}j[(i,j)=1]$，于是我们考虑求 $f(n)$ = 与小于 $n$ 且与其互质的数之和。

于是有，$f(n)=\sum _{i=1}^{n}i[(i,n)=1]=\sum _{i=1}^{n}i\sum _{k\mid (i,n)}\mu (k)=\sum _{k\mid n}\mu (k)\times k\sum _{i=1}^{n/k}i=\frac{n}{2}\left(\left(\sum _{k\mid n}\mu (k)\times n/k\right)+\left(\sum _{k\mid n}\mu (k)\right)\right)=\frac{n}{2}(\phi (n)+ϵ(n))$。

于是原式可化为 $\frac{n}{2}+\frac{1}{2}\left(\sum _{d=1}^n\sum _{i=1}^{n/d}i\times \phi (i)\right)$，观察到 $\text{id}(i)\ast i\times \phi (i)=\sum _{d\mid i}\phi (i)\times i/d=i^2$，套一个二次方和公式即可杜教筛了。logical & reasonable！

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define cmin(x, ...) x = min({x, __VA_ARGS__})
#define cmax(x, ...) x = max({x, __VA_ARGS__})
#define fors(i, l, r, ...) for(int i = (l), REP##i = (r), ##__VA_ARGS__; i <= REP##i; ++i)
#define dfors(i, r, l, ...) for(int i = (r), REP##i = (l), ##__VA_ARGS__; i >= REP##i; --i)
const int P = 1e9+7;
struct custom_hash {
    static uint64_t splitmix64(uint64_t x) {
        x += 0x9e3779b97f4a7c15;
        x = (x^(x>>30))*0xbf58476d1ce4e5b9;
        x = (x^(x>>27))*0x94d049bb133111eb;
        return x^(x>>31);
    }
    size_t operator()(uint64_t x) const {
        static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
        return splitmix64(x+FIXED_RANDOM);
    }
};
signed main() {
    ios::sync_with_stdio(0), cin.tie(0);
    // f(n) = sum(1~n) i*phi(i)
    // id*f = sum(1~n) i^2
    // S(n) = sum(1~n) i^2 - sum(2~n) S(n/i)*i
    const function<int(const ll)> p_f = [&](const ll x) {
        static const ll THRESHOLD = 1e7;
        static unordered_map<ll, int, custom_hash> g_f;
        static auto l_f = [&](const int n) {
            vector<int> f(n+1), prime;
            vector<bool> not_prime(n+1);
            f[1] = 1;
            fors(i, 2, n) {
                if(not_prime[i] == 0) f[i] = i-1,prime.emplace_back(i);
                for(const int p : prime) {
                    if(i > n/p) break;
                    not_prime[i*p] = 1;
                    if(i%p == 0) {
                        f[i*p] = ll(f[i])*p%P;
                        break;
                    }
                    f[i*p] = ll(f[i])*f[p]%P;
                }
            }
            fors(i, 1, n) f[i] = (f[i-1]+ll(f[i])*i%P)%P;
            return f;
        }(THRESHOLD);
        if(x <= THRESHOLD) return l_f[x];
        if(g_f.count(x)) return g_f[x];
        int _x = x%P, res = ll(_x)*((2*_x+1)%P)%P*(_x+1)%P*166666668%P; // sum(1~n) i^2
        for(ll l = 2, r; l <= x; l = r+1) {
            r = x/(x/l);
            res -= (ll)p_f(x/l)*(r-l+1)%P*(l+r)%P*(P+1)/2%P, res %= P;
        }
        return g_f[x] = (res+P)%P;
    };
    auto cal = [&](const ll x) {
        int res = x;
        for(ll l = 1, r; l <= x; l = r+1) {
            r = x/(x/l);
            res += (ll)p_f(x/l)*(r-l+1)%P, res %= P;
        }
        (res += P) %= P;
        return ll(res)*(P+1)/2%P;
    };
    ll a, b; cin >> a >> b;
    cout << (cal(b)-cal(a-1)+P)%P << "\n";
    return 0;
}