[BZOJ 3930] [CQOI 2015]选数(莫比乌斯反演+杜教筛)

[BZOJ 3930] [CQOI 2015]选数(莫比乌斯反演+杜教筛)

题面

我们知道，从区间([L,R])（L和R为整数）中选取N个整数，总共有((R-L+1)^N)种方案。求最大公约数刚好为K的选取方案有多少个。由于方案数较大，你只需要输出其除以1000000007的余数即可。

[N,K,L,H leq 10^9,H-L leq 10^5 ]

分析

(ecause gcd(ka,kb)=kgcd(a,b)),我们先把(L,R)除以(K),然后问题就变成了求gcd=1的方案数

设(f(x))表示区间[l,r]里选n个数，gcd为x的方案数

设(F(x))表示区间[l,r]里选n个数，gcd被x整除的方案数

(ecause x|gcd(i,j), herefore x|i,x|j)

[l,r]里被x整除的数有((lfloor frac{r}{x} floor-lfloor frac{l-1}{x} floor))个

因此(F(x)=(lfloor frac{r}{x} floor-lfloor frac{l-1}{x} floor)^n)

(F,f)显然满足莫比乌斯反演的第二种形式，(F(x)=sum_{d|x} f(d))

(f(x)=sum_{x|d} F(d) mu(frac{d}{x}))

我们要求的是

[f(1)=sum_{1|d} F(d) mu(d)=sum_{d=1}^r mu(d) (lfloor frac{r}{d} floor-lfloor frac{l-1}{d} floor)^n ]

后面的部分可以数论分块然后快速幂求解，但由于(r leq 10^9)，不能直接线性筛(mu)的前缀和，需要用杜教筛。

---

杜教筛：

套路公式:

我们要求(f)的前缀和,构造两个函数(g,h)满足(h=f*g), (F,G,H)为它们的前缀和

[g(1)F(n)=H(n)-sum_{d=2}^n g(d) F(frac{n}{d}) ]

如果(f=mu),注意到(mu*I=varepsilon),那么(g(n)=I(n)=1,h(n)=varepsilon(n),H(n)=varepsilon(1)=1)

代入得(F(n)=1-sum_{d=2}^n F(frac{n}{d}))

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<map>
#define maxn 2000000
#define mod 1000000007
using namespace std;
typedef long long ll;
int n,k;
ll A,B;

int cnt;
bool vis[maxn+5];
int prime[maxn+5];
int mu[maxn+5];
ll s_mu[maxn+5];
void sieve(int n){
	mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]){
			prime[++cnt]=i;
			mu[i]=-1;
		} 
		for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
			vis[i*prime[j]]=1;
			if(i%prime[j]==0){
				mu[i*prime[j]]=0;
				break;
			}else mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<=n;i++) s_mu[i]=(s_mu[i-1]+mu[i])%mod;
}

map<ll,ll>sum_mu;
ll dujiao_sieve(ll x){
	if(x<=maxn) return s_mu[x];
	if(sum_mu.count(x)) return sum_mu[x];
	ll ans=1;
	for(ll l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		ans-=(r-l+1)*dujiao_sieve(x/l)%mod;
		ans=(ans+mod)%mod;
	}
	sum_mu[x]=ans;
	return ans;
} 

inline ll fast_pow(ll x,ll k){
	ll ans=1;
	while(k){
		if(k&1) ans=ans*x%mod;
		x=x*x%mod;
		k>>=1;
	}
	return ans;
}

int main(){
	sieve(maxn);
	scanf("%d %d %lld %lld",&n,&k,&A,&B);
	A=(A-1)/k;
	B/=k;
	ll ans=0;
	for(ll l=1,r;l<=B;l=r+1){
		if(A/l) r=min(A/(A/l),B/(B/l));
		else r=B/(B/l);
//		printf("%d %d
",l,r);
		ans+=fast_pow(B/l-A/l,n)*(dujiao_sieve(r)-dujiao_sieve(l-1)+mod)%mod;
		ans%=mod;
	}
	printf("%lld
",ans);
}