题面:bzoj3930
题意:从 [ L , H ] [L,H] [L,H]选出 n n n个数,使它们的最大公约数为 k k k的方案数,对 1 e 9 + 7 1e9+7 1e9+7取模。
1 ≤ n , k , L , H ≤ 1 0 9 1 \le n, k, L, H \le 10^9 1≤n,k,L,H≤109, H − L ≤ 1 0 5 H-L \le 10^5 H−L≤105
题解:即求 ∑ a 1 = L H ∑ a 2 = L H … ∑ a n = L H [ g c d i = 1 n ( a i ) = k ] \sum _{a_1=L} ^H \sum_{a_2=L} ^H … \sum_{a_n=L} ^H [gcd _{i=1} ^n (a_i) = k] ∑a1=LH∑a2=LH…∑an=LH[gcdi=1n(ai)=k]
令 f ( x ) = ∑ a 1 = L H ∑ a 2 = L H … ∑ a n = L H [ g c d i = 1 n ( a i ) = x ] f(x)=\sum _{a_1=L} ^H \sum_{a_2=L} ^H … \sum_{a_n=L} ^H [gcd _{i=1} ^n (a_i) = x] f(x)=∑a1=LH∑a2=LH…∑an=LH[gcdi=1n(ai)=x]
g ( x ) = ∑ x ∣ d f ( d ) = ( ⌊ H x ⌋ − ⌊ L − 1 x ⌋ ) n g(x) = \sum _{x \mid d} f(d) = (\lfloor \frac H x \rfloor - \lfloor \frac {L-1} x \rfloor) ^n g(x)=∑x∣df(d)=(⌊xH⌋−⌊xL−1⌋)n
莫比乌斯反演: f ( x ) = ∑ x ∣ d μ ( d x ) g ( d ) = ∑ x ∣ d μ ( d x ) ( ⌊ H d ⌋ − ⌊ L − 1 d ⌋ ) n f(x) = \sum _{x \mid d} \mu (\frac d x) g(d) = \sum _{x \mid d} \mu(\frac d x) (\lfloor \frac H d \rfloor - \lfloor \frac {L-1} d \rfloor)^n f(x)=∑x∣dμ(xd)g(d)=∑x∣dμ(xd)(⌊dH⌋−⌊dL−1⌋)n
答案即为 f ( k ) = ∑ k ∣ d μ ( d k ) ( ⌊ H d ⌋ − ⌊ L − 1 d ⌋ ) n = ∑ i = 1 ⌊ H k ⌋ μ ( i ) ( ⌊ H k i ⌋ − ⌊ L − 1 k i ⌋ ) n f(k) = \sum _{k \mid d} \mu (\frac d k) (\lfloor \frac H d \rfloor - \lfloor \frac {L-1} d \rfloor)^n = \sum _{i=1} ^{\lfloor \frac H k \rfloor} \mu (i) (\lfloor \frac H {ki} \rfloor - \lfloor \frac {L-1} {ki} \rfloor)^n f(k)=∑k∣dμ(kd)(⌊dH⌋−⌊dL−1⌋)n=∑i=1⌊kH⌋μ(i)(⌊kiH⌋−⌊kiL−1⌋)n
明显可以数论分块。但是 ⌊ H k ⌋ \lfloor \frac H k \rfloor ⌊kH⌋可能达到 1 e 9 1e9 1e9,需要快速对 ∑ i = 1 n μ ( i ) \sum _{i=1} ^n \mu (i) ∑i=1nμ(i)求和。
考虑杜教筛:
令 S ( n ) = ∑ i = 1 n μ ( i ) S(n) = \sum _{i=1} ^n \mu (i) S(n)=∑i=1nμ(i)
g ( 1 ) S ( n ) = ∑ i = 1 n ( g ∗ μ ) ( i ) − ∑ i = 2 n g ( i ) S ( n i ) g(1)S(n) = \sum _{i=1} ^n (g * \mu) (i) - \sum _{i=2} ^n g(i)S(\frac n i) g(1)S(n)=∑i=1n(g∗μ)(i)−∑i=2ng(i)S(in)
令 g ( i ) = 1 g(i) = 1 g(i)=1
S ( n ) = 1 − ∑ i = 2 n S ( n i ) S(n)=1 - \sum _{i=2} ^n S(\frac n i) S(n)=1−∑i=2nS(in)
记忆化算出 S ( n ) S(n) S(n)
#include <bits/stdc++.h>
using namespace std;
typedef long long ll ;
const int maxn = 1e7 + 10, mod = 1e9 + 7 ;
const ll INF = 1e18 ;
int prime[maxn], tot ;
ll mu[maxn] ;
bool vis[maxn] ;
map<ll, ll> val ;
void init () {
mu[1] = 1 ;
for (int i = 2; i < maxn; i ++) {
if (!vis[i]) prime[++ tot] = i, mu[i] = -1 ;
for (int j = 1; j <= tot && i * prime[j] < maxn; j ++) {
vis[i * prime[j]] = 1 ;
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0; break ;
}
mu[i * prime[j]] = -mu[i] ;
}
}
for (int i = 1; i < maxn; i ++) mu[i] = (mu[i - 1] + mu[i]) % mod ;
}
inline ll power (ll x, ll y) {
ll res = 1 ;
while (y) {
if (y & 1) res = res * x % mod ;
x = x * x % mod; y >>= 1 ;
}
return res ;
}
ll cal (ll n) {
if (n < maxn) return mu[n] ;
if (val[n]) return val[n] ;
ll res = 1 ;
for (ll i = 2, nxt = 0; i <= n; i = nxt + 1) {
nxt = n / (n / i) ;
res = (res - (nxt - i + 1) * cal (n / i) % mod) % mod ;
}
return val[n] = (res + mod) % mod ;
}
int main() {
ll n, k, l, h ;
cin >> n >> k >> l >> h ;
init () ;
-- l /= k; h /= k ;
ll ans = 0 ;
for (ll i = 1, nxt = 0; i <= h; i = nxt + 1) {
nxt = min (l / i ? l / (l / i) : INF, h / (h / i)) ;
ans = (ans + power (h / i - l / i, n) * (cal (nxt) - cal (i - 1)) % mod) % mod ;
}
cout << (ans + mod) % mod << endl ;
return 0 ;
}