题目地址：洛谷：【P3172】[CQOI2015]选数 – 洛谷、BZOJ：Problem 3930. — [CQOI2015]选数

题目描述

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

题意简述

求区间[L, H]中选择N个整数（可重），求选出数组成的可重集gcd为K的方案数%1e9+7的值。

输入输出格式

输入格式：
输入一行，包含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

题解

参考资料：题解 P3172 【[CQOI2015]选数】 – xyz32768 的博客 – 洛谷博客
我们观察到L与H虽然很大，但是H-L并不大，这是一个我们能利用的地方。首先，可以先行缩小区间范围，即令 L' = \lceil \frac{L}{K} \rceil, H' = \lfloor \frac{H}{K} \rfloor ，问题转变为[L’, H’]区间内选出N个gcd为1的数字的方案数了。
我们考虑一个DP+去重的思想，令f[i]表示含公约数为i且选出的数不全相同的方案数，令[L’, H’]中i的倍数的个数为x，则f[i] = x^n - x。而我们想求的是g[i]表示最大公约数为i且选出的数不全相同的方案数中的g[1]，我们容易知道g[i] = f[i] - g[2i] - g[3i] - \cdots，因此可以通过反推求出g[i]，从而推得g[1]。f定义中，要求不全相同的原因是如果含相同数字，是因为两两互质肯定不存在选出的数都相同的情况。需要注意的是，如果此L’=1，还可以全选1，因此答案要+1。

代码

// Code by KSkun, 2018/6
#include <cstdio>
#include <cctype>  

#include <algorithm>

typedef long long LL;

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    register char c = fgc();
    while(!isdigit(c)) {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(isdigit(c)) {
        res = (res << 1) + (res << 3) + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 100005, MO = 1e9 + 7;

inline LL fpow(LL n, LL k) {
    LL res = 1;
    while(k) {
        if(k & 1) res = res * n % MO;
        n = n * n % MO;
        k >>= 1;
    }
    return res;
}

LL n, k, l, h, dp[MAXN];

int main() {
    n = readint(); k = readint(); l = readint(); h = readint();
    if(l % k) l = l / k + 1; else l /= k; h /= k;
    if(l > h) {
        puts("0"); return 0;
    }
    for(int i = 1; i <= h - l; i++) {
        LL ll = l, hh = h;
        if(ll % i) ll = ll / i + 1; else ll /= i; hh /= i;
        if(ll > hh) continue;
        dp[i] = ((fpow(hh - ll + 1, n) - (hh - ll + 1)) % MO + MO) % MO;
    }
    for(int i = h - l; i; i--) {
        for(int j = i << 1; j <= h - l; j += i) {
            dp[i] = ((dp[i] - dp[j]) % MO + MO) % MO;
        }
    }
    if(l == 1) dp[1] = (dp[1] + 1) % MO;
    printf("%lld", dp[1]);
    return 0;
}