@hdu

目录

• @description@
• @solution@
  • @part - 1@
  • @part - 2@
• @accepted code@
• @details@

---

@description@

将 n 本书放入 k 个书架，设第 i 个书架中有 (cnt_i)，再设 ({fi}) 表示斐波那契数列，其中 (f[0] = 0, f[1] = 1)。

一个放书方案的权值记为 (gcd(2^{f[cnt_1]} - 1, 2^{f[cnt_2]} - 1, ..., 2^{f[cnt_K]} - 1))，求在随机放书的前提下的期望权值。

Input
第一组包含一个数据组数 T (T≤10)。
接下来 T 行，每行两个整数 n,k(0<n,k≤10^6)。

Output
对于每组数据，输出答案模 10^9 + 7。

Sample Input
1
6 8

Sample Output
797202805

@solution@

不难用组合数算出总方案数为 (C_{n+k-1}^{k-1})，于是可以得到答案的表达式为：

[ans = frac{1}{C_{n+k-1}^{k-1}}*(sum_{(sum_{i=1}^{k}cnt_i) = n}gcd(2^{f[cnt_1]} - 1, 2^{f[cnt_2]} - 1, ..., 2^{f[cnt_K]} - 1)) ]

@part - 1@

首先让我们发现几个性质：
（1）通过百度可以发现 (gcd(2^M - 1, 2^N - 1) = 2^{gcd(M, N)} - 1)。
（2）依然是通过百度可以发现 (gcd(f[i], f[j]) = f[gcd(i, j)])。

简要地证明（跟题目本身无关）：
证明（1）：
使用辗转相减法（不妨假设 N > M）：(gcd(2^M - 1, 2^N - 1) = gcd(2^M - 1, 2^N - 2^M) = gcd(2^M - 1, (2^{N-M} - 1)*2^M))。
由于 2^M 与 2^M - 1 互质，所以 (gcd(2^M - 1, 2^N - 1) = gcd(2^M - 1, 2^{N - M} - 1))。可以发现这个等同于在指数上进行辗转相减，最后得到的结果必然是 (2^{gcd(N, M)}-1)。
证明（2）：
首先不难根据斐波那契的递推公式以及辗转相减得到 (gcd(f[i-1], f[i]) = 1)。
然后有 (f[i] = f[i-1] + f[i-2] = 2*f[i-2] + f[i-3] = ... = f[j]*f[i-j+1] + f[j-1]*f[i-j])。
依然假设 i > j，则 (gcd(f[i], f[j]) = gcd(f[j]*f[i-j+1] + f[j-1]*f[i-j], f[j]))。
对其使用辗转相减得到 (gcd(f[i], f[j]) = gcd(f[j-1]*f[i-j], f[j]))。又因 (gcd(f[j-1], f[j]) = 1)，所以有 (gcd(f[i], f[j]) = gcd(f[i-j], f[j]))。又变成了对下标进行辗转相减，因而得证。

于是我们只需要对于每一个 x 求满足 ((sum_{i=1}^{k}cnt_i) = n) 且 (gcd(cnt_1, cnt_2, ..., cnt_k) = x) 的方案数，不妨令其为 v[x]。
则答案的表达式变为：

[ans = frac{1}{C_{n+k-1}^{k-1}}*(sum_{i=1}^{n}(2^{f[i]} - 1)*v[i]) ]

@part - 2@

可以发现 v[x] 长得就非常莫比乌斯反演，于是我们考虑反演，得到：

[v[x] = sum_{x|d}C_{frac{n}{x}+k-1}^{k-1}*mu(d/x) ]

其中要求 x, d 都是 n 的因子。

因为我们相当于是要求 ((sum_{i=1}^{k}cnt_i) = n) 且 (x | gcd(cnt_1, cnt_2, ..., cnt_k)) 的方案数。
将所有 cnt 除以 x 转为求 ((sum_{i=1}^{k}cnt_i) = frac{n}{x}) 的方案数，就是如上的组合数。

因为 x, d 都是 n 的因子，枚举 x, d 是 (O(因子个数^2))，可以发现这个复杂度的一个上界为 (O(sqrt{n}^2) = O(n)) 的。于是直接暴力做。

@accepted code@

#include<cstdio>
const int MOD = int(1E9) + 7;
const int MAXN = 3000000;
int fct[MAXN + 5], ifct[MAXN + 5], f[MAXN + 5];
int pow_mod(int b, int p) {
	int ret = 1;
	while( p ) {
		if( p & 1 ) ret = 1LL*ret*b%MOD;
		b = 1LL*b*b%MOD;
		p >>= 1;
	}
	return ret;
}
int prm[MAXN + 5], miu[MAXN + 5], pcnt = 0;
bool nprm[MAXN + 5];
void sieve() {
	miu[1] = 1;
	for(int i=2;i<=MAXN;i++) {
		if( !nprm[i] )
			prm[++pcnt] = i, miu[i] = -1;
		for(int j=1;1LL*i*prm[j]<=MAXN;j++) {
			nprm[i*prm[j]] = true;
			if( i % prm[j] == 0 ) {
				miu[i*prm[j]] = 0;
				break;
			}
			else miu[i*prm[j]] = -miu[i];
		}
	}
}
void init() {
	f[0] = 1, f[1] = 2;
	for(int i=2;i<=MAXN+2;i++) {
		f[i] = 1LL*f[i-1]*f[i-2]%MOD;
		f[i-2] = (f[i-2] + MOD - 1)%MOD;
	}
	fct[0] = 1;
	for(int i=1;i<=MAXN;i++)
		fct[i] = 1LL*fct[i-1]*i%MOD;
	ifct[MAXN] = pow_mod(fct[MAXN], MOD-2);
	for(int i=MAXN-1;i>=0;i--)
		ifct[i] = 1LL*ifct[i+1]*(i+1)%MOD;
	sieve();
}
int comb(int n, int m) {return 1LL*fct[n]*ifct[m]%MOD*ifct[n-m]%MOD;}
int div[MAXN + 5];
void solve() {
	int n, k; scanf("%d%d", &n, &k);
	int ans = 0, cnt = 0;
	for(int i=1;i<=n;i++)
		if( n % i == 0 ) div[++cnt] = i;
	for(int i=cnt;i>=1;i--)
		for(int j=i;j<=cnt;j++)
			if( div[j] % div[i] == 0 )
				ans = (ans + 1LL*miu[div[j]/div[i]]*comb(n/div[j]+k-1, k-1)%MOD*f[div[i]]%MOD)%MOD;
	printf("%lld
", (1LL*ans*pow_mod(comb(n+k-1, k-1), MOD-2)%MOD + MOD)%MOD);
}
int main() {
	init();
	int T; scanf("%d", &T);
	for(int i=1;i<=T;i++) solve();
}

@details@

记得去年赛场做这道题的时候先是使用的 O(nlog n) 的算法，最后被卡了才被迫换成 O(n) 算法。
但是现在看来。。。好像自己根本 BB 不出来 O(nlogn) 的算法？
或许这就算是人的一种成长吧。