[CQOI2015]选数（数论分块+杜教筛）

problem

洛谷链接

solution

将 $L,H$ 的范围放缩 $\frac{1}{K}$，都除掉 $K$，特殊的 $L$ 边界注意一下。

$H\leftarrow H/K,L\leftarrow (L-1)/K+1$。

问题转化为 $[L,H]$ 中任选 $N$ 个数 $\text{gcd}=1$ 的方案数。

$T(i):$ 从 $[L,H]$ 中选 $N$ 个数 $i|\text{gcd}$ 的方案数，即 $(\lfloor \frac{H}{i}\rfloor -\lfloor \frac{L-1}{i}\rfloor {)}^N$。

$t(i):$ 从 $[L,H]$ 中选 $N$ 个数 $i=\text{gcd}$ 的方案数。

根据定义显然有，$T(i)={\sum }_{i|d}t(d)$

莫比乌斯反演得 $t(i)={\sum }_{i|d}\mu (\frac{d}{i})T(d)$

答案即为 $t(1)={\sum }_{i=1}^{inf}\mu (i)T(i)={\sum }_{i=1}^{inf}\mu (i)(\lfloor \frac{H}{i}\rfloor -\lfloor \frac{L-1}{i}\rfloor {)}^N$。

$H,L$ 非常大，不能直接线筛后整除分块。

但可杜教筛，设定一个阀值 $M$ 预处理出 $M$ 以内的 $\mu$，同样整除分块后杜教筛能做到 $O(H^{\frac{2}{3}})$。

按 $T(i)$ 分类，假设 $i\in [l,r]$ 的 $T(i)$ 都是一样的，那么对 ${\sum }_{i=l}^r\mu (i)$ 进行杜教筛迅速求和。

${\sum }_{i=l}^r\mu (i)={\sum }_{i=1}^r\mu (i)-{\sum }_{i=1}^{l-1}\mu (i)$。这样就化成了标准的杜教筛形式。

$\sum \mu (i):ϵ=\mu \ast I$。$h\leftrightarrow ϵ;f\leftrightarrow \mu ;g\leftrightarrow I$。

$s(n)={\sum }_{i=1}^{n}f(i)={\sum }_{i=1}^n\mu (i)$

$g(1)s(n)={\sum }_{i=1}^nϵ(i)-{\sum }_{i=2}^{n}g(i)s(\lfloor \frac{n}{i}\rfloor )\iff s(n)=1-{\sum }_{i=2}^{n}s(\lfloor \frac{n}{i}\rfloor )$。

对 $s$ 函数进行记忆化递归，以及同样的整除分块。

这样子连 $H-L\le 1e5$ 的性质都没有用上！^_^ 这性质好像是拿来容斥递推用的。

一些省掉的推导☞莫比乌斯反演，杜教筛。

code

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define mod 1000000007
#define maxn 1000005
int mu[maxn], prime[maxn];
bool vis[maxn];
unordered_map < int, int > mp;
int N, M, K, L, H, cnt;

int qkpow( int x, int y ) {
    int ans = 1;
    while( y ) {
        if( y & 1 ) ans = ans * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ans;
}

int solve( int n ) {
    if( n <= M ) return mu[n];
    if( mp.find( n ) != mp.end() ) return mp[n];
    int ans = 1;
    for( int l = 2, r;l <= n;l = r + 1 ) {
        r = n / ( n / l );
        ( ans -= solve( n / l ) * ( r - l + 1 ) ) %= mod;
    }
    return mp[n] = ans;
}

signed main() {
    scanf( "%lld %lld %lld %lld", &N, &K, &L, &H );
    H /= K, L = ( L - 1 ) / K + 1;
    M = min( (int)1e6, H );
    mu[1] = 1; 
    for( int i = 2;i <= M;i ++ ) {
        if( ! vis[i] ) prime[++ cnt] = i, mu[i] = -1;
        for( int j = 1;j <= cnt and i * prime[j] <= M;j ++ ) {
            vis[i * prime[j]] = 1;
            if( i % prime[j] == 0 ) { mu[i * prime[j]] = 0; break; }
            else mu[i * prime[j]] = -mu[i];
        }
    }
    for( int i = 1;i <= M;i ++ ) mu[i] += mu[i - 1];
    L --; int ans = 0;
    for( int l = 1, r;l <= H;l = r + 1 ) {
        if( ! ( L / l ) ) r = H / ( H / l );
        else r = min( H / ( H / l ), L / ( L / l ) );
        ( ans += qkpow( H / l - L / l, N ) * ( solve( r ) - solve( l - 1 ) ) ) %= mod;
    }
    printf( "%lld\n", ( ans + mod ) % mod );
    return 0;
}