洛谷P3768 简单的数学题题解

题意：

输入一个整数 \(n\) 和一个整数 \(p\)，你需要求出：

\[ \left(\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)\right) \bmod p \] 其中 \(\gcd(a,b)\) 表示 \(a\) 与 \(b\) 的最大公约数。

输入格式：

一行两个整数 \(p,n\)。

输出格式：

一行一个整数表示答案。

数据范围：

\(1\le n \le 10^{10},~5 \times 10^8 \leq p \leq 1.1 \times 10^9\) 且 \(p\) 为质数。时限 4s。

考虑欧拉反演。 \[ \begin{aligned} \sum_{i=1}^n \sum_{j=1}^n i j \operatorname{gcd}(i, j) & =\sum_{i=1}^n \sum_{j=1}^n i j \sum_{k\mid i,\, k\mid j} \varphi(k) \\[4pt]& =\sum_{k=1}^n \varphi(k) \sum_{k \mid i} \sum_{k \mid j} i j \\[4pt]& =\sum_{k=1}^n \varphi(k) \cdot k^2 \cdot \left(\sum_{i=1}^{\left\lfloor n/k \right\rfloor} i\right)^2 \\[4pt]& =\sum_{k=1}^n k^2\varphi(k) \sum_{i=1}^{\left\lfloor n/k\right\rfloor} i^3 \\[4pt]& =\sum_{k=1}^n k^2\varphi(k) S\left(\left\lfloor n/k\right\rfloor\right) && S(n) = \left(\frac{n(n+1)}{2}\right)^2 \end{aligned} \] 显然 \(S(n)\) 可以数论分块，因此我们只需要求出 \(k^2\varphi(k)\) 的前缀和。

考虑杜教筛。记 \(f=\varphi \cdot \mathrm{id}^2\) ，设 \(g=\mathrm{id}^2\) ，则 \[ \begin{aligned} (f*g)(n) &= \sum_{d \mid n}\varphi(d)d^2\cdot \frac{n^2}{d^2} \\[6pt] &= n^2\sum_{d \mid n}\varphi(d) \\[6pt] &= n^3 \end{aligned} \] 那么前缀和就是一个 \(\sum i^2\) 和 \(\sum i^3\) ，分别是 \(\frac{n(n+1)(2 n+1)}{6}\) 和 \(\left(\frac{n(n+1)}{2}\right)^2\) 。

时间复杂度 \(\mathcal{O}(n^{\frac{2}{3}})\) ，因为杜教筛有记忆化，所以筛一次就能筛出这 \(\mathcal{O}(\sqrt{n})\) 个值。

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x, int y) { x < y ? x = y : 0; }
void down(int &x, int y) { x > y ? x = y : 0; }
#define rep(i, a, b) for(int i = (a), i##END = (b); i <= i##END; i++)
#define Rep(i, a, b) for(int i = (a), i##END = (b); i >= i##END; i--)
#define N ((int)(5e6 + 15))

char ck[N]; int inv, pcnt, mod, p[N], phi[N], sum_phi[N];
unordered_map<int, int> vis;
void add(int &x, int y) { (x += y) >= mod ? x -= mod : 0; }
int qpow(int a, int b)
{
    int r = 1;
    for(a %= mod; b > 0; b >>= 1, a = a * a % mod)
        if(b & 1) r = r * a % mod;
    return r;
}
void init(const int _ = N - 5)
{
    phi[1] = 1;
    rep(i, 2, _)
    {
        if(!ck[i]) { p[++pcnt] = i; phi[i] = i - 1; }
        for(int j = 1; j <= pcnt && i * p[j] <= _; j++)
        {
            const int pos = i * p[j]; ck[pos] = true;
            if(i % p[j]) { phi[pos] = phi[i] * phi[p[j]]; }
            else { phi[pos] = phi[i] * p[j]; break; }
        }
    }
    rep(i, 1, _) add(sum_phi[i], (sum_phi[i - 1] + i * i % mod * phi[i] % mod) % mod);
}
int fun2(int x) { x %= mod; return x * (x + 1) % mod * (2 * x % mod + 1) % mod * inv % mod; }
int fun3(int x) { x %= mod; const int y = (x * (x + 1) / 2) % mod; return y * y % mod; }
int getsum(int n)
{
    if(n <= N - 5) return sum_phi[n];
    if(vis[n]) return vis[n];
    int res = fun3(n);
    for(int l = 2, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        add(res, mod - (fun2(r) - fun2(l - 1) + mod) % mod * getsum(n / l) % mod);
    }
    return vis[n] = res;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n; cin >> mod >> n; init();
    int res = 0; inv = qpow(6ll, mod - 2);
    for(int l = 1, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        add(res, (getsum(r) - getsum(l - 1) + mod) % mod * fun3(n / l) % mod);
    }
    cout << res << '\n';
    return 0;
}