這題做的時候接連想錯了好多次……但是回到正軌上之後依然是一個套路題。（不過這題好像有比莫比烏斯反演更好的做法，莫比烏斯反演貌似是某種能過的暴力ヽ（´ー｀）┌）不過能過也就行了吧哈哈。

　　首先我們把數字的範圍要進行縮小：最大公約數為 \(K\) 那自然所有選出來的數都必須是 \(K\) 的倍數。所以我們改選數為選擇是 \(K\) 的多少倍。然後由於是最大公約數，所以選出來的這些數必須最大公約數等於\(1\)。實際上多個數的最大公約數\( = 1\)完全可以和兩個數的最大公約數 \( = 1\) 用一樣的方法去反演。只不過這題由於數據範圍非常的大，所以處理 \(\mu\) 的前綴和必須要使用杜教篩。

#include <bits/stdc++.h>
using namespace std;
#define maxn 1000300
#define db double
#define int long long
int maxx = maxn - 1e2, mod = 1e9 + 7;
int N, K, L, H, ans, Sum[maxn];
int tot, pri[maxn];
map <int, int> Map;
bitset <maxn> is_prime;

int read()
{
    int x = 0, k = 1
 
;
    char c;
    c = getchar();
    while(c < ‘0‘ || c > ‘9‘) { if(c == ‘-‘) k = -1; c = getchar(); }
    while(c >= ‘0‘ && c <= ‘9‘) x = x * 10 + c - ‘0‘, c = getchar();
    return x * k;
}

int qpow(int x, int times)
{
    int base = 1; x %= mod;
    for(; times; times >>= 1
 
, x = (x * x) % mod)
        if(times & 1) base = (base * x) % mod;
    return base;
}

void Get_Mu()
{
    Sum[1] = 1;
    for(int i = 2; i <= maxx; i ++)
    {
        if(!is_prime[i]) pri[++ tot] = i, Sum[i] = -1;
        for(int j = 1; j <= tot; j ++)
        {
            int tem = i * pri[j];
            if(tem > maxx) break;
            is_prime[tem] = 1;
            if(!(i % pri[j])) { Sum[tem] = 0; break; }
            else Sum[tem] = - Sum[i];
        }
    }
    for(int i = 1; i <= maxx; i ++) Sum[i] = (Sum[i] + Sum[i - 1]) % mod;
}

int Mu(int x)
{
    if(x <= maxx) return Sum[x];
    if(Map[x]) return Map[x];
    int ret = 0;
    for(int l = 2, r; l <= x; l = r + 1)
    {
        r = x / (x / l);
        ret = (ret + (r - (l - 1)) * Mu(x / l) % mod) % mod; 
    }
    return Map[x] = (1 - ret + mod) % mod;
}

int Solve(int n, int m)
{
    int ret = 0;
    for(int l = 1, r; l <= m; l = r + 1)
    {
        if(n / l) r = min((n / (n / l)), (m / (m / l)));
        else r = (m / (m / l));
        ret += qpow(m / l - n / l, N) % mod * (Mu(r) - Mu(l - 1)) % mod;
        ret %= mod;
    }
    return ret;
}

signed main()
{
    N = read(), K = read(), L = read(), H = read(), ans = 0;
    Get_Mu();
    int l = floor((db) (L - 1) / (db) K), r = floor((db) H / (db) K);
    ans = Solve(l, r);
    printf("%lld\n", (ans + mod) % mod);
    return 0;
}

【題解】CQOI2015選數