bookshelf：
又又又刷新了我对数学的感觉（痛并快乐着
首先要对题意进行转化：
有gcd(2 ^ x - 1,2 ^ y - 1) = min(2 ^ x - 1,2 ^ y - 1) = 2 ^ gcd(x,y) - 1。
这里证明也比较简单 ：我们对2 ^ x - 1,2 ^ y - 1二进制展开，那么他们的gcd就是公共的那段1。也就是min(2 ^ x - 1,2 ^ y - 1)。也就是2 ^ gcd(x,y) - 1就是公共长度.
然后又因为gcd(fbi[x],fbi[y] = f[gcd(x,y)]。斐波那契数列的性质，证明可以百度。
那么对题目式子进行变形。
有gcd(beau1,beau2,beau3 ....) = gcd(2 ^ fbi[cnt1] - 1,2 ^ fbi[cnt2] - 1,....) = 2 ^ (gcd(fbi[cnt1],fbi[cnt2]...) - 1 = 2 ^ fbi[gcd(cnt1,cnt2...) - 1
首先扩充一下，N个数，放入K个盒子的方案数：
隔板法 = C(n - 1,k - 1).n 个 1放k - 1个板子。
这里是不能为空的方案数。
如果能为空，那么就默认对每个盒子预先都放了1个。
那就变成了,n + k个数，放入K个盒子的方案数：= C(n + k - 1,k - 1)
那么对于这里的总方案数就是C(n + k - 1,k - 1)。
从上面变形后的式子可以看出，很显然就是有关每个gcd的作为最小值的贡献。
我一开始想的对每个gcd做最小值来讨论，但是这样复杂度N ^ 2了。
这里考虑直接枚举gcd作为最小的那堆，那么这样每堆就要满足是gcd的倍数。
那么我们可以把n堆分成n / gcd堆(因为整除了)。
然后把每堆看成整体隔板法：对n / gcd个堆，放入K个盒子即C(n / gcd + k - 1,k - 1)的方案数。
但是这里可以发现，可能存在最小堆的大小是gcd的倍数，不是gcd的情况，这样的话，在我们后面枚举gcd的倍数作为gcd的时候，就会重复。
所以要容斥，这里的容斥，用到了莫比乌斯系数，因为mu[i] = 1，因子偶数个，mu[i] = -1，因子奇数个。
满足容斥的思想。
gcd的贡献就是$ans = sum_{gcd | n}^{}sum_{p | gcd}^{}mu (p / gcd) * 2^{f[gcd]} - 1 * C(n / p + k - 1,k - 1)$
预处理逆元加速。
fbi数组写成了f数组害我调了好久。。
// Author: levil #include<bits/stdc++.h> using namespace std; typedef long long LL; typedef pair<int,int> pii; const int N = 2e6 + 5; const int M = 1e6 + 5; const LL Mod = 1e9 + 7; #define pi acos(-1) #define INF 1e18 #define dbg(ax) cout << "now this num is " << ax << endl; #define reads(n) FastIO::read(n) namespace FastIO { const int SIZE = 1 << 16; char buf[SIZE], obuf[SIZE], str[60]; int bi = SIZE, bn = SIZE, opt; int read(char *s) { while (bn) { for (; bi < bn && buf[bi] <= ' '; bi++); if (bi < bn) break; bn = fread(buf, 1, SIZE, stdin); bi = 0; } int sn = 0; while (bn) { for (; bi < bn && buf[bi] > ' '; bi++) s[sn++] = buf[bi]; if (bi < bn) break; bn = fread(buf, 1, SIZE, stdin); bi = 0; } s[sn] = 0; return sn; } bool read(int& x) { int n = read(str), bf; if (!n) return 0; int i = 0; if (str[i] == '-') bf = -1, i++; else bf = 1; for (x = 0; i < n; i++) x = x * 10 + str[i] - '0'; if (bf < 0) x = -x; return 1; } }; LL f[N << 1],inv[N << 1],fbi[M],mu[M]; int prime[M],tot = 0; bool vis[N]; LL quick_mi(LL a,LL b) { LL re = 1; while(b) { if(b & 1) re = re * a % Mod; a = a * a % Mod; b >>= 1; } return re; } void init() { f[0] = 1; for(int i = 1;i < N;++i) f[i] = f[i - 1] * i % Mod; inv[N - 1] = quick_mi(f[N - 1],Mod - 2); for(int i = N - 2;i >= 0;--i) inv[i] = inv[i + 1] * (i + 1) % Mod; fbi[0] = 0,fbi[1] = 1; for(int i = 2;i < M;++i) fbi[i] = (fbi[i - 1] + fbi[i - 2]) % (Mod - 1); mu[1] = 1; for(int i = 2;i < M;++i) { if(!vis[i]) { prime[++tot] = i; mu[i] = -1; } for(int j = 1;j <= tot && prime[j] * i < M;++j) { vis[i * prime[j]] = 1; if(i % prime[j] == 0) break; else mu[i * prime[j]] = -mu[i]; } } } LL C(LL n,LL m) { return f[n] * inv[m] % Mod * inv[n - m] % Mod; } int main() { init(); int ca;reads(ca); while(ca--) { int n,k;reads(n),reads(k); LL sum = C(n + k - 1,k - 1); LL ans = 0; for(int i = 1;i <= n;++i) {//gcd if(n % i) continue; LL ff = (quick_mi(2,fbi[i]) - 1 + Mod) % Mod; for(int j = i;j <= n;j += i) {//倍数 if(n % j != 0) continue; LL tmp = mu[j / i] * ff % Mod * C(n / j + k - 1,k - 1) % Mod; ans = ((ans + tmp) % Mod + Mod) % Mod; } } ans = ans * quick_mi(sum,Mod - 2) % Mod; printf("%lld ",ans); } // system("pause"); return 0; }