[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;
}

偏偏在最后出现的补充说明

这题的反演做法好像并没有用到$H-L \leq 10^5$这个限制？但是因为区间范围比较小似乎有一种代码长度短一些的容斥做法。。