洛谷 3768 简单的数学题

zzq 的题啊……

直接开推 n∑i=1n∑j=1ijgcd(ij)=n∑d=1d∑d|i∑d|j[gcd(i,j)=d]ij=n∑d=1d3∑id≤n∑jd≤m[gcd(i,j)=1]ij=n∑d=1d3∑id≤ni∑jd≤nj∑k|gcd(i,j)μ(k)=n∑d=1d3⌊nd⌋∑k=1μ(k)k2(∑idk≤ni)2

令 T=dk，更改枚举方式，有 n∑d=1d3⌊nd⌋∑k=1μ(k)k2(∑idk≤ni)2=n∑T=1(∑iT≤ni)2∑d|Td3μ(Td)T2d2=n∑T=1(∑iT≤ni)2T2∑d|Tdμ(Td)=n∑T=1(∑iT≤ni)2T2φ(T)

杜教筛时，考虑设 f(i)=i2φ(i),g(i)=i2,h(i)=(f∗g)(i)=∑d|id2φ(d)i2d2=d3。
同时要知道 n∑i=1i2=n(n+1)(2n+1)6,n∑i=1i3=n2(n+1)24。

代码：

#include <cstdio>
#include <tr1/unordered_map>
using namespace std;
using namespace tr1;
const long long N = 1e10;
const int CNT = 1e7;
long long n,mod,inv2,inv6;
int vis[CNT + 5],prime[CNT + 5],cnt;
long long phi[CNT + 5];
unordered_map<long long,long long> w;
inline long long fpow(long long a,long long b)
{
    long long ret = 1;
    for(;b;b >>= 1)
        (b & 1) && (ret = ret * a % mod),a = a * a % mod;
    return ret;
}
inline long long sum(long long n)
{
    n %= mod;
    return n * (n + 1) % mod * inv2 % mod;
}
inline long long sum1(long long n)
{
    n %= mod;
    return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
}
long long calc(long long n)
{
    if(n <= CNT)
        return phi[n];
    if(w.count(n))
        return w[n];
    long long ret = sum(n) * sum(n) % mod;
    for(long long l = 2,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ret = (ret - (sum1(r) - sum1(l - 1)) * calc(n / l) % mod) % mod;
    }
    return w[n] = (ret + mod) % mod;
}
long long ans;
int main()
{
    scanf("%lld%lld",&mod,&n),inv2 = fpow(2,mod - 2),inv6 = fpow(6,mod - 2);
    phi[1] = 1;
    for(register int i = 2;i <= CNT;++i)
    {
        if(!vis[i])
            phi[prime[++cnt] = i] = i - 1;
        for(register int j = 1;j <= cnt && i * prime[j] <= CNT;++j)
        {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j]))
            {
                phi[i * prime[j]] = phi[i] * prime[j] % mod;
                break;
            }
            else
                phi[i * prime[j]] = phi[i] * (prime[j] - 1) % mod;
        }
    }
    for(register int i = 1;i <= CNT;++i)
        phi[i] = (phi[i - 1] + phi[i] * i % mod * i % mod) % mod;
    for(register long long l = 1,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ans = (ans + sum(n / l) * sum(n / l) % mod * (calc(r) - calc(l - 1)) % mod) % mod;
    }
    printf("%lld\n",(ans + mod) % mod);
}