[CQOI2015]选数

学了杜教筛后做的第一道题。

传送门

求 $\sum_{i_1=l}^h…\sum_{i_n=l}^h\ [gcd(i_1…i_n)=k]$

为了看起来好看偷懒，下面把$\sum_{i_1=l}^{h}…\sum_{i_n=l}^{h}$这个东西写成n=2的情况：$\sum_{i=l}^h\sum_{j=l}^k$

先提出k，枚举k的倍数
$$
\sum_{i=(l-1)/k+1}^{h/k}\sum_{j=(l-1)/k}^{h/k}[gcd(i,j)=1]
$$
套路的反演一下
$$
\sum_{i=(l-1)/k}^{h/k}\sum_{j=(l-1)/k}^{h/k}\sum_{d|i,d|j}\mu(d)
$$
改为枚举d和他的倍数
$$
\sum_{d=1}^{h/k}\mu(d)\sum_{i=\frac{l-1}{kd}}^{h/kd}\sum_{j=\frac{l-1}{kd}}^{h/kd}
$$
这时候把后面的改成应有的形式：$(\frac{h}{kd}-\frac{l-1}{kd})^n$

这样只要求出$\mu$的前缀和就好了，杜教筛盘他。

#include<map>
#include<cstdio>
#include<iostream>
#define ll long long
#define mod 1000000007
using namespace std;
int o[2000001],prime[1000010],sum[2000001];
int p,k,n,m,blo,mu[2000001];
ll ans;
inline ll ksm(ll a,ll b){
    ll num=1;
    while (b){
        if (b&1) num=num*a%mod;
        a=a*a%mod; b>>=1;
    }
    return num;
}
inline void euler(){
    o[1]=mu[1]=1;
    for (int i=2; i<=blo; i++){
        if (!o[i]){
            o[i]=1; prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for (int j=1; j<=prime[0]&&i*prime[j]<=blo; j++){
            o[i*prime[j]]=1;
            if (i%prime[j]==0) break;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for (int i=1; i<=blo; i++) sum[i]=sum[i-1]+mu[i];
}
map<int,int> nub;
int get_mu(int x){
    if (x<=blo) return sum[x];
    if (nub[x]) return nub[x];
    int num=1;
    for (int l=2,r; l<=x; l=r+1){
        r=x/(x/l);
        num-=(r-l+1)*get_mu(x/l);
    }
    return nub[x]=num;
    // return num;
}
int main(){
    scanf("%d%d%d%d",&p,&k,&n,&m);
    n=(n-1)/k; m=m/k;// if (!n) n++;
    blo=2e6; euler(); //printf("%d\n",blo);
    // puts("ok");
    for (int l=1,r; l<=m; l=r+1){
        if (n/l==0) r=m/(m/l);
        else r=min(n/(n/l),m/(m/l));
        // r=min(n/l?n/(n/l):(int)2e9,m/(m/l));
        ans=(ans+1LL*(get_mu(r)-get_mu(l-1))*ksm(m/l-n/l,p)%mod)%mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}