问题分析
题目就是要求
[Ans=sum_{i=1}^nsum_{j=1}^nijgcd(i,j)
]
首先乱搞一波：
[egin{aligned}
Ans&=sum_{t=1}^nsum_{i=1}^nsum_{j=1}^nijt[gcd(i,j)==t]\
&=sum_{t=1}^nt^3sum_{i=1}^{lfloorfrac{n}{t}
floor}sum_{j=1}^{lfloorfrac{n}{t}
floor}ij[gcd(i,j)==1]
end{aligned}
]
然后对后面那一部做莫比乌斯反演：
[f(d)=sum_{i=1}^nsum_{j=1}^nij[gcd(i,j)==1]\
F(d)=sum_{i=1}^nsum_{j=1}^nij[d|gcd(i,j)]=d^2sum_{i=1}^{lfloorfrac{n}{d}
floor}sum_{j=1}^{lfloorfrac{n}{d}
floor}ij\
F(d)=sum_{d|k}f(k)\
f(d)=sum_{d|k}mu(frac{k}{d})F(k)=sum_{d|k}mu(frac{k}{d})k^2sum_{i=1}^{lfloorfrac{n}{k}
floor}sum_{j=1}^{lfloorfrac{n}{k}
floor}ij\
f(1)=sum_{k=1}^nmu(k)k^2sum_{i=1}^{lfloorfrac{n}{k}
floor}sum_{j=1}^{lfloorfrac{n}{k}
floor}ij=sum_{k=1}^nmu(k)k^2(sum_{i=1}^{lfloorfrac{n}{k}
floor}i)^2
]
那么就有：
[Ans=sum_{t=1}^{n}t^3sum_{k=1}^{n/t}mu(k)k^2(sum_{i=1}^{lfloorfrac{n}{tk}
floor}i)^2
]
令(K=tk)，
[egin{aligned}
Ans&=sum_{K=1}^nsum_{t|K}tmu(frac{K}{t})K^2(sum_{i=1}^{lfloorfrac{n}{K}
floor}i)^2\
&=sum_{K=1}^nK^2(sum_{i=1}^{lfloorfrac{n}{K}
floor}i)^2sum_{t|K}tmu(frac{K}{t})
end{aligned}
]
因为(phi*1=id_1)，(1*mu=e)，所以(id_1*mu=phi*1*mu=phi*e=phi)。（其中(*)表示狄利克雷卷积）
那么
[sum_{t|K}frac{K}{t}mu(t)=phi(K)\
Ans=sum_{K=1}^nK^2(sum_{i=1}^{lfloorfrac{n}{K}
floor}i)^2sum_{t|K}frac{K}{t}mu(t)=sum_{K=1}^n(sum_{i=1}^{lfloorfrac{n}{K}
floor}i)^2K^2phi(K)
]
其中((sum_{i=1}^{lfloorfrac{n}{K} floor}i)^2)可以整除分块，后面我们考虑杜教筛：
[f(n)=n^2phi(n)\
S(n)=sum_{i=1}^nf(i)\
g(1)S(n)=sum_{i=1}^n(g*f)(i)-sum_{i=2}^ng(i)S(lfloorfrac{n}{i}
floor)\
g=id_2\
S(n)=sum_{i=1}^ni^3-sum_{i=2}^ni^2S(lfloorfrac{n}{i}
floor)
]
那么就可以啦！时间复杂度(O(n^{frac{2}{3}}))。
参考程序
这个写法不优良，时间复杂度是 (O(n^{frac{3}{4}})) 的。还要带上一个 map 的大常数 log ?
#include <bits/stdc++.h>
#define LL long long
using namespace std;
const LL Maxn066 = 6000000;
LL Phi[ Maxn066 + 10 ], Vis[ Maxn066 + 10 ], Prime[ Maxn066 ], Count;
LL Mod, n, Ans, Inv6;
void EulerSieve();
void Cal();
LL Inv( LL a );
int main() {
scanf( "%lld%lld", &Mod, &n );
Inv6 = Inv( 6 );
EulerSieve();
Cal();
return 0;
}
void EulerSieve() {
Phi[ 1 ] = 1;
for( LL i = 2; i <= Maxn066; ++i ) {
if( !Vis[ i ] ) {
Phi[ i ] = i - 1;
Prime[ ++Count ] = i;
}
for( LL j = 1; j <= Count && Prime[ j ] * i <= Maxn066; ++j ) {
Vis[ Prime[ j ] * i ] = 1;
if( !( i % Prime[ j ] ) ) {
Phi[ Prime[ j ] * i ] = Phi[ i ] * Prime[ j ] % Mod;
break;
}
Phi[ Prime[ j ] * i ] = Phi[ i ] * Phi[ Prime[ j ] ] % Mod;
}
}
for( LL i = 2; i <= Maxn066; ++i ) Phi[ i ] = ( Phi[ i - 1 ] + Phi[ i ] * i % Mod * i % Mod ) % Mod;
return;
}
map< LL, LL > Map;
LL S( LL Key ) {
if( Key <= Maxn066 ) return Phi[ Key ];
if( Map.find( Key ) != Map.end() ) return Map[ Key ];
LL Ans = Key % Mod * ( Key % Mod + 1 ) / 2 % Mod; Ans = Ans * Ans % Mod;
for( LL i = 2, j; i <= Key; i = j + 1 ) {
j = Key / ( Key / i );
LL T = j % Mod * ( j % Mod + 1 ) % Mod * ( 2 * j % Mod + 1 ) % Mod * Inv6 % Mod;
T = T - ( i % Mod - 1 ) * ( i % Mod ) % Mod * ( 2 * i % Mod - 1 ) % Mod * Inv6 % Mod;
T = ( T % Mod + Mod ) % Mod;
T = T * S( Key / i ) % Mod;
Ans = ( ( Ans - T ) % Mod + Mod ) % Mod;
}
Map[ Key ] = Ans;
return Ans;
}
void Cal() {
for( LL i = 1, j; i <= n; i = j + 1 ) {
j = n / ( n / i );
LL T = ( 1 + n / i % Mod ) * ( n / i % Mod ) / 2 % Mod; T = T * T % Mod;
LL Delta = ( ( S( j ) - S( i - 1 ) ) % Mod + Mod ) % Mod;
Ans = ( Ans + T * Delta % Mod ) % Mod;
}
printf( "%lld
", Ans );
return;
}
void Exgcd( LL a, LL b, LL &x, LL &y ) {
if( b == 0 ) {
x = 1; y = 0; return;
}
Exgcd( b, a % b, y, x );
y -= a / b * x;
return;
}
LL Inv( LL a ) {
LL x, y;
Exgcd( a, Mod, x, y );
if( x < 0 ) x += Mod;
return x;
}