标签：find ... check code red pre show line 使用
从[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];
}
标签：find ... check code red pre show line 使用
原文地址：https://www.cnblogs.com/a1b3c7d9/p/10795462.html
