Grisaia (推公式) 2018 第十届四川省程序设计竞赛_2018四川省大学生程序设计题解-优快云博客

本文链接：https://blog.youkuaiyun.com/KXL5180/article/details/99073728

本文深入解析一道数论竞赛题，通过巧妙的数学推导，将复杂的问题简化为可计算的公式。文章详细介绍了如何利用积性函数和线性筛法优化计算过程，最终实现了高效求解。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

https://www.oj.swust.edu.cn/problem/show/2810

题意：很简单就是求题目上面那个式子的和。

做法：直接开始推公式。

$ans=\sum_{i=1}^{n}\sum_{j=1}^{i}(n\ mod (i \times j))=\sum_{i=1}^{n}\sum_{j=1}^{i}(n-\left \lfloor \frac{n}{ij} \right \rfloor \times ij)$

然后，将前后两部分分开

$ans=\sum_{i=1}^{n}\sum_{j=1}^{i}n-\sum_{i=1}^{n}\sum_{j=1}^{i}\left \lfloor \frac{n}{ij} \right \rfloor \times ij$

然后对于前面的一部分直接可以由公式计算得到： $\sum_{i=1}^{n}\sum_{j=1}^{i}n=\frac{n^{2}(n+1)}{2}$

我们现在考虑后面一部分怎么计算。我刚刚拿到的时候一点头绪都没有。

但不过我们做题做多了，经常可以发现这种式子一般可以化成，两个上界为n的和式：

$last1=\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}\left \lfloor \frac{n}{ij} \right \rfloor \times ij$

这样是对的吗，很可惜，我打了一下表发现并不对。。。。。

发现是漏掉一些情况，其实我们可以这样想一个对称矩阵，你现在想用一个矩阵的所有元素的值，加起来求下三角的值

这是你会发现对角线的值会有问题，这时你会把对角线的值乘上2，在平方。自己yy一下吧。

正确公式这样的 $\frac{1}{2}(\sum_{i}^{n}\left \lfloor \frac{n}{i^{2}} \right \rfloor \times i^{2}+\sum_{i=1}^{n}\sum_{j=1}^{n}\left \lfloor \frac{n}{ij} \right \rfloor \times ij)$

~~其实我是打表结合瞎猜的~~

关于前面一部分可以直接 $O(\sqrt n)$ ，求出来。

我们重点讨论后面一部分怎么做，打了一下表发现根本没有规律可循，好像也不可以莫比乌斯和欧拉。。。

还是首先先进行变形看看吧(其实是乱搞)：

$\sum_{i=1}^{n}\sum_{j=1}^{n}\left \lfloor \frac{n}{ij} \right \rfloor \times ij=\sum_{i=1}^{n}\sum_{j=1}^{\left \lfloor \frac{n}{i} \right \rfloor} \left \lfloor \frac{\left \lfloor \frac{n}{i} \right \rfloor}{j} \right \rfloor \times ij$ 我觉得这一步有一点点不好想到，弄了好久，这样i和j就可以分开。

我们令： $h(x)=\sum_{i=1}^{x} \left \lfloor \frac{x}{i} \right \rfloor \times i$

然后发现原式可以这样变化：

$\sum_{i=1}^{n}\sum_{j=1}^{\left \lfloor \frac{n}{i} \right \rfloor} \left \lfloor \frac{\left \lfloor \frac{n}{i} \right \rfloor}{j} \right \rfloor \times ij=\sum_{i=1}^{n}i\times h(\left \lfloor \frac{n}{i} \right \rfloor)$

这样以及很不错了以及可以分块求这一部分了，但是，如果直接这样前面的i虽然可以用公式，而后面的则是 $O(\sum \sqrt {\left \lfloor \frac{n}{i} \right \rfloor})$

这个复杂度，虽然到后面会慢慢的降低，但不过前面的复杂的会很大，已经T了。

然后我们观察 $i\times h(\left \lfloor \frac{n}{i} \right \rfloor)$ ，或者 h(x) ，好像不能筛吧。。。看看有没有什么性质，比如说什么积性函数啊，用杜教筛等等。很遗憾打了一天表毛关系都没有找到。

我们还是来进行一些稍微严格~~(乱搞)~~的推导吧：

我们令： g(n)=h(n)-h(n-1) ，我们发现g(n)的前缀和就是h(n)

$g(n)=\sum_{i=1}^{n}\left \lfloor \frac{n}{i} \right \rfloor \times i - \sum_{i=1}^{n-1}\left \lfloor \frac{n-1}{i} \right \rfloor \times i=\sum_{i=1}^{n}(\left \lfloor \frac{n}{i} \right \rfloor -\left \lfloor \frac{n-1}{i} \right \rfloor) \times i$ 这个式子除了i等于n的约数情况下其他全部为零。

并且等于约数情况下，就等于i 所以就是个约数和： $g(n)=\sum_{d|n}d$

~~不信？你可以打打表。~~

这肯定是一个积性函数，然后我们就可以线性筛了。预处理，他的前面的一部分的和。这样时间复杂度就可以得到优化了。

其实，我真正的经历是这样的：确实弄了好久发现h(n)的差是一个积性函数，但不过并不知道是约数和，然后就在哪里线性筛。

筛了大半天，才筛出来和打表一样，然后想着进行狄利克雷卷积，用杜教筛，然后卷了大半天，终于才发现就是一个约数和。

然后就百度一下，没有找到。然后我寻思这一想，直接对剩下的h(n)进行分块暴力搞，这样复杂度也不高吧，好像和杜教筛差不多耶。

~~我记住了约数和的前缀和可以这样求了。。。~~

~~我就接着上面那样随便搞了一下就过了。。。。~~

我们来分析一下复杂度~~(瞎说)~~： $O(1+\sqrt{n} + p+q\sqrt{\left \lfloor \frac{n}{i} \right \rfloor})$

其实后面一部分的p是预处理的大小，q也是和p有关系是这样的 $\frac{n}{p}$ ，我们近似处理一下：

$O(p+\frac{n^{3/2}}{p})$ ,我们根据均值不等式，得出了是n的0.75次方差不多最合适，但不过由于内存，加上答案会爆long long 我们用int128来处理，会爆内存，因此我们开2000000万差不多了。

~~题解貌似是n^3/2...~~

#include "bits/stdc++.h"

using namespace std;
///typedef long long ll;
typedef __int128 ll;
const int N = 20000000 + 10;
int vis[N], cnt = 0;
ll f[N], h[N];
long long n, pri[N], num[N];

void print(__int128 x) {
    if (!x) {
        puts("0");
        return;
    }
    string ret = "";
    while (x) {
        ret += x % 10 + '0';
        x /= 10;
    }
    reverse(ret.begin(), ret.end());
    cout << ret << endl;
}

void init() {
    f[1] = vis[1] = 1;
    for (int i = 2; i < N; i++) {
        if (!vis[i]) {
            pri[++cnt] = i;
            f[i] = i + 1;
            num[i] = i;
        }
        for (int j = 1; j <= cnt && i * pri[j] < N; j++) {
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0) {
                if (i / num[i] != 1) f[i * pri[j]] = f[num[i] * pri[j]] * f[i / num[i]];
                else f[i * pri[j]] = 1LL * i * pri[j] + f[i];
                num[i * pri[j]] = num[i] * pri[j];
                break;
            }
            f[i * pri[j]] = f[i] * f[pri[j]];
            num[i * pri[j]] = pri[j];
        }
    }
    for (int i = 1; i < N; i++) {
        h[i] = h[i - 1] + f[i];
    }
}

ll get_last(ll n) {
    ll ret = 0;
    for (ll i = 1; i * i <= n; i++) {
        ret += (n / (i * i)) * i * i;
    }
    return ret;
}

ll get_s1(ll x) {
    ll ret = x * (x + 1) / 2;
    return ret;
}

ll fun(ll x) {
    if (x < N) return h[x];
    ll ret = 0;
    for (ll l = 1, r; l <= x; l = r + 1) {
        r = x / (x / l);
        ret += (x / l) * (get_s1(r) - get_s1(l - 1));
    }
    return ret;
}

ll solve(ll n) {
    ll ans = n * n * (n + 1) / 2;
    ll pre = 0;
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        pre += (get_s1(r) - get_s1(l - 1)) * fun(n / l);
    }
    ans = ans - (pre + get_last(n)) / 2;
    return ans;
}

int main() {
    int T;
    init();
    scanf("%d", &T);
    while (T--) {
        scanf("%lld", &n);
        print(solve(n));
    }
    return 0;
}