bookshelf HDU - 6363(数论结论+莫比乌斯反演)

bookshelf HDU - 6363

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≤106)$
.
Output
For each test case, output the answer as a value of a rational number modulo 109+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 109+7
.
Sample Input

1
6 8

Sample Output

797202805

题意：

把N本书放到K层的书架上，每一层的美丽值为$b_{i}=2^{f[cnt]}-1$，其中cnt是这一层书的数量，f[x]为斐波那契数列，整个书架的美丽值为$gcd(b_{1},b_{2},...,b_{k})$，问整个书架的美丽值的期望

分析：

首先需要利用数论中的两个结论进行化简

$gcd(x^{a}-1,x^{b}-1)=x^{gcd(a,b)}-1$

$gcd(f[x],f[y])=f[gcd(x,y)]$($f[]$是斐波那契数列)

于是可以进一步化简公式为$2^{f[gcd(cnt1,cnt2,...,cntk)]}-1$

这样只需要枚举gcd，计算gcd出现次数就可以知道gcd的贡献值了

设$f(g)表示gcd(x_1,x_2,...x_k) =g的个数$

$F(p)表示gcd(x_1,x_2,...x_k) = g的倍数p的个数$

我们发现对于$f(g)$来说不是很好求，但是$F(p)$好求

如果想知道n内有多少个g的倍数的话那么应该共有$\frac{n}{p} 个$

那么将这些数分成k份的话，利用隔板法

$|*|*|...|**|**|$

$k+1$个板子可以得到k个盒子，中间盒子中任意放$\frac{n}{p} 个$数

那么相当于中间$\frac{n}{p}+k-1$个位置的组合，其中选取k-1个位置，那么其他位置便确定了

因此$F(p) = C_{\frac{n}{p}+k-1}^{k-1}$

而gcd的贡献值应该为$f(gcd)\times (2^{f[gcd]}-1)$

这样我们知道了$F()$就可以通过莫比乌斯反演计算$f()$

因此可以先在$n的范围内枚举gcd（n | gcd），然后在从gcd为基础枚举gcd的倍数p(p|n)$

$f(gcd) = \sum_{gcd|p}\mu (\frac{p}{gcd})F(p) = \sum_{gcd|p}\mu(\frac{p}{gcd})C_{\frac{n}{p}+k-1}^{k-1}$

这样求出$f()$后可以乘上值$(2^{f[gcd]}-1)$

所以答案为

$ans = \sum_{gcd|n}2^{f[gcd]}-1 \sum_{gcd|p}\mu(\frac{p}{gcd})C_{\frac{n}{p}+k-1}^{k-1}$

因为要求期望，所以再除去总的可能数

$ans = \frac{ans}{C_{n+k-1}^{k-1}}$

code:

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int maxn = 1e6 + 5;
const LL mod = 1e9 + 7;
const LL mm = 1e9 + 6;

LL inv[maxn * 2 + 5],fac[maxn * 2 + 5];
bool check[maxn];
LL prime[maxn],mu[maxn],f[maxn],n,k;

LL exgcd(LL a,LL b,LL& x,LL& y){
    if(a == 0 && b == 0) return -1;
    if(b == 0){
        x = 1,y = 0;
        return a;
    }
    LL d = exgcd(b,a % b,y,x);
    y -= a / b * x;
    return d;
}
LL rev(LL a,LL n){
    LL x,y;
    LL d = exgcd(a,n,x,y);
    return (x % n + n) % n;
}
void init_inv(){
    fac[0] = 1;
    for(LL i = 1;i < 2 * maxn;++i) fac[i] = fac[i - 1] * i % mod;
    inv[2 * maxn - 1] = rev(fac[2 * maxn - 1],mod);
    for(LL i = 2 * maxn - 2;i >= 0;--i){
        inv[i] = (i + 1) * inv[i + 1] % mod;
    }
}
LL C(LL a,LL b){
    LL m = a,n = b;
    return ((fac[n] * inv[m]) % mod * inv[n - m]) % mod;
}
void init(){
    f[0] = 0,f[1] = 1;
    for(int i = 2;i < maxn;++i){
        f[i] = (f[i - 1] + f[i - 2]) % mm;
    }
}
void mobius(){
    memset(check,false,sizeof(check));
    mu[1] = 1;
    int tot = 0;
    for(LL i = 2;i < maxn;++i){
        if(!check[i]){
            prime[tot++] = i;
            mu[i] = -1;
        }
        for(int j = 0;j < tot;++j){
            if(i * prime[j] > maxn) break;
            check[i * prime[j]] = true;
            if(i % prime[j] == 0){
                mu[i * prime[j]] = 0;
                break;
            }
            else mu[i * prime[j]] = -mu[i];
        }
    }
}
LL qmod(LL a,LL n){
    LL ret = 1;
    while(n){
        if(n & 1) ret = ret * a % mod;
        a = a * a % mod;
        n = n>>1;
    }
    return ret;
}

int main(){
    init_inv();
    mobius();
    init();
    int T;
    scanf("%d",&T);
    while(T--){
        scanf("%lld%lld",&n,&k);
        LL denominator = C(k-1,n + k - 1);
        denominator = rev(denominator,mod);
        LL ans = 0;
        for(int i = 1;i <= n;++i){
            if(n % i) continue;
            LL div = i;
            LL contri = (qmod(2,f[div]) - 1 + mod) % mod;
            LL cnt = 0;
            for(LL j = div;j <= n;j += div){
                if(n % j) continue;
                LL d = j;
                cnt += mu[d / i] * C(k-1,n / d + k - 1) % mod;
                cnt = (cnt % mod + mod) % mod;
            }
            ans += contri * cnt % mod;
            ans %= mod;
        }
        ans = ans * denominator % mod;
        printf("%lld\n",ans);
    }
    return 0;
}