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\n",ans);
return 0;
}
上一篇:SQL Server 索引和表体系结构(包含列索引)


下一篇:Android笔记---常用控件以及用法