[CQOI2015]选数

Posted 2021-02-03 wolfycz

Description
我们知道，从区间[L,H]（L和H为整数）中选取N个整数，总共有(H-L+1)^N种方案。小z很好奇这样选出的数的最大公约数的规律，他决定对每种方案选出的N个整数都求一次最大公约数，以便进一步研究。然而他很快发现工作量太大了，于是向你寻求帮助。你的任务很简单，小z会告诉你一个整数K，你需要回答他最大公约数刚好为K的选取方案有多少个。由于方案数较大，你只需要输出其除以1000000007的余数即可。

Input
输入一行，包含4个空格分开的正整数，依次为N，K，L和H。

Output
输出一个整数，为所求方案数。

Sample Input
2 2 2 4

Sample Output
3

HINT
(1leqslant N,Kleqslant10^9,1leqslant Lleqslant Hleqslant10^9,H-Lleqslant10^5)

题目要求（令(H)变为(R)）

[sumlimits_{a_1=L}^Rsumlimits_{a_2=L}^R...sumlimits_{a_n=L}^R[gcdlimits_{i=1}^n(a_i)=K]]

同样的，我们令所有数整除一个(K)，注意到(Lleqslant a_ileqslant R)，那么再次枚举的(a_i')要满足(Lleqslant a_i'×Kleqslant R)，所以有

[sumlimits_{a_1=lfloorfrac{L-1}{d} floor+1}^{lfloorfrac{R}{d} floor}...sumlimits_{a_n=lfloorfrac{L-1}{d} floor+1}^{lfloorfrac{R}{d} floor}[gcdlimits_{i=1}^n(a_i)=1]]

对最后一部分反演，并且将枚举(x)提前，得到

[sumlimits_{x=1}^{lfloorfrac{R}{d} floor}mu(x)sumlimits_{a_1=lfloorfrac{L-1}{d} floor+1}^{lfloorfrac{R}{d} floor}...sumlimits_{a_n=lfloorfrac{L-1}{d} floor+1}^{lfloorfrac{R}{d} floor}1]

后面那一坨长的都一个样子，而且可以给他们容斥一下，得到

[sumlimits_{x=1}^{lfloorfrac{R}{d} floor}mu(x)(lfloordfrac{R}{d} floor-lfloordfrac{L-1}{d} floor)^n]

如果我们能够求出(sumlimits_{i=1}^nmu(i))，那么我们就可以进行整除分块了

线筛？死了……(nleqslant 10^9)；于是我们杜教筛一波，可以参考浅谈算法——杜教筛(这篇文章还没写，处于咕咕咕的状态……)

然后这题就做完了

/*program from Wolfycz*/
#include<map>
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define Fi first
#define Se second
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline char gc(){
    static char buf[1000000],*p1=buf,*p2=buf;
    return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
inline int frd(){
    int x=0,f=1; char ch=gc();
    for (;ch<'0'||ch>'9';ch=gc())   if (ch=='-')    f=-1;
    for (;ch>='0'&&ch<='9';ch=gc()) x=(x<<3)+(x<<1)+ch-'0';
    return x*f;
}
inline int read(){
    int x=0,f=1; char ch=getchar();
    for (;ch<'0'||ch>'9';ch=getchar())  if (ch=='-')    f=-1;
    for (;ch>='0'&&ch<='9';ch=getchar())    x=(x<<3)+(x<<1)+ch-'0';
    return x*f;
}
inline void print(int x){
    if (x<0)    putchar('-'),x=-x;
    if (x>9)    print(x/10);
    putchar(x%10+'0');
}
const int N=1e6,p=1e9+7;
map<int,int>Sum;
int prime[N+10],miu[N+10],sum[N+10];
bool inprime[N+10];
void prepare(){
    miu[1]=sum[1]=1; int tot=0;
    for (int i=2;i<=N;i++){
        if (!inprime[i])    miu[prime[++tot]=i]=-1;
        for (int j=1;j<=tot&&i*prime[j]<=N;j++){
            inprime[i*prime[j]]=1;
            if (i%prime[j]==0){
                miu[i*prime[j]]=0;
                break;
            }miu[i*prime[j]]=-miu[i];
        }
        sum[i]=sum[i-1]+miu[i];
    }
}
int S(int n){
    if (n<=N)   return sum[n];
    map<int,int>::iterator it=Sum.find(n);
    if (it!=Sum.end())  return it->Se;
    int res=1;
    for (int i=2,pos;i<=n;i=pos+1){
        pos=n/(n/i);
        res=(res-1ll*(pos-i+1)*S(n/i)%p)%p;
    }
    res=(res+p)%p;
    Sum.insert(map<int,int>::value_type(n,res));
    return res;
}
int mlt(int a,int b){
    int res=1;
    for (;b;b>>=1,a=1ll*a*a%p)  if (b&1)    res=1ll*res*a%p;
    return res;
}
int main(){
    prepare();
    int n=read(),k=read(),L=read(),R=read(),Ans=0;
    --L/=k,R/=k;
    for (int i=1,pos;i<=R;i=pos+1){
        pos=min(L/i?L/(L/i):inf,R/(R/i));
        Ans=(Ans+1ll*(S(pos)-S(i-1))*mlt(R/i-L/i,n)%p)%p;
    }
    printf("%d
",(Ans+p)%p);
    return 0;
}