【题解】CQOI2015选数

　　这题做的时候接连想错了好多次……但是回到正轨上之后依然是一个套路题。（不过这题好像有比莫比乌斯反演更好的做法，莫比乌斯反演貌似是某种能过的暴力ヽ（´ー｀）┌）不过能过也就行了吧哈哈。

　　首先我们把数字的范围要进行缩小：最大公约数为 \(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 + ;

int N, K, L, H, ans, Sum[maxn];

int tot, pri[maxn];

map <int, int> Map;

bitset <maxn> is_prime;

int read()

{

    int x = , k = ;

    char c;

    c = getchar();

    while(c < '' || c > '') { if(c == '-') k = -; c = getchar(); }

    while(c >= '' && c <= '') x = x *  + c - '', c = getchar();

    return x * k;

}

int qpow(int x, int times)

{

    int base = ; x %= mod;

    for(; times; times >>= , x = (x * x) % mod)

        if(times & ) base = (base * x) % mod;

    return base;

}

void Get_Mu()

{

    Sum[] = ;

    for(int i = ; i <= maxx; i ++)

    {

        if(!is_prime[i]) pri[++ tot] = i, Sum[i] = -;

        for(int j = ; j <= tot; j ++)

        {

            int tem = i * pri[j];

            if(tem > maxx) break;

            is_prime[tem] = ;

            if(!(i % pri[j])) { Sum[tem] = ; break; }

            else Sum[tem] = - Sum[i];

        }

    }

    for(int i = ; i <= maxx; i ++) Sum[i] = (Sum[i] + Sum[i - ]) % mod;

}

int Mu(int x)

{

    if(x <= maxx) return Sum[x];

    if(Map[x]) return Map[x];

    int ret = ;

    for(int l = , r; l <= x; l = r + )

    {

        r = x / (x / l);

        ret = (ret + (r - (l - )) * Mu(x / l) % mod) % mod;

    }

    return Map[x] = ( - ret + mod) % mod;

}

int Solve(int n, int m)

{

    int ret = ;

    for(int l = , r; l <= m; l = r + )

    {

        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 - )) % mod;

        ret %= mod;

    }

    return ret;

}

signed main()

{

    N = read(), K = read(), L = read(), H = read(), ans = ;

    Get_Mu();

    int l = floor((db) (L - ) / (db) K), r = floor((db) H / (db) K);

    ans = Solve(l, r);

    printf("%lld\n", (ans + mod) % mod);

    return ;

}