简单的数学题

求

∑i=1n∑j=1nij gcd(i,j) ∑ i = 1 n ∑ j = 1 n i j   g c d ( i , j )

n≤1010 n ≤ 10 10

---

ans=∑d=1nd∑i=1n∑j=1nij [gcd(i,j)=d] a n s = ∑ d = 1 n d ∑ i = 1 n ∑ j = 1 n i j   [ g c d ( i , j ) = d ]

=∑d=1nd3∑i=1n/d∑j=1n/dij [gcd(i,j)=1] = ∑ d = 1 n d 3 ∑ i = 1 n / d ∑ j = 1 n / d i j   [ g c d ( i , j ) = 1 ]

=∑d=1nd3∑i=1n/d∑j=1n/dij∑e|gcd(i,j)μ(e) = ∑ d = 1 n d 3 ∑ i = 1 n / d ∑ j = 1 n / d i j ∑ e | g c d ( i , j ) μ ( e )

=∑d=1nd3∑i=1n/d∑j=1n/dij∑e|i且e|jμ(e) = ∑ d = 1 n d 3 ∑ i = 1 n / d ∑ j = 1 n / d i j ∑ e | i 且 e | j μ ( e )

=∑d=1nd3∑e=1n/dμ(e)e2∑i=1n/de∑j=1n/deij = ∑ d = 1 n d 3 ∑ e = 1 n / d μ ( e ) e 2 ∑ i = 1 n / d e ∑ j = 1 n / d e i j

设 f(x)=∑xi=1i f ( x ) = ∑ i = 1 x i

ans=∑d=1nd3∑e=1n/dμ(e)e2f2(nde) a n s = ∑ d = 1 n d 3 ∑ e = 1 n / d μ ( e ) e 2 f 2 ( n d e )

我们知道 de≤n d e ≤ n 设 T=de T = d e

ans=∑T=1n∑d|Td3μ(Td)(Td)2f2(nT) a n s = ∑ T = 1 n ∑ d | T d 3 μ ( T d ) ( T d ) 2 f 2 ( n T )

=∑T=1nf2(nT)∑d|Td3μ(Td)(Td)2 = ∑ T = 1 n f 2 ( n T ) ∑ d | T d 3 μ ( T d ) ( T d ) 2

=∑T=1nT2f2(nT)∑d|Td μ(Td) = ∑ T = 1 n T 2 f 2 ( n T ) ∑ d | T d   μ ( T d )

我们知道

∑d|nμ(d)d=ϕ(n)n ∑ d | n μ ( d ) d = ϕ ( n ) n

∑d|nμ(d)nd=ϕ(n) ∑ d | n μ ( d ) n d = ϕ ( n )

∴ans=∑T=1nT2f2(nT)∑d|TTd μ(d) ∴ a n s = ∑ T = 1 n T 2 f 2 ( n T ) ∑ d | T T d   μ ( d )

=∑T=1nT2f2(nT)ϕ(T) = ∑ T = 1 n T 2 f 2 ( n T ) ϕ ( T )

我们知道

∑i=1ni2=n(n+1)(2n+1)6 ∑ i = 1 n i 2 = n ( n + 1 ) ( 2 n + 1 ) 6

f(x)=∑i=1x=x(x+1)2 f ( x ) = ∑ i = 1 x = x ( x + 1 ) 2

所以我们只需杜教筛出

S(n)=∑i=1ni2ϕ(i) S ( n ) = ∑ i = 1 n i 2 ϕ ( i ) 即可在 O(n23) O ( n 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; }