求在区间[L,H][L,H]中任意取nn个数，它们的最大公约数恰好为kk的方案数mod 109+7\mathrm{mod}~10^9+7的值

n,k≤109,L≤H≤109,H−L≤105n,k\le 10^9, L\le H\le10^9, H-L\le 10^5

Luogu P3172

Solution

设f(k)f(k)表示区间内选nn个恰好gcd\gcd为kk的数的方案数，F(k)F(k)表示区间内选nn个k∣gcdk|\gcd（即gcd\gcd为kk的倍数）的数的方案数

由题意很容易得到F(k)=(⌊Hk⌋−⌊L−1k⌋)nF(k)=(\lfloor \frac{H}{k}\rfloor - \lfloor\frac{L-1}{k}\rfloor)^n

又因为

F(k)=∑k∣df(d) F(k)=\sum_{k|d}f(d)

由莫比乌斯反演可知

f(k)=∑k∣dF(d)μ(dk)=∑i=1⌊Hk⌋F(ik)μ(i)=∑i=1⌊Hk⌋μ(i)(⌊Hik⌋−⌊L−1ik⌋)n \begin{aligned} f(k)=&\sum_{k|d}F(d)\mu(\frac{d}{k})\\ =&\sum_{i=1}^{\lfloor\frac{H}{k}\rfloor}F(ik)\mu(i)\\ =&\sum_{i=1}^{\lfloor\frac{H}{k}\rfloor}\mu(i)(\lfloor \frac{H}{ik}\rfloor - \lfloor\frac{L-1}{ik}\rfloor)^n \end{aligned}

换元，令H=⌊Hk⌋,L=⌊L−1k⌋H=\lfloor\frac{H}{k}\rfloor,L=\lfloor\frac{L-1}{k}\rfloor，则有

ans=∑i=1Hμ(i)(⌊Hi⌋−⌊Li⌋)n ans=\sum_{i=1}^{H}\mu(i)(\lfloor\frac{H}{i}\rfloor-\lfloor\frac{L}{i}\rfloor)^n

很显然，直接杜教筛+整除分块即可

Summary

这个第一步容斥反演的套路很常见，见到类似的这种问题时要往容斥反演上想

Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>

#define x first
#define y second
#define y1 Y1
#define y2 Y2
#define mp make_pair
#define pb push_back

using namespace std;

typedef long long LL;
typedef pair<int, int> pii;

template <typename T> inline int Chkmax (T &a, T b) { return a < b ? a = b, 1 : 0; }
template <typename T> inline int Chkmin (T &a, T b) { return a > b ? a = b, 1 : 0; }

inline void proc_status()
{
	ifstream t ("/proc/self/status");
	cerr << string (istreambuf_iterator <char> (t), istreambuf_iterator <char> ()) <<endl;
}
template <typename T> T read ()
{
	T sum = 0, fl = 1; char ch = getchar();
	for (; !isdigit(ch); ch = getchar()) if (ch == '-') fl = -1;
	for (; isdigit(ch); ch = getchar()) sum = (sum << 3) + (sum << 1) + ch - '0';
	return sum * fl;
}

const int Maxn = 5e6 + 100, Mod = 1000000007;

int N, K, L, H;
int Prime[Maxn], Not_Prime[Maxn], prime_cnt;
LL Mu[Maxn];
__gnu_pbds :: gp_hash_table <int, LL> FMu;

inline void Add (LL &a, int b) { a += b; if (a >= Mod) a -= Mod; }

inline int Pow (int a, int b)
{
	int ans = 1;
	while (b)
	{
		if (b & 1) ans = 1ll * ans * a % Mod;
		a = 1ll * a * a % Mod;
		b >>= 1;
	}
	return ans;
}

inline LL Calc_Mu (int n)
{
	if (n <= Maxn - 5) return Mu[n];
	if (FMu[n]) return FMu[n];
	LL ans = 1;
	for (int l = 2, r; l <= n; l = r + 1)
	{
		r = n / (n / l);
		Add (ans, Mod - 1ll * (r - l + 1) * Calc_Mu(n / l) % Mod);
	}
	return FMu[n] = ans;
}

inline void Solve ()
{
	H /= K, L = (L - 1) / K;
	LL ans = 0;
	for (int l = 1, r; l <= H; l = r + 1)
	{
/**/	if (L && L / l) r = min(H / (H / l), L / (L / l));
		else r = H / (H / l);
		Add (ans, 1ll * (Calc_Mu(r) - Calc_Mu(l - 1) + Mod) % Mod * Pow((H / l - L / l), N) % Mod);
	}
	printf("%lld\n", ans);
}

inline void Input ()
{
	N = read<int>(), K = read<int>(), L = read<int>(), H = read<int>();
}

inline void Init ()
{
	Mu[1] = Not_Prime[1] = 1;
	for (int i = 2; i <= Maxn - 5; ++i)
	{
		if (!Not_Prime[i]) Prime[++prime_cnt] = i, Mu[i] = Mod - 1;
		for (int j = 1; j <= prime_cnt && i * Prime[j] <= Maxn - 5; ++j)
		{
			Not_Prime[i * Prime[j]] = 1;
			if (!(i % Prime[j])) { Mu[i * Prime[j]] = 0; break; }
			Mu[i * Prime[j]] = Mod - Mu[i];
		}
	}

	for (int i = 1; i <= Maxn - 5; ++i) Add (Mu[i], Mu[i - 1]);
}

int main()
{
#ifndef ONLINE_JUDGE
	freopen("A.in", "r", stdin);
	freopen("A.out", "w", stdout);
#endif
	Init();
	Input();
	Solve();
	return 0;
}

Debug

  • 73L: 要特判一下L=0,⌊Ll⌋=0L=0,\lfloor \frac{L}{l}\rfloor=0的情况