简单的数学题

求

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

$n \leq 10^{10}$

---

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

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

$$= \sum_{d=1}^n d^3 \sum_{i=1}^{n/d} \sum_{j=1}^{n/d} ij \sum_{e|gcd(i, j)} \mu(e)$$

$$= \sum_{d=1}^n d^3 \sum_{i=1}^{n/d} \sum_{j=1}^{n/d} ij \sum_{e|i 且 e|j} \mu(e)$$

$$= \sum_{d=1}^n d^3 \sum_{e=1}^{n/d} \mu(e) e^2 \sum_{i=1}^{n/de} \sum_{j=1}^{n/de} ij$$

设$f(x) = \sum_{i=1}^x i$

$$ans = \sum_{d=1}^n d^3 \sum_{e=1}^{n/d} \mu(e)e^2 f^2(\dfrac{n}{de})$$

我们知道$de \leq n$
设$T = de$

$$ans = \sum_{T=1}^n \sum_{d|T} d^3\mu(\dfrac{T}{d}) (\dfrac{T}{d})^2 f^2(\dfrac{n}{T})$$

$$= \sum_{T=1}^n f^2(\dfrac{n}{T}) \sum_{d|T} d^3\mu(\dfrac{T}{d}) (\dfrac{T}{d})^2$$

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

我们知道

$$\sum_{d|n} \dfrac{\mu(d)}{d} = \dfrac{\phi(n)}{n}$$

$$\sum_{d|n} \mu(d) \dfrac{n}{d} = \phi(n)$$

$$\therefore ans = \sum_{T=1}^n T^2 f^2(\dfrac{n}{T}) \sum_{d|T} \dfrac{T}{d} \ \mu(d)$$

$$= \sum_{T=1}^n T^2 f^2(\dfrac{n}{T}) \phi(T)$$

我们知道

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

$$f(x) = \sum_{i=1}^x = \dfrac{x(x+1)}{2}$$

所以我们只需杜教筛出

$$S(n) = \sum_{i=1}^n i^2 \phi(i)$$
即可在 $O(n^{\frac{2}{3}})$内解决此题

code：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5000000;
const int LIM = 4000000;

ll MOD, INV_6, INV_2;

map <ll, ll> mp; 

bool is_prime[N]; 
int tot = 0, prime[N]; 
ll n, s[N], euler[N];

inline ll ksm( ll a, ll b )
{
    ll ret = 1; a %= MOD;
    for( ; b; b >>= 1, a = a * a %MOD )
        if( b & 1 ) ret = ret * a %MOD;
    return ret;
}

inline ll calc_sum( ll x ) {
    x %= MOD; return x * ( x + 1 ) %MOD * INV_2 %MOD;
}

inline ll calc_poi( ll x ) {
    x %= MOD; return x * ( x + 1 ) %MOD * ( 2*x + 1 ) %MOD * INV_6 %MOD;
}

inline void init( int size )
{   
    memset( is_prime, true, sizeof( is_prime ) );
    is_prime[1] = false; euler[1] = 1;

    for( int i = 2; i <= size; i ++ )
    {
        if( is_prime[i] )
            prime[++tot] = i, euler[i] = i-1;

        for( int j = 1; j <= tot; j ++ )
        {
            int p = prime[j];
            if( i * p > size ) break;

            is_prime[i*p] = false;
            if( i % p == 0 )
            {
                euler[i*p] = euler[i] * p;
                break;
            }

            euler[i*p] = euler[i] * euler[p];
        }
    }

    for( int i = 1; i <= size; i ++ )
        s[i] = ( s[i-1] + (ll)euler[i] * i %MOD * i %MOD ) %MOD;
}

inline ll calc_euler( ll x )
{
    if( x <= LIM ) return s[x];
    if( mp[x] ) return mp[x];

    ll next = 0, ret = calc_sum( x ) * calc_sum( x ) %MOD;
    for( ll i = 2; i <= x; i = next + 1 )
    {
        next = x / ( x / i );
        ll tmp = ( calc_poi( next ) - calc_poi( i - 1 ) ) %MOD;
        if( tmp < 0 ) tmp += MOD;

        ret -= calc_euler( x / i ) * tmp %MOD;
        while( ret < 0 ) ret += MOD;
    }

    mp[x] = ret; 
    return ret;
}

int main()
{ 
    scanf( "%lld%lld", &MOD, &n );

    INV_2 = ksm( 2, MOD-2 );
    INV_6 = ksm( 6, MOD-2 );
    init( LIM );

    ll next = 0, ans = 0;
    for( ll i = 1; i <= n; i = next + 1 )
    {
        next = n / ( n / i );
        ll tmp = calc_euler( next ) - calc_euler( i - 1 );
        while( tmp < 0 ) tmp += MOD;
        ll _tmp = calc_sum( n / i ) * calc_sum( n / i ) %MOD; 

        ans += _tmp * tmp %MOD;
        while( ans >= MOD ) ans -= MOD;
    }

    printf( "%lld\n", ans );

    return 0;
}