CQOI2015 选数

标签： 数学 莫比乌斯反演 递推 容斥原理

---

题目描述

我们知道，从区间[L,H]（L和H为整数）中选取N个整数，总共有$(H-L+1)^N$种方案。小z很好奇这样选出的数的最大公约数的规律，他决定对每种方案选出的N个整数都求一次最大公约数，以便进一步研究。然而他很快发现工作量太大了，于是向你寻求帮助。你的任务很简单，小z会告诉你一个整数K，你需要回答他最大公约数刚好为K的选取方案有多少个。由于方案数较大，你只需要输出其除以1000000007的余数即可。

输入输出格式

输入格式：
输入一行，包含4个空格分开的正整数，依次为N，K，L和H。

输出格式：
输出一个整数，为所求方案数。

输入输出样例

输入样例#1：
2 2 2 4

输出样例#1：
3

说明

样例解释:
所有可能的选择方案：$(2, 2), (2, 3), (2, 4), (3, 2), (3, 3), (3, 4), (4, 2), (4, 3), (4, 4)$
其中最大公约数等于2的只有3组：$(2, 2), (2, 4), (4, 2)$

对于100%的数据，1<=N,K<=10^9，1<=L<=H<=10^9，H-L<=10^5

题解

方法一（莫比乌斯反演）

题意一句话描述:
$$\sum_{a1=L}^{H}\sum_{a2=L}^{H}…\sum_{an=L}^{H}[gcd(a1,a2,…,an)=k]$$
看到这个格式就比较熟悉了，我们不妨直接设:
$$f(x)=\sum_{a1=L}^{H}\sum_{a2=L}^{H}…\sum_{an=L}^{H}[gcd(a1,a2,…,an)=x]$$
$$F(x)=\sum_{x|d}f(d)=(\frac{H}{x}-\frac{L-1}{x})^{n}$$
（上界是$\frac{H}{x}$，下界是$\frac{L-1}{x}+1$）
那么
$$Ans=f(k)=\sum_{k|d}u(\frac{d}{k})F(d)$$
考虑枚举k的倍数:
$$Ans=\sum_{d=1}^{\frac{H}{k}}u(d)(\frac{H}{dk}-\frac{L-1}{dk})^{n}$$
做到这里就可以整除分块了，不过考虑到H达到$10^9$因此我们需要用杜教筛

方法二（容斥+递推）

那么注意到这道题有一个特殊的条件:$H-L<=10^5$
因此我们可以考虑到从这方面入手
一个简单的处理是先缩小L和H的范围:
令$l=\frac{L-1}{k}+1$, $r=\frac{R}{k}$那么问题转换成求$[l,r]$内$gcd(a1,a2,…,an)=1$的数对
然后这里存在一条性质:
在$[l,r]$中任意取不相同的两个数,这两个数的$gcd$不会超过$R-L$
证明：
设x,y两个数的gcd为d，则

$$x-y=d\ast p1-d\ast p2=d\ast (p1-p2)>d$$

证明完毕
有了上述的铺垫，我们可以想到设$f(x)$表示在$[l,r]$中选n个数且这n个数不完全相同且公约数为x的倍数的对数，那么
$$f(x)=d(x)^N-d(x)$$
其中$d(x)$表示区间内x的倍数个数
但是我们怎样求出公约数刚好为x的对数呢，我们发现，$f(x)$是包含了$f(2x),f(3x),….,f(kx)$的个数的，因此我们只要减去x的倍数的$f$就可以了,此时$f(x)$的含义就变成了gcd为x的个数了，因此我们可以倒推，最后答案为$f(1)$
注意，如果$l=1$的话最后我们答案要+1，因为存在$(1,1,1,1,…,1)$这个数对

代码(方法一)

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<tr1/unordered_map>
#define ll long long
#define N 2000100
using namespace std;

const ll p = 1000000007;
int n, k, L, H, cnt = 0;
int pri[N], u[N], sum[N];
bool vis[N];
tr1::unordered_map<int, int> tu;

void get()
{
    u[1] = 1; u[0] = 0;
    int t = N - 100;
    for(int i = 2; i <= N; i ++)
    {
        if(!vis[i]) pri[++ cnt] = i, u[i] = -1;
        for(int j = 1; j <= cnt && i * pri[j] <= t; j ++)
        {
            vis[i * pri[j]] = true;
            if(i % pri[j]) u[i * pri[j]] = -u[i];
            else break; 
        }
    }
    for(int i = 1; i <= t; i ++) sum[i] = sum[i - 1] + u[i];
}

int djsu(int x)
{
    if(x <= 2e6) return sum[x];
    if(tu[x]) return tu[x];
    int ans = 1;
    for(int l = 2, r; l <= x; l = r + 1)
    {
        r = x / (x / l);
        ans -= (r - l + 1) * djsu(x / l);
    }
    return tu[x] = ans;
}

ll pow_mul(ll a, ll b, ll p)
{
    ll ans = 1;
    while(b)
    {
        if(b & 1) ans = ans * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ans;
}

int main()
{
    scanf("%d%d%d%d", &n, &k, &L, &H);
    get();
    ll ans = 0;
    int le = (L - 1) / k, re = H / k;
    for(int l = 1, r; l <= re; l = r + 1) 
    {
        if(le / l == 0) r = re / (re / l); //除数不能为0 
        else r = min(le / (le / l), re / (re / l));
        ans += 1ll * (djsu(r) - djsu(l - 1)) * pow_mul(re / l - le / l, n, p);
        ans %= p;
    }
    while(ans < 0) ans += p; 
    printf("%lld", ans);
    return 0;
} 
/*
3 5 214147541 321411411
*/

心得

莫比乌斯函数本质是一个容斥函数，下界是要从1开始的，否则就不满足莫比乌斯函数的某些性质了，感谢peng-ym大佬的讲解