HDU-6363 bookshelf（莫比乌斯反演）

bookshelf

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others)
Total Submission(s): 208    Accepted Submission(s): 96
 

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≤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层的书架上，每一层的美丽值为，其中cnt是这一层书的数量，f[x]为斐波那契数列，整个书架的美丽值为，问整个书架的美丽值的期望

题解：数学选手不在的第N天，，，莫比乌斯都要自己敲了。。。结果还敲残，改成容斥才过。。。

首先有个规律：

这样就可以把化简为

打表还可以知道一个规律：

于是可以进一步化简公式为

这样就只要枚举gcd，计算gcd出现的次数，就能知道gcd的贡献了，这个显然可以用莫比乌斯反演来做。

首先有隔板法：将N分解成k个非负整数的和有种方法。并且考虑到k个最大公因数的为gcd的数之和，一定能整除gcd。

那么，将N分解成k个gcd的倍数的和有种方法

此时gcd的贡献为

当然这里有重复情况，只要枚举gcd的倍数，用莫比乌斯系数容斥一下即可：

#include <bits/stdc++.h>
#define lson l,m,rt<<1
#define rson m+1,r,rt<<1|1
#define x first
#define y second
#define rep(i,a,b) for(int i=a;i<b;++i)
#define per(i,a,b) for(int i=a-1;i>=b;--i)
#define fuck(x) cout<<'['<<#x<<' '<<(x)<<']'
#define FIN freopen("in.txt","r",stdin);
using namespace std;
typedef long long ll;
typedef vector<int> VI;
typedef pair<int, int> PII;

const int INF = 0x3f3f3f3f;
const ll INFLL = 0x3f3f3f3f3f3f3f3fll;
const int mod = 1e9 + 7;
const int MX = 2e6 + 5;

ll F[MX], invF[MX];
ll power(ll a, ll b) {
    ll ret = 1;
    while(b) {
        if(b & 1) ret = (ret * a) % mod;
        a = (a * a) % mod;
        b >>= 1;
    }
    return ret;
}
void initc() {
    F[0] = 1;
    for(int i = 1; i < MX; i++) {
        F[i] = (F[i - 1] * i) % mod;
    }
    invF[MX - 1] = power(F[MX - 1], mod - 2);
    for(int i = MX - 2; i >= 0; i--) {
        invF[i] = invF[i + 1] * (i + 1) % mod; //invF[i]*i!=1,invF[i+1]*i!*(i+1)=1
    }
}

ll C(int n, int m) {
    if(n < 0 || m < 0 || m > n) return 0;
    if(m == 0 || m == n)    return 1;
    return F[n] * invF[n - m] % mod * invF[m] % mod;
}
int prime[MX], mob[MX], p[MX];
void init() {
    int pnum = 0;
    memset (prime, true, sizeof (prime) );
    mob[1] = 1;
    for (int i = 2; i < MX; i++) {
        if (prime[i]) {
            p[pnum ++] = i;
            mob[i] = -1;
        }
        for (int j = 0; j < pnum && (ll) i * p[j] < MX; j++) {
            prime[i * p[j]] = false;
            if (i % p[j] == 0) {
                mob[i * p[j]] = 0;
                break;
            }
            mob[i * p[j]] = -mob[i];
        }
    }
}
ll fob[MX];
int main() {
    fob[0] = 0; fob[1] = 1;
    rep(i, 2, MX) fob[i] = (fob[i - 1] + fob[i - 2]) % (mod - 1);
    init(); initc();
    int T;
    //FIN;
    cin >> T;
    while(T--) {
        int n, k; scanf("%d%d", &n, &k);
        ll ans = 0;
        rep(i, 1, n + 1) {
            if(n % i) continue;
            ll t = 0;
            for(int j = i; j <= n; j += i)
                if(n % j == 0) t = (t + mob[j / i] * C(n / j + k - 1, k - 1) + mod) % mod;
            ans = (ans + t * (power(2ll, fob[i]) - 1 + mod)) % mod;
        }
        ll dx = C(n + k - 1, k - 1);
        dx = power(dx, (ll)mod - 2);
        printf("%lld\n", ans * dx % mod);
    }
    return 0;
}