P3172 [CQOI2015] 选数

题意

在 \([L,H]\) 中去 \(n\) 个数，问祂们的 \(\gcd\) 为 \(k\) 的方案数。
\(n,k\le10^9,L\le H\le10^9,H-L\le10^5\)。

思路

设 \(f(i)\) 表示选 \(n\) 个数，有 \(i\) 这个公约数的方案数，\(F(i)\) 表示选 \(n\) 个数，\(GCD\) 是 \(i\) 的方案数。有 \(f(i)=\sum_{i|d}F(d)\)，反演，变成 \(F_d=\sum_{d|i}\mu(\frac id)f(i)\)。显然， \(f(d)=(\lfloor\frac Hd\rfloor-\lfloor\frac {L-1}d\rfloor)^n\)。
考虑用整除分块优化，并使用杜教筛求出 \(\mu\) 的前缀和。
我的杜教筛学习笔记

代码

#include<bits/stdc++.h>
using namespace std;
namespace IO{
    template<typename T>
    inline void read(T&x){
        x=0;char c=getchar();bool f=0;
        while(!isdigit(c)) c=='-'?f=1:0,c=getchar();
        while(isdigit(c)) x=x*10+c-'0',c=getchar();
        f?x=-x:0;
    }
    template<typename T>
    inline void write(T x){
        if(x==0){putchar('0');return ;}
        x<0?x=-x,putchar('-'):0;short st[50],top=0;
        while(x) st[++top]=x%10,x/=10;
        while(top) putchar(st[top--]+'0');
    }
    inline void read(char&c){c=getchar();while(isspace(c)) c=getchar();}
    inline void write(char c){putchar(c);}
    inline void read(string&s){s.clear();char c;read(c);while(!isspace(c)&&~c) s+=c,c=getchar();}
    inline void write(string s){for(int i=0,len=s.size();i<len;i++) putchar(s[i]);}
    template<typename T>inline void write(T*x){while(*x) putchar(*(x++));}
    template<typename T,typename...T2> inline void read(T&x,T2&...y){read(x),read(y...);}
    template<typename T,typename...T2> inline void write(const T x,const T2...y){write(x),putchar(' '),write(y...),sizeof...(y)==1?putchar('\n'):0;}
}using namespace IO;
template<int mod>struct Modint{
    int z;
    Modint(){z=0;}
    Modint(int x){x%=mod;z=x<0?x+mod:x;}
    Modint(long long x){x%=mod;z=x<0?x+mod:x;}
    Modint(short x){x%=mod;z=x<0?x+mod:x;}
    Modint(char x){x%=mod;z=x<0?x+mod:x;}
    Modint(bool x){x%=mod;z=x<0?x+mod:x;}
    friend Modint operator+(Modint t,Modint t2){Modint ans;ans.z=(t.z+t2.z)%mod;return ans;}
    friend Modint operator*(Modint t,Modint t2){Modint ans;ans.z=1ll*t.z*t2.z%mod;return ans;}
    friend Modint operator-(Modint t,Modint t2){Modint ans;ans.z=(t.z-t2.z)%mod;return ans;}
    Modint operator-()const{return (Modint){-z};}
    Modint operator<<(const int t)const{Modint ans;ans.z=(z<<t)%mod;return ans;}
    Modint operator>>(const int t)const{Modint ans;ans.z=(z>>t)%mod;return ans;}
    Modint&operator+=(const Modint t){z=(z+t.z)%mod;return *this;}
    Modint&operator*=(const Modint t){z=1ll*z*t.z%mod;return *this;}
    Modint&operator-=(const Modint t){z=(z-t.z)%mod;return *this;}
    Modint&operator<<=(const int t){z=(z<<t)%mod;return *this;}
    Modint&operator>>=(const int t){z=(z>>t)%mod;return *this;}
    Modint&operator++(){z++,z%=mod;return *this;}
    Modint&operator--(){z--,z%=mod;return *this;}
    Modint operator++(int){Modint ls=*this;z++,z%=mod;return ls;}
    Modint operator--(int){Modint ls=*this;z--,z%=mod;return ls;}
    friend Modint ksm(Modint a,int b){
        Modint ans=1;
        while(b){if(b&1) ans=ans*a;a=a*a,b>>=1;}
        return ans;
    }
    friend void read(Modint&z){
        int x=0;char c=getchar();bool f=0;
        while(!isdigit(c)) c=='-'?f=1:0,c=getchar();
        while(isdigit(c)) x=(x*10ll+c-'0')%mod,c=getchar();
        f?x=-x:0;
        z.z=x;
    }
    friend void write(Modint x){x.z<0?x.z+=mod:0;write(x.z);}
};
const int maxn=1000010,mod=1000000007;
#define M Modint<mod>
int mu[maxn],sum[maxn],prim[maxn],n,k,H,L;
map<int,M>mp;
bool flag[maxn];
void init(int n){
    int cnt=0;
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(flag[i]==0) prim[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt;j++){
            if(prim[j]*i>n) break;
            flag[prim[j]*i]=1;
            if(i%prim[j]==0){mu[i*prim[j]]=0;break;}
            mu[i*prim[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;i++) sum[i]=sum[i-1]+mu[i];
}
M calc(int n){
    if(n<=1000000) return sum[n];
    if(mp.find(n)!=mp.end()) return mp[n];
    M ans=0;
    for(int i=2,nxt;i<=n;i=nxt+1){
        nxt=n/(n/i);
        ans-=M(nxt-i+1)*calc(n/i);
    }
    return mp[n]=ans+1;
}
M calc_f(int d){return ksm(M(H/d-(L-1)/d),n);}
signed main(){
    init(1000000);
    read(n,k,L,H);
    M ans=0;
    L=(L-1)/k+1,H/=k;
    for(int i=1,nxt;i<=H;i=nxt+1){
        nxt=H/(H/i);
        if((L-1)/i) nxt=min(nxt,(L-1)/((L-1)/i));
        ans+=(calc(nxt)-calc(i-1))*calc_f(i);
    }
    write(ans);
    return 0;
}