【bzoj3930】[CQOI2015]选数 【莫比乌斯反演】【杜教筛】

题目传送门
题意：求从区间[L,H]$[L,H]$（L$L$和H$H$为整数）中选取N$N$个整数，使它们的gcd$gcd$为K$K$的方案总数模1000000007$1000000007$的值。
题解：我们令l=⌊L−1K⌋$l=\lfloor \frac{L-1}{K}\rfloor$，r=⌊HK⌋$r=\lfloor \frac{H}{K}\rfloor$。则原问题等价于求从区间[l,r]$[l,r]$中选取N$N$个整数，使它们的gcd$gcd$为1$1$的方案总数模1000000007$1000000007$的值。
我们令F(i)$F(i)$为选数的gcd为i的倍数的方案总数，则显然F(i)=(⌊ri⌋−⌊li⌋)N$F(i)=(\lfloor \frac{r}{i}\rfloor -\lfloor \frac{l}{i}\rfloor {)}^N$。
我们令f(i)$f(i)$为选数的gcd为i的方案总数，则显然F(i)=∑i|df(d)$F(i)=\sum _{i|d}f(d)$。
莫比乌斯反演可得f(i)=∑i|dμ(di)F(d)$f(i)=\sum _{i|d}\mu (\frac{d}{i})F(d)$。
把i=1带入，得f(1)=∑rd=1μ(d)(⌊rd⌋−⌊ld⌋)N$f(1)=\sum _{d=1}^r\mu (d)(\lfloor \frac{r}{d}\rfloor -\lfloor \frac{l}{d}\rfloor {)}^N$。
于是这个式子我们可以先杜教筛求出μ$\mu$的前缀和，再分块计算。
杜教筛是什么？
杜教筛点这里
代码：

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N=1000005,M=1000033;
const ll mod=1000000007;
int n,k,l,r,p[N];
ll ans,mu[N];
bool vis[N];
ll fastpow(ll a,ll x){
    a%=mod;
    ll res=1;
    while(x){
        if(x&1){
            res=res*a%mod;
        }
        x>>=1;
        a=a*a%mod;
    }
    return res;
}
struct Hashset{
    int cnt,head[M],to[M],nxt[M];
    ll v[M];
    bool count(int f){
        int x=f%M;
        for(int i=head[x];i;i=nxt[i]){
            if(to[i]==f){
                return true;
            }
        }
        return false;
    }
    ll get(int f){
        int x=f%M;
        for(int i=head[x];i;i=nxt[i]){
            if(to[i]==f){
                return v[i];
            }
        }
    }
    void set(int f,ll val){
        int x=f%M;
        to[++cnt]=f;
        nxt[cnt]=head[x];
        v[cnt]=val;
        head[x]=cnt;
    }
}s;
ll solve(int n){
    if(n<=1000000){
        return mu[n];
    }
    if(s.count(n)){
        return s.get(n);
    }
    ll res=1;
    for(int i=2,last;i<=n;i=last+1){
        last=n/(n/i);
        res-=solve(n/i)*(last-i+1)%mod;
        res%=mod;
    }
    s.set(n,res);
    return res;
}
int main(){
    scanf("%d%d%d%d",&n,&k,&l,&r);
    l=(l-1)/k;
    r/=k;
    if(!r){
        puts("0");
        return 0;
    }
    mu[1]=1;
    for(int i=2;i<=1000000;i++){
        if(!vis[i]){
            p[++p[0]]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=p[0]&&i*p[j]<=1000000;j++){
            vis[i*p[j]]=true;
            if(i%p[j]){
                mu[i*p[j]]=-mu[i];
            }else{
                break;
            }
        }
    }
    for(int i=2;i<=1000000;i++){
        mu[i]+=mu[i-1];
        mu[i]%=mod;
    }
    for(int i=1,last;i<=r;i=last+1){
        if(l/i){
            last=min(r/(r/i),l/(l/i));
        }else{
            last=r/(r/i);
        }
        ans+=fastpow(r/i-l/i,n)*(solve(last)-solve(i-1))%mod;
        ans%=mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}