【2018 Multi-University Training Contest 6】bookshelf（莫比乌斯反演）

Problem Description
Patrick Star bought a bookshelf, he named it ZYG !!

Patrick Star has N book .

The ZYG has K layers (count from 1 to K) and there is no limit on the capacity of each layer !

Now Patrick want to put all N books on ZYG :

1. Assume that the i-th layer has cnti(0≤cnti≤N) books finally.

2. Assume that f[i] is the i-th fibonacci number (f[0]=0,f[1]=1,f[2]=1,f[i]=f[i−2]+f[i−1]).

3. Define the stable value of i-th layers stablei=f[cnti].

4. Define the beauty value of i-th layers beautyi=2stablei−1.

5. Define the whole beauty value of ZYG score=gcd(beauty1,beauty2,…,beautyk)(Note: gcd(0,x)=x).

Patrick Star wants to know the expected value of score if Patrick choose a distribute method randomly !

Input
The first line contain a integer T (no morn than 10), the following is T test case, for each test case :

Each line contains contains three integer $n,k(0<n,k≤10^6).$

Output
For each test case, output the answer as a value of a rational number modulo $10^9+7.$

Formally, it is guaranteed that under given constraints the probability is always a rational number pq (p and q are integer and coprime, q is positive), such that q is not divisible by 109+7. Output such integer a between 0 and 109+6 that p−aq is divisible by $10^9+7.$

Sample Input
1
6 8

Sample Output
797202805

Source
2018 Multi-University Training Contest 6

前置知识：

$$gcd(2^a-1,2^b-1)=2^{gcd(a,b)}-1$$
$$gcd(F[a],F[b])=F[gcd(a,b)]$$

然后用组合数算分配数量（隔板法）
在莫比乌斯反演算枚举每个gcd的贡献数量
代码：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#define maxx 1000005
#define _maxx 2000005
#define mod 1000000007
#define ll long long
using namespace std;
bool isP[maxx];
int prime[maxx];
int mu[maxx];
ll f[maxx];
ll fac[maxx<<1];
ll inv[maxx<<1];
int cnt;
ll P(ll a,ll b)
{
    ll ans=1;
    while(b)
    {
        if(b&1)ans=ans*a%mod;
        b>>=1;
        a=a*a%mod;
    }
    return ans;
}
void init()
{
    mu[1]=1;
    for(int i=2;i<maxx;i++)
    {
        if(!isP[i])
        {
            prime[cnt++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<cnt&&(ll)i*prime[j]<maxx;j++)
        {
            isP[i*prime[j]]=true;
            if(i%prime[j])
                mu[i*prime[j]]=-mu[i];
            else
            {
                mu[i*prime[j]]=0;
                break;
            }
        }
    }

    f[0]=0;f[1]=1;
    for(int i=2;i<=maxx;i++)f[i]=(f[i-1]+f[i-2])%(mod-1);

    fac[0]=1;
    for(int i=1;i<_maxx;i++)fac[i]=fac[i-1]*i%mod;
    inv[_maxx-1]=P(fac[_maxx-1],mod-2);
    for(int i=_maxx-2;i>=0;i--)inv[i]=inv[i+1]*(i+1)%mod;
}
ll C(int n,int m)
{
    return fac[n]*inv[m]%mod*inv[n-m]%mod;
}
int main()
{
    init();
    int n,k;
    int t;
    cin>>t;
    while(t--)
    {
        scanf("%d%d",&n,&k);
        ll ans=0;
        for(int i=1;i<=n;i++)
        {
            if(n%i)continue;
            ll _count=0;
            for(int j=i;j<=n;j+=i)
                if(n%j==0)_count=(_count+mu[j/i]*C(n/j+k-1,k-1))%mod;
            _count=(_count+mod)%mod;
            ans=(ans+_count*(P(2,f[i])-1)%mod)%mod;
            //cout<<P(2,f[i])<<endl;
        }
        ans=(ans+mod)%mod;
        //cout<<ans<<endl;
        ll total=C(n+k-1,k-1);
        //cout<<total<<endl;
        printf("%lld\n",ans*P(total,mod-2)%mod);
    }
    return 0;
}