[BZOJ3930][CQOI2015]选数（莫比乌斯反演+杜教筛）-优快云博客

=== ===

题解

设 $F(k)$ 为在范围内选取n个数字，最大公因数为 $k$ 的方案数， $f(k)$ 为在范围内选取n个数字，公因数为 $k$ 的方案数。那么显然 $f(k)=\sum\limits_{k|d}F(d)$ ，并且 $f(k)=({\lfloor{H \over k}\rfloor}-{\lfloor{L-1 \over k}\rfloor})^n$ 。
由反演公式可得 $F(k)=\sum\limits_{k|d}\mu({d \over k})f(d)$ ，可以发现 $d$ 的枚举上限最大到 $H$ 并且必须是 $k$ 的倍数。那么设 $d=kd'$ ，用 $d'$ 替换 $d$ 可得 $F(k)=\sum\limits_{d=1}^{\lfloor{H \over k}\rfloor}\mu(d)f(kd)=\sum\limits_{d=1}^{\lfloor{H \over k}\rfloor}\mu(d)({\lfloor{H \over {kd}}\rfloor}-{\lfloor{L-1 \over {kd}}\rfloor})^n$
那么我们做一个预处理令 $H={\lfloor{H \over k}\rfloor}$ ， $L={\lfloor{L-1 \over k}\rfloor}$ ，于是 $F(k)=\sum\limits_{d=1}^H\mu(d)({\lfloor{H \over d}\rfloor}-{\lfloor{L \over d}\rfloor})^n$
写成这个样的式子就比较明确了，利用分块枚举除法+杜教筛就可以了。
有一个问题是因为这里枚举上限是到较大的 $H$ 而不是较小的 $L$ ，所以如果 $\lfloor{L \over d}\rfloor$ 变成0了就不要用它做除法了，需要特判一下

代码

#include<map>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=1000000;
const long long Mod=1000000007;
int n,k,L,H,prm[1000010],mu[1000010],tail;
bool ext[1000010];
long long ans;
map<int,int> rec;
long long powww(long long a,int t){
    long long ans=1;a%=Mod;
    while (t!=0){
        if (t&1) ans=(ans*a)%Mod;
        a=(a*a)%Mod;t>>=1;
    }
    return ans;
}
void get_prime(int N){
    mu[1]=1;
    for (int i=2;i<=N;i++){
        if (ext[i]==false){
            prm[++prm[0]]=i;mu[i]=-1;
        }
        for (int j=1;j<=prm[0];j++){
            if ((long long)i*prm[j]>N) break;
            ext[i*prm[j]]=true;
            if (i%prm[j]==0) break;
            else mu[i*prm[j]]=-mu[i];
        }
    }
    for (int i=1;i<=N;i++) mu[i]+=mu[i-1];
}
long long get_mu(int n){
    if (n<N) return mu[n];
    if (rec[n]!=0) return rec[n];
    long long sum=1;
    for (int i=2,tail;i<=n;i=tail+1){
        tail=n/(n/i);
        sum-=get_mu(n/i)*(long long)(tail-i+1)%Mod;
        sum%=Mod;
    }
    return rec[n]=sum;
}
int main()
{
    scanf("%d%d%d%d",&n,&k,&L,&H);
    L--;L/=k;H/=k;
    get_prime(N);
    for (int i=1;i<=H;i=tail+1){
        if (L/i!=0) tail=min(H/(H/i),L/(L/i));
        else tail=H/(H/i);
        ans=(ans+powww((H/i-L/i),n)*(get_mu(tail)-get_mu(i-1))%Mod)%Mod;
    }
    ans=(ans+Mod)%Mod;
    printf("%I64d\n",ans);
    return 0;
}