HDU-6363 bookshelf丨排列组合丨数论丨容斥丨附两个结论来源

HDU-6363 bookshelf
给出n本书，和k层书架，cnti表示i层书架放的书本数量，求下式的期望：

$score=gcd(2^{f[cnt_{1}]}-1,2^{f[cnt_{2}]}-1,\cdots,2^{f[cnt_{k}]}-1)$
其中有2个结论 【1式来源】 【2式证明】:

$gcd(2^a-1,2^b-1,\cdots)=2^{gcd(a,b,\cdots)}-1$

$gcd(f[cnt1],f[cnt2])=f[gcd(cnt1,cnt2)]$

以上两式分别可以在英文维基的gcd词条与Fibonacci number词条找到，但显然都不是一般人随随便便能推出来的，所以只能靠一定的知识储备了 （是魔鬼吧？）

问题转化为求下式的期望：

$score=2^{f[gcd(cnt_{1},cnt_{2},\cdots,cnt_{k})]}-1$
  到这里问题还只解决到一半。因为gcd(cnti)的结果肯定是n的 某个因子，我们设为g(i)，故枚举因子数求和，再除以方案总数就是答案。
  下面的问题是，我们如何计算n本书，放到k层书架的方案总和？可以利用隔板法的思维，书架是隔板(k-1)，首先隔板与书进行全排列，因为书架是无序的，除去书架的全排列，又因为书本是无序的，再除去书本的全排列。可以得到下面的式子：

$\frac{(n+k-1)!}{n!\cdot (k-1)!}=C_{n+k-1}^{k-1}$

  这也就是所有书放入书架的组合数。下面由于每种g(i)都会存在一定的方案数，这个方案数也可以类似的求得到，我们把n/g(i)本书看作一组，我们按组放入书架，那么同理，书架作为隔板，去重复上述操作，对于每个g(i)，我们会得到下面这样的公式：

$\frac{(n/g(i)+k-1)!}{(n/g(i))!\cdot (k-1)!}=C_{n/g(i)+k-1}^{k-1}$
  机智的读者有可能会感到疑惑，如果有一个方案，完全没有按照组合放入（比如，{1,2,3,7,11,17}），那不就不满足这个算法了吗？其实不然，因为g(i)是n的所有因子，包括1。若n>k，我们按一本为一组放入书架，不就包含了 所有的方案了吗？由此，机智的读者你肯定会发现，这个总和并不简单，显然这样直接求和得到的方案数存在大量的重复！
  现在的情况是我们可以求出每种g(i)情况下的方案数，但是直接求和加起来，会 存在重复，那么我们需要开始除重，也就是容斥原理的应用。首先对于g(i)=8的情况，在g(i)=4与g(i)=2与g(i)=1中重复计数过。如果还不理解可以这样想，按2本书一组插入k个书架，肯定会包括gcd等于8的情况（题目交代了gcd(x,0)=x），数量足够的话16，32···也都会有包括。那么我们的做法可以是，从最大的向下逐项对g(j)减去g(i)的方案数，这样就完成了除重。

完成了上面这些之后，答案的计算公式就是：

$ans=\frac{\sum_{i=1}^{cnt}rav(i)\cdot (2^{f[div(i)]}-1)}{C_{n+k-1}^{k-1}}$

注1 命名同下面代码，div(i)表示第i个因子（也就是g(i)），rav(i)表示div(i)因子所对应的方案数（要经过去重）
注2 考虑到2的指数有可能会非常大，结果又要需要取模，所以计算时我们可以利用欧拉降幂公式来处理（还有一种办法时预处理2^fi，但是学个新东西也不吃亏 。

$A^BmodC=A^{Bmod\phi(C)+\phi(C)}modC$

其中φ是欧拉函数（表示n以内所有与n互质的数的个数）显然φ(1e9+7)=1e9+6，因为1e9+7是质数。具体的欧拉函数的知识需要去别的地方了解。

#include<cstdio>
typedef long long ll;
const int maxn = 2e6+5;
const long long mod = 1e9+7;
ll fib[maxn];
ll inv[maxn],fac[maxn],invf[maxn];
ll rav[maxn];
int div[maxn];

ll fpow(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1)ans=(ans*a)%mod;
        b>>=1;
        a=(a*a)%mod;
    }
    return ans;
}
ll calinv(ll x)
{
    return fpow(x,mod-2);
}
void Init()
{
    fib[0]=0,fib[1]=fib[2]=1;
    for(int i=3;i<maxn;i++)//fib数列
        fib[i]=(fib[i-1]+fib[i-2])%(mod-1);
        //欧拉降幂，其中mod-1是φ(mod)的值

    inv[1]=1;//逆元
    for(int i=2;i<maxn;i++)
        inv[i]=((mod-mod/i)*inv[mod%i])%mod;

    fac[0]=1;//累乘
    for(int i=1;i<maxn;i++)
        fac[i]=(fac[i-1]*i)%mod;

    invf[0]=1;//累乘逆元
    for(int i=1;i<maxn;i++)
        invf[i]=(invf[i-1]*inv[i])%mod;
}

ll C(ll a,ll b)
{
    return fac[a]*invf[b]%mod*invf[a-b]%mod;
}
int main()
{
    int Kase,n,k,bok,cnt;
    ll sum,tot;
    Init();
    scanf("%d",&Kase);
    while(Kase--)
    {
        scanf("%d %d",&n,&k);
        bok=n,cnt=0;

        for(int i=1;i<=bok;i++)
            if(bok%i==0)
                div[cnt++]=i;

        for(int i=0;i<cnt;i++)//公式
            rav[i]=C(n/div[i]+k-1,k-1);

        for(int i=cnt-2;i>=0;i--)//去重
        for(int j=i+1;j<cnt;j++)
            if(div[j]%div[i]==0)
                rav[i]=(rav[i]-rav[j]+mod)%mod;

        sum=0;
        for(int i=0;i<cnt;i++)//公式，求和
            sum=(sum+(fpow(2,fib[div[i]])+mod-1)%mod*rav[i]%mod)%mod;
                            //其中mod-1又是φ(mod)

        tot=calinv(C(n+k-1,k-1));
        //C(n+k-1,k-1)方案总数
        printf("%lld\n",(sum*tot)%mod);
    }
    return 0;
}