《2018 Multi-University Training Contest 6》

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;
}