「CQOI2015」选数

题目描述

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

由于方案数较大，你只需要输出其除以 (10^9+7) 的余数即可。

输入格式

输入一行，包含四个空格分开的正整数，依次为 (N,K,L,H)。

输出格式

输出一个整数，为所求方案数除以 (10^9+7) 的余数。

输入输出样例

输入

2 2 2 4

输出

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\%) 的数据，(1le N,Kle 10^9)，(1le Lle Hle 10^9)，(H-Lle 10^5)。

题解

问题可以直接转化为一个式子，然后进行化简

[egin{aligned} sum_{i=L}^Hsum_{j=L}^H[gcd(i,j)=k] &=sum_{i=lfloor frac{L}{k} floor}^{lfloor frac{H}{k} floor}sum_{j=lfloor frac{L}{k} floor}^{lfloor frac{H}{k} floor}[gcd(i,j)=1]\ &=sum_{i=lfloor frac{L}{k} floor}^{lfloor frac{H}{k} floor}sum_{j=lfloor frac{L}{k} floor}^{lfloor frac{H}{k} floor}sum_{d|i,d|j}mu(d)\ &=sum_{i=lfloor frac{L}{k} floor}^{lfloor frac{H}{k} floor}sum_{j=lfloor frac{L}{k} floor}^{lfloor frac{H}{k} floor}sum_{d|i,d|j}mu(d)\ &=sum_{d=1}^{lfloor frac{H}{k} floor}mu(d) sum_{i=lfloor frac{L}{k} floor}^{lfloor frac{H}{k} floor}sum_{j=lfloor frac{L}{k} floor}^{lfloor frac{H}{k} floor}[dmid i][dmid j]\ &=sum_{d=1}^{lfloor frac{H}{k} floor}mu(d) (left lfloor frac{lfloor frac{H}{k} floor}{d} ight floor - left lfloor frac{lfloor frac{L-1}{k} floor}{d} ight floor)^2 end{aligned} ]

前面的 (mu) 杜教筛掉就行了，外层的循环整除分块去做，总复杂度是 (mathcal{O}(n^{frac{2}{3}})) 的

代码

#include <cstdio>
#include <iostream>
#include <unordered_map>

typedef long long ll;

using namespace std;

const int mod = 1e9 + 7;

inline int addmod (register int a, register int b) {
	return a += b, a >= mod ? a - mod : a;
}

inline int delmod (register int a, register int b) {
	return a -= b, a < 0 ? a + mod : a;
}

inline ll mulmod (register ll a, register int b) {
	return a *= b, a >= mod ? a % mod : a;
}

inline int qpow (register int a, register int b, register int ans = 1) {
	for (; b; b >>= 1, a = mulmod (a, a))
		if (b & 1) ans = mulmod (ans, a);
	return ans;
}

int n, m, L, R, cnt;
int prime[3000005], mul[3000005], sum[3000005];
bool vis[3000005];
unordered_map <int, int> f;
	
inline void xxs () {
	mul[1] = sum[1] = 1;
	for (register int i = 2; i <= 3e6; i ++) {
		if (! vis[i]) prime[++ cnt] = i, mul[i] = mod - 1;
		for (register int j = 1; j <= cnt && i * prime[j] <= 3e6; j ++) {
			vis[i * prime[j]] = 1;
			if (i % prime[j] == 0) break;
			mul[i * prime[j]] = mod - mul[i];
		}
		sum[i] = addmod (sum[i - 1], mul[i]);
	}
}

inline int F (register int n, register int ans = 1) {
	if (n <= 3e6) return sum[n];
	if (f[n]) return f[n];
	for (register int l = 2, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = delmod (ans, mulmod (r - l + 1, F (n / l)));
	return f[n] = ans;
}

inline int Calc (register int n, register int m, register int k, register int ans = 0) {
	for (register int l = 1, r; l <= m; l = r + 1) 
		r = n < l ? m / (m / l) : min (n / (n / l), m / (m / l)), 
		ans = addmod (ans, mulmod (delmod (F (r), F (l - 1)), qpow (delmod (m / l, n / l), k)));
	return ans;
}

int main () {
	scanf ("%d%d%d%d", &n, &m, &L, &R), L = (L - 1) / m, R = R / m, xxs ();
	return printf ("%d
", Calc (L, R, n)), 0;
}