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 :
Assume that the i-th layer has cnti(0≤cnti≤N) books finally.
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]).
Define the stable value of i-th layers stablei=f[cnti].
Define the beauty value of i-th layers beautyi=2stablei−1.
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)
(
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层的书架上,每一层的美丽值为 bi=2f[cnt]−1 b i = 2 f [ c n t ] − 1 ,其中cnt是这一层书的数量,f[x]为斐波那契数列,整个书架的美丽值为 gcd(b1,b2,...,bk) g c d ( b 1 , b 2 , . . . , b k ) ,问整个书架的美丽值的期望
分析:
首先需要利用数论中的两个结论进行化简
gcd(xa−1,xb−1)=xgcd(a,b)−1 g c d ( x a − 1 , x b − 1 ) = x g c d ( a , b ) − 1
gcd(f[x],f[y])=f[gcd(x,y)] g c d ( f [ x ] , f [ y ] ) = f [ g c d ( x , y ) ] ( f[] f [ ] 是斐波那契数列)
于是可以进一步化简公式为 2f[gcd(cnt1,cnt2,...,cntk)]−1 2 f [ g c d ( c n t 1 , c n t 2 , . . . , c n t k ) ] − 1
这样只需要枚举gcd,计算gcd出现次数就可以知道gcd的贡献值了
设 f(g)表示gcd(x1,x2,...xk)=g的个数 f ( g ) 表 示 g c d ( x 1 , x 2 , . . . x k ) = g 的 个 数
F(p)表示gcd(x1,x2,...xk)=g的倍数p的个数 F ( p ) 表 示 g c d ( x 1 , x 2 , . . . x k ) = g 的 倍 数 p 的 个 数
我们发现对于 f(g) f ( g ) 来说不是很好求,但是 F(p) F ( p ) 好求
如果想知道n内有多少个g的倍数的话那么应该共有 np个 n p 个
那么将这些数分成k份的话,利用隔板法
|∗|∗|...|∗∗|∗∗| | ∗ | ∗ | . . . | ∗ ∗ | ∗ ∗ |
k+1 k + 1 个板子可以得到k个盒子,中间盒子中任意放 np个 n p 个 数
那么相当于中间 np+k−1 n p + k − 1 个位置的组合,其中选取k-1个位置,那么其他位置便确定了
因此 F(p)=Ck−1np+k−1 F ( p ) = C n p + k − 1 k − 1
而gcd的贡献值应该为 f(gcd)×(2f[gcd]−1) f ( g c d ) × ( 2 f [ g c d ] − 1 )
这样我们知道了 F() F ( ) 就可以通过莫比乌斯反演计算 f() f ( )
因此可以先在 n的范围内枚举gcd(n|gcd),然后在从gcd为基础枚举gcd的倍数p(p|n) n 的 范 围 内 枚 举 g c d ( n | g c d ) , 然 后 在 从 g c d 为 基 础 枚 举 g c d 的 倍 数 p ( p | n )
f(gcd)=∑gcd|pμ(pgcd)F(p)=∑gcd|pμ(pgcd)Ck−1np+k−1 f ( g c d ) = ∑ g c d | p μ ( p g c d ) F ( p ) = ∑ g c d | p μ ( p g c d ) C n p + k − 1 k − 1
这样求出 f() f ( ) 后可以乘上值 (2f[gcd]−1) ( 2 f [ g c d ] − 1 )
所以答案为
ans=∑gcd|n2f[gcd]−1∑gcd|pμ(pgcd)Ck−1np+k−1 a n s = ∑ g c d | n 2 f [ g c d ] − 1 ∑ g c d | p μ ( p g c d ) C n p + k − 1 k − 1
因为要求期望,所以再除去总的可能数
ans=ansCk−1n+k−1 a n s = a n s 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;
}