「CQOI2015」选数

「CQOI2015」选数

求在区间$[L, R]$内选取$n$个数(可重复)，最大公约数为$k$的方案。答案对$10^9+7$取模。

$n, k <= 10^9$
$L, R <= 10^9$

---

我们发现如果要使选出的数的最大公约数为$k$，那选出的这些数必定都是$k$的倍数。
再观察一下，我们会发现这些$k$的倍数的最大公约数也只可能是$k$的倍数。

所以我们设
$f(x) = [L, R]$内选$n$个数, gcd为$k$的方案数
$g(x) = [L, R]$内选$n$个数，gcd为$k$的倍数的方案数

$[L, R]$内$k$的倍数的个数为

$$\lfloor \dfrac{R}{k} \rfloor - \lfloor \dfrac{L-1}{k}\rfloor$$
那么

$$g(x) =\bigg(\lfloor \dfrac{R}{k} \rfloor - \lfloor \dfrac{L-1}{k}\rfloor\bigg)^n$$

$$g(x) = \sum_{x|d} f(d)$$
$$\Rightarrow f(x) = \sum_{x|d} \mu(\dfrac{d}{x})g(d)$$
$$\therefore f(x) = \sum_{x|d} \mu(\dfrac{d}{x})\bigg(\lfloor \dfrac{R}{d} \rfloor - \lfloor \dfrac{L-1}{d}\rfloor\bigg)^n$$

<script type="math/tex; mode=display" id="MathJax-Element-257"></script>

<script type="math/tex; mode=display" id="MathJax-Element-258"></script>

我们要求的就是$f(k)$

$$ans = \sum_{x|k} \mu(\dfrac{d}{x})\bigg(\lfloor \dfrac{R}{k} \rfloor - \lfloor \dfrac{L-1}{k}\rfloor\bigg)^n$$

$$ans = \sum_{x=1} \mu(x) \bigg(\lfloor \dfrac{R}{xk} \rfloor - \lfloor \dfrac{L-1}{xk}\rfloor\bigg)^n$$

<script type="math/tex; mode=display" id="MathJax-Element-944"></script>

<script type="math/tex; mode=display" id="MathJax-Element-945"></script>

我们发现，当$\lfloor \dfrac{R}{xk} \rfloor,\lfloor \dfrac{L-1}{xk}\rfloor$相同时
$\bigg(\lfloor \dfrac{R}{xk} \rfloor - \lfloor \dfrac{L-1}{xk}\rfloor\bigg)^n$也相同

所以我们可以在$O(\sqrt{n}logn)$内求出ans 的值，用杜教筛求出$\mu(x)$的前缀和即可。

要注意取模的数作差的时候要加上模数！
$(a - b) \% MOD \Rightarrow (a - b + MOD) \% MOD$

code：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD = 1000000007;
const int N = 5000000;
const int LIM = 1000000;

ll ksm( ll a, int b )
{
    ll ret = 1; 
    while( b ) 
    {
        if( b & 1 ) ret = (ret * a) %MOD;
        a = (a * a) %MOD; b >>= 1; 
    }
    return ret;
}

map <int, ll> m;

bool is_prime[N];
int tot = 0, prime[N], mobius[N];

int n, k, L, R;

void init( int size )
{
    memset( is_prime, true, sizeof( is_prime ) );
    is_prime[1] = 1; mobius[1] = 1;

    for( int i = 2; i <= size; i ++ )
    {
        if( is_prime[i] )
            prime[++tot] = i, mobius[i] = -1;

        for( int j = 1; j <= tot; j ++ )
        {
            int k = i * prime[j];
            if( k > size ) break;
            is_prime[k] = false;

            if( i % prime[j] == 0 )
                mobius[k] = 0;
            else mobius[k] = - mobius[i];   
        }
    }

    for( int i = 2; i <= size; i ++ )
        mobius[i] = (mobius[i] + mobius[i-1]) %MOD;
}

ll calc_mobius( int x )
{
    if( x <= LIM ) return mobius[x];
    if( m[x] ) return m[x];

    int next = 0; ll ret = 0;
    for( int i = 2; i <= x; i = next + 1 )
    {
        next = x / ( x / i );
        ret -= calc_mobius(x/i) * (next-i+1);
    }

    m[x] = ret %MOD;
    return ret %MOD;
}


int main()
{
    init( LIM ); 

    scanf( "%d%d%d%d", &n, &k, &L, &R );

    R = R / k; L = ( L - 1 ) / k;

    int next = 0; ll ans = 0;
    for( int i = 1; i <= R; i = next + 1 )   
    {
//      printf( "i:%d, L:%d, R:%d\n", i, L, R );
        if( L/i ) next = min( R/(R/i), L/(L/i) ); else next = R / (R/i);
        ans += ( ksm( (R/i)-(L/i), n ) * ( calc_mobius(next) - calc_mobius(i-1) + MOD ) %MOD);
        ans %= MOD;
//      printf( "next:%d, mobsum(next):%d\n", next, calc_mobius(next) );
    }

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

    return 0;
}

// 2 2 2 4