【题解】简单的数学题 杜教筛 ——quhengyi11

喵链

题目喵述

求
$$(\sum_{i=1}^n \sum_{j=1}^nijgcd(i,j)) \% p$$

题解

杜教筛模板题

原式=

$$\sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)==d]ij$$

显然有

$$\sum_{d=1}^n d^3\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac n d}[gcd(i,j)==1]ij$$

令

$$f(d)=\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)==d]$$

$$sum(x)=\frac {(1+x) \times x} 2$$

$$
\begin{align}
g(d)&=\sum_{d|n}f(n)\\
&=\sum_{i=1}^n\sum_{j=1}^nij[d|gcd(i,j)]\\
&=\sum_{i=1}^{\frac n d}\sum_{j=1}^{\frac n d}d^2ij\\
&=d^2sum^2(\frac n d)
\end{align}
$$

所以

$$f(1)=\sum_{d=1}^n \mu(d)g(d)$$

把 $g$展开并代回原式

$$ans = \sum_{d=1}^n d^3\sum_{p=1}^{\frac n d}\mu(p)p^2sum^2(\frac n {pd})$$

我们可以交换枚举顺序，令 $T=pd$

$$\sum_{T=1}^n sum^2(\frac n {pd})\sum_{d|T} d^3 (\frac T d)^2 \mu(\frac T d)$$

化简顺便把 $T$提出去

$$\sum_{T=1}^n sum^2(\frac n T)T^2\sum_{d|T} d \mu(\frac T d)$$

由 $(Id * \mu)(n) = \phi(n)$

$$\sum_{T=1}^n sum^2(\frac n T)T^2\phi(T)$$

观察这个式子，前面的部分可以数论分块 $\sqrt n$求，后面那一部分就是杜教筛基本套路了

令

$$f(n)=n^2\phi(n)$$

$$s(n)=\sum_{i=1}^n f(i)$$

因为 $(1 * \phi)(n) = Id(n)$

令

$$g(n)=n^2$$

则

$$
\begin{align}
g(n)&=(f*g)(n) – \sum_{d=2}^n f(d)g(\frac n d)\\
&= \sum_{i=1}^n i^3 – g(d)s(\frac n d)
\end{align}
$$

然后用数论分块+递归求这个式子即可，

还有一些辅助式子

$$\sum_{i=1}^n i^3 = \frac {i^2 (i+1)^2}{4}$$

$$\sum_{i=1}^n i^2 = \frac {(2n+1)(n+1)n} 6$$

讲真这个自然数幂求和有没有什么有效的推公式法呀我基本都是记的（平方和有一个著名的三角构造，立方和就不知道了 QAQ）

注意取模，试一下极限数据即可

#include<bits/stdc++.h>
#define fo(i, a, b) for (int i = (a); i <= (b); ++i)

#define N 10000005
#define ll long long
std::bitset<N> vis;
ll phi[N], p[N], tot, mod, up;
std::map<ll, ll> mp;
inline void init ()
{
    up = N - 3;
    phi[1] = 1;
    fo (i, 2, up)
    {
        if (!vis[i])
            p[++tot] = i, phi[i] = i - 1;
        vis[i] = 1;
        for (int j = 1; j <= tot && i * p[j] <= up; ++j)
        {
            vis[i * p[j]] = 1;
            if (!(i % p[j])) 
            {
                phi[i * p[j]] = phi[i] * p[j];
                break;
            }
            phi[i * p[j]] = phi[i] * phi[p[j]];
        }
    }
    fo (i, 2, up) phi[i] = (phi[i] * i % mod * i + phi[i - 1]) % mod;
}
ll n, m, inv6;
inline ll pow (ll x, ll y = mod - 2)
{
    ll ret = 1;
    while (y)
    {
        if (y & 1) ret = ret * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ret;
}
inline ll sqr (ll x) {return x * x % mod;}
inline ll sum (ll x) {x %= mod; return ((x + 1) * x >> 1) % mod;}
inline ll md (ll x) {x %= mod; return (x < 0) ? x + mod : x;}
inline ll cube (ll x) {x %= mod; return (x << 1 | 1) % mod * (x + 1) % mod * x % mod * inv6 % mod;}
inline ll s (ll x)
{
    if (x <= up) return phi[x];
    if (mp[x]) return mp[x];
    ll ret = sqr(sum(x));
    for (ll i = 2, j; i <= x; i = j + 1)
    {
        j = x / (x / i);
        ret = (ret - (cube(j) - cube(i - 1)) % mod * s(x / i)) % mod;
    }
    return mp[x] = md(ret);
}
int main()
{
    scanf("%lld %lld", &mod, &n);
    init(); inv6 = pow(6);
    ll ans = 0;
    for (ll i = 1, j; i <= n; i = j + 1)
    {
        j = n / (n / i);
        ans = (ans + sqr(sum(n / i)) * (s(j) - s(i - 1))) % mod;
    }
    printf("%lld\n", md(ans));
    return 0;
}