@hdu - 6363@ bookshelf
@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}\)，于是可以得到答案的表达式为：
@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]。
则答案的表达式变为：
@part - 2@
可以发现 v[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\n", (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) 的算法？
或许这就算是人的一种成长吧。
