[CQOI2015]选数

Posted 2021-11-05 a1b3c7d9

[CQOI2015]选数

从[L,H]中可重复地选出N个数，问其gcd等于K的方案\(mod\ 10^9+7\)，1<=N,K<=10^9，1<=L<=H<=10^9，H-L<=10^5。

解

法一：直接

转换为Mobius问题，不难有

\[ans=\sum_{i_1=L}^H\sum_{i_2=L}^H...\sum_{i_N=L}^H(gcd(i_1,i_2,...,i_N)==k)\]

于是设

\[f(k)=\sum_{i_1=L}^H\sum_{i_2=L}^H...\sum_{i_n=L}^H(gcd(i_1,i_2,...,i_n)==k)\]
\[F(k)=\sum_{i_1=L}^H\sum_{i_2=L}^H...\sum_{i_N=L}^H(k|gcd(i_1,i_2,...,i_N))=([\frac{H}{k}]-[\frac{L-1}{k}])^N\]

由Mobius反演定理代入有

\[ans=\sum_{k|d}([\frac{H}{d}]-[\frac{L-1}{d}])^N\mu(d/k)=\]
\[\sum_{d=1}^{[N/k]}([\frac{H}{dk}]-[\frac{L-1}{dk}])^N\mu(d)\]

显然左式可以整除分块，而至于后式只需要快速维护其前缀和，注意到N很大，于是用杜教筛优化即可。

参考代码：

#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define ll long long
#define control 6000000
#define yyb 1000000007
using namespace std;
template<class f1,class f2>
struct chain{
    chain*next;f1 x;f2 y;
};
template<class f1,class f2,class f3,f3 key>
struct unordered_map{
    chain<f1,f2>*head[key],*p;
    il void insert(f1 x,f2 y){
        p=new chain<f1,f2>,p->x=x,p->y=y;
        x%=key,p->next=head[x],head[x]=p;
    }
    il chain<f1,f2>* find(f1 x){
        for(p=head[x%key];p!=NULL;p=p->next)
            if(p->x==x)return p;
        return NULL;
    }
};
bool check[control+1];
int prime[600000],pt,mb[control+1];
il void prepare(int);
il int djs(int),pow(int,int),min(int,int);
int main(){
    int n,k,l,h,i,j,ans(0);
    scanf("%d%d%d%d",&n,&k,&l,&h);
    h/=k,(--l)/=k,prepare(control);
    for(i=1;i<=h;i=j+1){
        j=h/(h/i);if(l/i)j=min(j,l/(l/i));
        (ans+=(ll)(djs(j)-djs(i-1))*pow(h/i-l/i,n)%yyb)%=yyb;
    }printf("%d",(ans+yyb)%yyb);
    return 0;
}
il int min(int a,int b){
    return a<b?a:b;
}
il int pow(int x,int y){
    int ans(1);while(y){
        if(y&1)ans=(ll)ans*x%yyb;
        x=(ll)x*x%yyb,y>>=1;
    }return ans;
}
unordered_map<int,int,int,999983>m;
il int djs(int n){
    if(n<=control)return mb[n];
    if(m.find(n)!=NULL)return m.find(n)->y;int i,j,ans(1);
    for(i=2;i<=n;i=j+1)j=n/(n/i),ans-=(j-i+1)*djs(n/j);
    return m.insert(n,ans),ans;
}
il void prepare(int n){
    int i,j;check[1]|=mb[1]|=true;
    for(i=2;i<=n;++i){
        if(!check[i])prime[++pt]=i,mb[i]=-1;
        for(j=1;j<=pt&&prime[j]<=n/i;++j){
            check[i*prime[j]]|=true;
            if(!(i%prime[j]))break;
            mb[i*prime[j]]=~mb[i]+1;
        }
    }for(i=1;i<=n;++i)mb[i]+=mb[i-1];
}

法二：优化

思路一：

---

\[gcd(a,b)\leq|a-b|(a\neq b)\]

证明：
\(|x-y|=|k_1gcd(x,y)-k_2gcd(x,y)|=gcd(x,y)|k_1-k_2|\)

故得证。

---

注意到\(H-L\leq10^5\),于是猜测gcd范围应在里面，于是我们有结论\(gcd(a,b)\leq|a-b|(a\neq b)\),所以只要对a,b相同的情况单独讨论，而显然这种方案有且仅有一种，只要看其是否落在\([L,H]\)内即可，于是接下来就可以直接对法一进行优化了,但是注意该结论一定是不存在全相同对数。

思路二：

而注意到问题比较简单，而有可以使用Mobius反演，于是我们考虑容斥，于是令\(L=[L-1]/k+1,H/=k\),接下来问题变成只要对互质对数考虑即可，所以设\(D(n)=([\frac{H}{n}]-[\frac{L-1}{n}])^N-[\frac{H}{n}]+[\frac{L-1}{n}]\)，表示有公约数n的互质的不全相同的对数，根据容斥原理，我们有

\[ans=\sum_{d=1}^{|H-L|}\mu(d)D(d)\]

进行整除分块即可，注意后面全相同对数的情况，因为这也是一种方案。

参考代码：

#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define ll long long
#define yyb 1000000007
using namespace std;
bool check[100001];
int prime[2500],pt,mb[100001];
il void prepare(int);
il int pow(int,int),min(int,int);
int main(){
    int n,k,l,h,l_,i,j,cjx,li,ans(0);
    scanf("%d%d%d%d",&n,&k,&l,&h);
    ++(--l/=k),h/=k,l_=l-1;
    li=h-l,prepare(li);
    for(i=1;i<=li;i=cjx+1){
        j=h/i-l_/i,cjx=h/(h/i);
        if(l_/i)cjx=min(cjx,l_/(l_/i));
        (ans+=(ll)(mb[cjx]-mb[i-1])*(pow(j,n)-j)%yyb)%=yyb;
    }if(l==1)++ans;printf("%d",(ans+yyb)%yyb);
    return 0;
}
il int min(int a,int b){
    return a<b?a:b;
}
il int pow(int x,int y){
    int ans(1);while(y){
        if(y&1)ans=(ll)ans*x%yyb;
        x=(ll)x*x%yyb,y>>=1;
    }return ans;
}
il void prepare(int n){
    int i,j;check[1]|=mb[1]|=true;
    for(i=2;i<=n;++i){
        if(!check[i])prime[++pt]=i,mb[i]=-1;
        for(j=1;j<=pt&&prime[j]<=n/i;++j){
            check[i*prime[j]]|=true;
            if(!(i%prime[j]))break;
            mb[i*prime[j]]=~mb[i]+1;
        }
    }for(i=1;i<=n;++i)mb[i]+=mb[i-1];
}