BZOJ 3930 Luogu P3172 选数 (莫比乌斯反演)

手动博客搬家：本文发表于20180310 11:46:11, 原地址https://blog.csdn.net/suncongbo/article/details/79506484

题目链接:
(Luogu)https://www.luogu.org/problemnew/show/P3172
(BZOJ)http://www.lydsy.com/JudgeOnline/problem.php?id=3930

题目大意:
给定N,M,L,R,从区间[L,R]内选出N个整数使得它们的gcd恰好为m，求合法的选数方案数对1e9+7取模的值。1<=N,M,L,R<=1e9, R-L<=1e5.

思路分析:
gcd? 那就莫比乌斯反演好了。
令F(m)表示从[L,R]中选出N个数，其gcd为m的倍数的方案数。
f(m)表示从[L,R]中选出N个数，其gcd 恰好为m方案数。(莫比乌斯反演常见做法)
我们要求的是f(m),为了简化运算，我们令l等于大于等于L的最小的m的倍数，r等于小于等于L的最大的m的倍数。然后l/=m,r/=m,问题转化为求f(1). (莫比乌斯反演常见做法)
根据莫比乌斯反演公式$$F(n)=sum_{n|d} f(d), f(n)=sum_{n|d}mu (frac{d}{n})F(d)$$, F(n)可以O(1)求得，直接反演即可。
现在面临两个问题:

1. F(x)和f(x)的定义域是什么?
2. 如何O(1)求F(x)?

先来解决第二个问题:
F(x)其实就是[l,r]内是x的倍数的数的个数的N次方，可以用快速幂求得。具体见代码getF函数。

难点在于第一个问题:
首先我们知道，定义域不超过r. 而r=R/M是1e9级别的，因此必须优化，发现更多的性质。
F(x)既然表示选出N个数gcd为x的方案数，那我们观察以下式子$$gcd (x,y)le y-x (x<y)$$如果选的数不全相等，那它们的gcd一定不会超过r-l, 也就是F(x)和f(x)的定义域就会缩小到r-l, 而r-l是1e5级别的！这就很美妙了！
现在只要处理一下选出的所有数全相等的情况了。
为了缩小定义域，我们给F(x)和f(x)分别添加一个条件: F(x)表示表示从[L,R]中选出不全相等的 N个数，其gcd为x的倍数的方案数，f(x)表示表示从[L,R]中选出不全相等的 N个数，其gcd 恰好为x的方案数，枚举定义域[1,r-l]莫比乌斯反演求出f(1)即可。
而定义变了以后，O(1)计算F(x)的方法也出现了变动: $$F(x)=a^N-a$$其中a为[l,r]内是x的倍数的数的个数。公式解释: 如果是随意选，共有(a^N)种选法，然后去掉全部相等的选法，选N个全部相等的数就相当于只选一个数，因此有a种选法，从(a^N)中扣除。
以上是计算f(1)的方法。
f(1)算完后，还要加上从[l,r]中选N个全相等的数使得gcd为1的方案数。那显然唯一方案就是全选1，如果1被包含在区间[l,r]中答案就是f(1)+1,否则答案为f(1).

代码实现:

#include<cstdio>
using namespace std;

const int N = 1e5+1;
const long long P = 1e9+7;
long long n,m,lb,rb;
int mu[N+4];
long long p[N+4];
bool f[N+4];
int pn;

void Mobius()
{
	mu[1] = 1; pn = 0;
	for(int i=2; i<=N; i++)
	{
		if(!f[i]) {pn++; p[pn] = i; mu[i] = -1;}
		for(int j=1; j<=pn && i*p[j]<=N; j++)
		{
			f[p[j]*i] = true;
			if(i%p[j]==0) {mu[i*p[j]] = 0; break;}
			else mu[i*p[j]] = -mu[i];
		}
	}
}

long long quickpow(long long a,long long b)
{
	a %= P;
	long long cur = a,ret = 1ll;
	for(int i=0; b; i++)
	{
		if(b&(1ll<<i)) {ret *= cur; ret %= P; b-=(1ll<<i);}
		cur *= cur; cur %= P;
	}
	return ret;
}

long long getF(long long a)
{
	long long lt,rt;
	if(lb%a>0ll) lt = lb/a+1;
	else lt = lb/a;
	rt = rb/a;
	return (quickpow(rt-lt+1,n)-(rt-lt+1)+P)%P;
}

int main()
{
	Mobius();
	scanf("%lld%lld%lld%lld",&n,&m,&lb,&rb);
	if(lb%m>0ll) lb = lb/m+1;
	else lb = lb/m;
	rb/=m;
	long long nn = rb-lb,ans = 0ll;
	for(int i=1; i<=nn; i++)
	{
		ans += mu[i]*getF(i);
		ans = (ans+P)%P;
	}
	if(lb<=1 && 1<=rb) {ans++; ans%=P;}
	printf("%lld
",ans);
	return 0;
}