bzoj3930 CQOI2015 選數

題面：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∣d​f(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=2n​g(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=2n​S(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 ;
}