问题分析
题意还是很好理解的，就是求
[Ans=f(k)=sum_{a_1=L}^Hsum_{a_2=L}^Hcdotssum_{a_N=L}^H[gcd(a_1,a_2,cdots,a_N)=K]
]
我们令
[F(k)=sum_{a_1=L}^Hsum_{a_2=L}^Hcdotssum_{a_N=L}^H[K|gcd(a_1,a_2,cdots,a_N)]=Num^N(k)\
(Num(k)=lfloorfrac{H}{k}
floor-max(0,lfloorfrac{L-1}{k}
floor))
]
就有
[F(k)=sum_{k|d}f(d)
]
所以
[f(k)=sum_{k|d}mu(frac{d}{k})F(d)
]
[Ans=f(k)=sum_{i=1}^{lfloorfrac{H}{k}
floor}mu(i)Num^N(ik)
]
(Num)可以整除分块，而(mu)由于数据范围大需要使用杜教筛。
考虑
[f(n)=mu(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_0)，那么
[S(n)=1-sum_{i=2}^nS(lfloorfrac{n}{i}
floor)
]
其中(*)代表狄利克雷卷积，(id_0(k)=k^0=1)。
这样就好了，时间复杂度(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 INF = 10000000000000000;
const LL Threshold = 2000000;
const LL Mod = 1000000007;
LL Mu[ Threshold + 10 ], Vis[ Threshold + 10 ], Prime[ Threshold + 10 ], Num, Sum[ Threshold + 10 ];
LL N, K, L, H, Ans;
map< LL, LL > Map;
void EulerSieve();
void Cal();
int main() {
EulerSieve();
scanf( "%lld%lld%lld%lld", &N, &K, &L, &H );
Cal();
printf( "%lld
", Ans );
return 0;
}
void EulerSieve() {
Mu[ 1 ] = 1;
for( LL i = 2; i <= Threshold; ++i ) {
if( !Vis[ i ] ) {
Mu[ i ] = -1;
Prime[ ++Num ] = i;
}
for( LL j = 1; j <= Num && i * Prime[ j ] <= Threshold; ++j ) {
Vis[ i * Prime[ j ] ] = 1;
if( i % Prime[ j ] == 0 ) break;
Mu[ i * Prime[ j ] ] = -Mu[ i ];
}
}
for( LL i = 1; i <= Threshold; ++i ) Sum[ i ] = ( Sum[ i - 1 ] + Mu[ i ] ) % Mod;
return;
}
LL Power( LL k, LL t ) {
if( t == 0 ) return 1;
LL T = Power( k, t >> 1 );
T = T * T % Mod;
if( t & 1 ) T = k % Mod * T % Mod;
return T;
}
LL S( LL t ) {
if( t <= Threshold ) return Sum[ t ];
if( Map.find( t ) != Map.end() ) return Map[ t ];
LL Ans = 1;
for( LL i = 2, j; i <= t; i = j + 1 ) {
j = t / ( t / i );
Ans = ( ( Ans - ( j - i + 1 ) % Mod * S( t / i ) % Mod ) % Mod + Mod ) % Mod;
}
Map[ t ] = Ans;
return Ans;
}
void Cal() {
LL l = ( L - 1 ) / K, r = H / K;
for( LL i = 1, j; i <= r; i = j + 1 ) {
j = min( ( l / i ) ? l / ( l / i ) : INF, r / ( r / i ) );
Ans = ( Ans + ( ( S( j ) - S( i - 1 ) ) % Mod + Mod ) % Mod * Power( r / i - l / i, N ) % Mod ) % Mod;
}
return;
}