<题解>简单的数学题

洛谷的题目链接

完全没有任何掩饰，求\Big(\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\Big)\bmod p

化简柿子：

\begin{aligned} &\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\\ =&\sum_{i=1}^n\sum_{j=1}^nij\sum_{d|\gcd(i,j)}\varphi(d)\\ =&\sum_{i=1}^n\sum_{j=1}^nij\sum_{d=1}^n[d|i][d|j]\varphi(d)\\ =&\sum_{d=1}^n\varphi(d)\sum_{i=1}^{\lfloor n/d\rfloor}\sum_{j=1}^{\lfloor n/d\rfloor}ijd^2\\ =&\sum_{d=1}^n\varphi(d)d^2(\sum_{i=1}^{\lfloor n/d\rfloor}i)^2\\ \end{aligned}

令h(n) = (\sum_{i=1}^ni)^2，常数时间算出，那么原式=\sum_{d=1}^n\varphi(d)d^2h(\lfloor\dfrac n d\rfloor)，前面杜教筛，后面整除分块。

那么我们就要用杜教筛算\sum_{d=1}^n\varphi(d)d^2：
f=\varphi\cdot\operatorname{ID}^2
g=?
(f\times g)(n)=\sum_{d|n}\varphi(d)d^2g(\dfrac n d)

令g=\operatorname{ID}^2，(f\times g)(n)=\sum_{d|n}\varphi(d)d^2\dfrac {n^2} {d^2}=\sum_{d|n}\varphi(d)n^2 = n^3（Tips:\varphi\times1=\operatorname{ID}）

1^2+2^2+3^2+\dots+n^2=\dfrac {n(n+1)(2n+1)} 6
1^3+2^3+3^3+\dots+n^3=\dfrac {n^2(n+1)^2} 4

然后就是背公式了：g(1)S(n)=\sum_{i=1}^n(f\times g)(i)-\sum_{i=2}^ng(i)S(\lfloor \dfrac n i\rfloor)

#include <cstdio>

const int PRE = 5000000, OUT = 3000;

long long n, p, small[PRE], big[OUT], inv6, inv2;
bool np[PRE];
int ps[PRE];
long long phi[PRE];

inline void pre(void)
{//线筛
    phi[1] = np[1] = small[1] = true;
    for (int i = 2; i != PRE; ++i)
    {
        if (np[i] == false)
        {
            ps[++ps[0]] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= ps[0] && i * ps[j] < PRE; ++j)
        {
            np[i * ps[j]] = true;
            if (i % ps[j] == 0)
            {
                phi[i * ps[j]] = phi[i] * ps[j] % p;
                break;
            }
            phi[i * ps[j]] = phi[i] * (ps[j] - 1) % p;
        }
        small[i] = phi[i] * i % p * i % p;
        small[i] = (small[i] + small[i - 1]) % p;
    }
}

inline long long sum2(long long v)
{//g
    v %= p;
    return inv6 * v % p * (v + 1) % p * (v + v + 1) % p;
}

inline long long sum3(long long v)
{//f*g
    v %= p;
    long long ret = v * (v + 1) % p * inv2 % p;
    return ret * ret % p;
}

inline long long sum1(long long v)
{//h
    v %= p;
    long long ret = v * (v + 1) % p * inv2 % p;
    return ret * ret % p;
}

long long S(long long c)
{
    if (c < PRE)
        return small[c];
    int x = n / c;//不用map的做法
    if (big[x])
        return big[x];
    long long ret = sum3(c);
    for (long long i = 2, j; i <= c; i = j + 1)
    {
        j = c / (c / i);
        ret -= (sum2(j) - sum2(i - 1) + p) % p * S(c / i);
        ret = (ret + p) % p;
    }
    return big[x] = ret;
}

inline int power(long long n, int m)
{//用于算逆元
    long long ret = 1;
    while (m)
    {
        if (m & 1)
        {
            ret *= n;
            ret %= p;
        }
        n *= n;
        n %= p;
        m >>= 1;
    }
    return ret;
}

int main(void)
{
    scanf("%lld %lld", &p, &n);
    pre();
    inv6 = power(6, p - 2);
    inv2 = power(2, p - 2);
    long long ans = 0;
    for (long long i = 1, j; i <= n; i = j + 1)
    {
        j = n / (n / i);
        long long S1 = S(j);
        long long S2 = S(i - 1);
        ans += (S1 - S2 + p) % p * sum1(n / i);
        ans %= p;
    }
    printf("%lld", ans);
}