HDU-6363 bookshelf 容斥原理

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\le10^{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 $10^{9}+6.$ that p−aq is divisible by $10^{9}+7.$

Sample Input

1 6 8

Sample Output

797202805

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

题解:考试的时候因为本人实在是太菜了。对这个题目完全没有一点思路,考完试以后听了dls直播算是对这个题目有了一点思路代码很短,但是从学习各种姿势到写完用了一天,感觉一下子掌握了好多神奇的姿势,开心。

经典定理: $$\large gcd(x^{a}-1,x^{b}-1)=x^{gcd(a,b)}-1$$

$$\large gcd(fib[a],fib[b])=fib[gcd(a,b)]$$
那么对于这个书架的美丽值可以做如下变形:
QQ截图20181001102636.png
于是我们设$$\large g=gcd(x_{1},x_{2},……,x_{k})$$
那么整个书架的美丽值变为$$\large 2^{fib[g]}……①$$
又因为
QQ截图20181001102720.png
显然g为N的约数而根据①式可得书架美丽值得期望仅与g有关
于是我们枚举N的每一个约数g并统计约数为g的方案数,可以由简单的期望公式就算出答案。
对于g的方案数我们可以对上述②继续变形
QQ截图20181001102739.png

所以方案数就等于③式非负整数解的组数。根据组合数知识我们可以很容易得到③式非负整数解的组数为
$$\Large C^{K-1}_{\frac{N}{g}+K-1}$$
需要注意的是若x*g(x为大于1的整数)也为N的约数那么③式存在这样的解
QQ截图20181001102751.png

所以N的x*g这个因子被重复计算了多次。我们可以简单的用莫比乌斯反演来做计算总方案数
因为我不会莫比乌斯反演所以我用的相对简单的容斥原理.

设约数g的方案数为$$\large f[g]$$
注意到我们在计算g的方案数的时候,将所有g的倍数且为N的约数也计算了。那么我们只要按照从大到小计算N的约数g的方案数,同时减去g的倍数且为N的约数的方案数。用式子表达为:
$$\large f[g] - \sum_{x|N\&g|x}f[x](x != g)$$
上式得到的就是严格是g的方案数

这样我们就得到了美丽值为$2^{fib[g]}$的方案数为$ f[i]$
理论上我们只要将两者相乘,然后除于总共方案数就可以了。但是$fib[g]$是一个指数级增长的函数式,所以我们需要利用欧拉函数对这个式子进行优化
欧拉函数
$$若(a,m)=1则a^{\phi(m)} \equiv1(m) $$
所以在本题中$m=10^{9}+7$
$$2^{m-1} \equiv1(mod\ m) $$
那么$$\large2^{fib[i]} \equiv2^{fib[i]\%(m-1)} (mod\ m) $$
本题中有关数学变化我已经全部做了较为详细的解释。我们接下来只要按照上面的思路先预处理出所有要用的值,然后分别算出来方案数和该方案对应的美丽值,然后除以总方案数就可以了

#include <bits/stdc++.h>

const long long MAXN = 1000000,MOD = 1000000007; 

using namespace std;
long long T,ans,N,K;
long long fib[MAXN+5],fnv[2*MAXN+5],fac[2*MAXN+5],inv[MAXN],f[MAXN];

vector<long long>d;

inline long long powmod(long long a,long long b){
    long long res=1;
    a%=MOD;
    for(;b;b>>=1){
    if(b&1)res=res*a%MOD;
    a=a*a%MOD;
    }
    return res;
}

inline void pre(){
    fib[0] = 0;fib[1] = 1;inv[1] = 1;
    fac[0] = fnv[0] = 1;
    for(int i = 2;i <= MAXN;i++) fib[i] = (fib[i-1] + fib[i-2]) % (MOD - 1);
    for(int i = 2;i <= MAXN;i++) inv[i] = (MOD-MOD/i)*inv[MOD%i]%MOD;
    for(int i = 1;i<=2*MAXN;i++) {
        fac[i] = fac[i-1] * i % MOD;
        fnv[i] = fnv[i-1] * inv[i] % MOD;        
    }
}

inline long long comb(long long a,long long b){
    return ((fac[a]*fnv[b])%MOD*fnv[a-b])%MOD;
}

int main()
{
    pre();
    scanf("%lld",&T);
    while(T--){
        scanf("%lld%lld",&N,&K);
        d.clear();
        for(int i = 1;i <= N;i++) if(N%i == 0) d.push_back(i);
        for(int i = 0;i < d.size();i++) f[i] = comb(N/d[i]+K-1,K-1);
        
        for(int i = d.size()-1;i >= 0;i--)
            for(int j = i + 1;j < d.size();j++)
                if(d[j]%d[i]==0)
                    f[i] = (f[i]-f[j]+MOD)%MOD;
        
        ans = 0;
            
        for(int i = 0;i < d.size();i++) ans=(ans+f[i]*(powmod(2,fib[d[i]])-1))%MOD;
        
        ans = ans * powmod(comb(N+K-1,K-1),MOD-2) % MOD;
        if(ans < 0) ans += MOD;
        printf("%lld\n",ans);
    }
    return 0;
} 
Last modification:October 1st, 2018 at 10:35 am
If you think my article is useful to you, please feel free to appreciate

Leave a Comment