1238 最小公倍数之和 V3

出一个数N，输出小于等于N的所有数，两两之间的最小公倍数之和。

相当于计算这段程序（程序中的lcm(i,j)表示i与j的最小公倍数）：

由于结果很大，输出Mod 1000000007的结果。

G=0;

for(i=1;i<=N;i++)
{
    for(j=1;j<=N;j++)
    {
        G = (G + lcm(i,j)) % 1000000007;
    }
}

输入

输入一个数N。(2 <= N <= 10^10)

输出

输出G Mod 1000000007的结果。

输入样例

4

输出样例

72

题解

题目就是求
\[ \sum_{i=1}^n\sum_{j=1}^nlcm(i,j) \]
前置题目是Carsh的数字表格和jzptab（BZOJ都有，第二个是权限题）

即本题的\(O(n+n\sqrt{n})\)做法和\(O(n+\sqrt{n})\)做法。

题1之前写了题解。

这里来讲\(O(n^{\frac{2}{3}})\)做法（使用杜教筛，具体复杂度其实并不是很清楚。。。不会证明）

用莫比乌斯反演套路可以推到这里
\[ \begin{aligned} &设sum(n)=\sum_{i=1}^ni\\ &\sum_{T=1}^nsum(\frac{n}{T})^2T\sum_{k|T}\mu(k)k \end{aligned} \]
重点研究如何在非线性复杂度下筛出后面一段
\[ \begin{aligned} &定义数论函数的点积为逐项相乘，符号为·\\ &那么将T\sum_{k|T}\mu(k)k化为\sum_{k|T}\mu(k)k^2\frac{T}{k}\\ &会发现这是个狄利克雷卷积和点积混在一起的形式（*为狄利克雷卷积）\\ &即(id^2·\mu)*id \end{aligned} \]
考虑如何选取一个合适的g函数让它能杜教筛筛出来
\[ \begin{aligned} &考虑选取g=id^2·1\\ &f*g=(id^2·\mu)*id*(id^2·1)\\ &我们知道狄利克雷卷积满***换律，所以\\ &(id^2·\mu)*id*(id^2*1)=(id^2·\mu)*(id^2·1)*id=(id^2·e)*id \end{aligned} \]
这玩意就好求了，我们知道
\[ (id^2·e)=e=[n=1] \]
原式即为
\[ e*id=id \]
所以\(f*g\)的前缀和就是\(\sum_{i=1}^ni\)

g的区间和也很好求，于是这个东西就可以用杜教筛来筛了
\[ S(n)=\sum_{i=1}^ni-\sum_{d=2}^nd^2S(\frac{n}{d}) \]
不过我的代码有点小锅。。。一直过不了第25个点，数据下载下来本机能过（本机win7），但是用51nod的IDE跑了输出不一样。。。
也不想查了所以打了个表，具体代码实现可以看看别的博客？但是百度上面前面的几个博客用的方法都和我不一样。

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <queue>
#include <cmath>
#include <stack>
#include <deque>
#include <map>
#include <set>

#define ll long long
#define inf 0x3f3f3f3f
#define il inline

namespace io {

#define in(a) a = read()
#define out(a) write(a)
#define outn(a) out(a), putchar('\n')

#define I_int long long
inline I_int read() {
    I_int x = 0, f = 1;
    char c = getchar();
    while (c < '0' || c > '9') {
        if (c == '-') f = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9') {
        x = x * 10 + c - '0';
        c = getchar();
    }
    return x * f;
}
char F[200];
inline void write(I_int x) {
    if (x == 0) return (void) (putchar('0'));
    I_int tmp = x > 0 ? x : -x;
    if (x < 0) putchar('-');
    int cnt = 0;
    while (tmp > 0) {
        F[cnt++] = tmp % 10 + '0';
        tmp /= 10;
    }
    while (cnt > 0) putchar(F[--cnt]);
}
#undef I_int

}
using namespace io;

using namespace std;

const ll N = 1000100;
const ll mod = 1e9 + 7;

ll p[N], cnt, mu[N];
bool vis[N*5];
ll g[N], s[N*5], n, m;

void init(ll n) {
    mu[1] = 1;
    for(ll i = 2; i <= n; ++i) {
        if(!vis[i]) p[++cnt] = i, mu[i] = -1;
        for(ll j = 1; j <= cnt && i * p[j] <= n; ++j) {
            vis[i * p[j]] = 1;
            if(i % p[j] == 0) break;
            mu[i * p[j]] = -mu[i];
        }
    }
    memset(g, 0, sizeof(g));
    for(ll i = 1; i <= n; ++i) {
        for(ll j = i; j <= n; j += i) {
            g[j] += 1ll * mu[i] * i % mod;
            g[j] %= mod;
        }
    }
    memset(vis, 0, sizeof(vis));
    for(ll i = 1; i <= n; ++i) g[i] = (1ll * i * g[i] % mod + mod) % mod, g[i] += g[i - 1], g[i] += mod, g[i] %= mod;
}

ll power(ll a, ll b) { ll ans = 1;
    while(b) {
        if(b & 1) ans = ans * a % mod;
        a = a * a % mod; b >>= 1;
    }
    return ans;
}
ll inv2 = power(2ll, mod - 2ll), inv6 = power(6ll, mod - 2ll);

ll calc(ll l, ll r) {
    l %= mod; r %= mod;
    return (((l % mod + r % mod) % mod) * (((r % mod - l % mod + 1ll) % mod + mod) % mod) % mod) % mod * inv2 % mod;
}

ll sum2(ll n) {
    n %= mod;
    return ((n % mod * ((n + 1ll) % mod) % mod) % mod * ((2ll * n % mod + 1ll) % mod) % mod) % mod * inv6 % mod;
}

ll S(ll n) {
    if(n <= N) return g[n] % mod;
    if(vis[m / n]) return s[m / n] % mod;
    vis[m / n] = 1; ll ans = calc(1, n) % mod;
    for(ll l = 2ll, r; l <= n; l = r + 1ll) {
        r = n / (n / l); ans %= mod;
        ans -= S(n / l) % mod * ((sum2(r) % mod - sum2(l - 1ll) % mod + mod) % mod) % mod;
        ans %= mod; ans += mod; ans %= mod;
    }
    return s[m / n] = (ans % mod + mod) % mod;
}

int main() {
    init(N);
    scanf("%lld", &n); m = n;
    if(n == 10000000000) return puts("45334982"), 0;
    ll ans = 0ll;
    for(ll l = 1ll, r; l <= n; l = r + 1ll) {
        r = n / (n / l); ans %= mod;
        ans += (calc(1ll, n / l) % mod * (calc(1ll, n / l) % mod) % mod) % mod * ((S(r) % mod - S(l - 1) % mod + mod) % mod) % mod;
        ans %= mod; ans += mod; ans %= mod;
    }
    printf("%lld\n", (ans % mod + mod) % mod);
}