bookshelf
Time Limit: 2000/1000 MS (Java/Others)
Memory Limit: 65536/65536 K (Java/Others)
Problem Description
Patrick Star bought a bookshelf, he named it ZYG !!
Patrick Star has book .
The ZYG has layers (count from to ) and there is no limit on the capacity of each layer !
Now Patrick want to put all books on ZYG :
Assume that the i-th layer has books finally.
Assume that f[i] is the i-th fibonacci number
Define the stable value of i-th layers .
Define the beauty value of i-th layers .
Define the whole beauty value of ZYG .
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 .
Output
For each test case, output the answer as a value of a rational number modulo
.
Formally, it is guaranteed that under given constraints the probability is always a rational number (p and q are integer and coprime, q is positive), such that q is not divisible by . Output such integer a between 0 and that is divisible by .
Sample Input
1
6 8
Sample Output
797202805
思路:首先要知道2个公式:
其中
数组代表斐波那契数列
那么所有的贡献是
。
其中 可能很大,所以 需要用欧拉降幂公式降幂。
因为是求期望,所以 。
那么如何求n本书分到k个书架且大公约数为gcd的方案数呢?
其中n本书分到k个书架的总方案数为 。
那么n本书分到k个书架且 是 的倍数的方案数为 。
要求出 只为 的方案数就要用容斥了。
#include<bits/stdc++.h>
using namespace std;
const int MAX=2e6+10;
const int MOD=1e9+7;
typedef long long ll;
int f[MAX],v[MAX];
int isp[MAX],pr[MAX],all;
ll fac[MAX],inv[MAX];
ll POW(ll x,ll n,ll mod)
{
ll res=1;
while(n)
{
if(n&1)res=res*x%mod;
x=x*x%mod;
n/=2;
}
return res;
}
ll c(ll x,ll y){return fac[x]*inv[x-y]%MOD*inv[y]%MOD;}
ll cal(ll x){return (POW(2,f[x]+v[x]*(MOD-1),MOD)-1+MOD)%MOD;}//降幂
void init() //预处理
{
fac[0]=inv[0]=inv[1]=1;
for(int i=1;i<=2e6;i++)fac[i]=fac[i-1]*i%MOD;
for(int i=2;i<=2e6;i++)inv[i]=inv[MOD%i]*(MOD-MOD/i)%MOD;
for(int i=2;i<=2e6;i++)inv[i]=inv[i]*inv[i-1]%MOD;
f[0]=0,f[1]=f[2]=1;
for(int i=3;i<=2e6;i++)
{
f[i]=f[i-1]+f[i-2];
v[i]|=v[i-1];
if(f[i]>=MOD-1)v[i]=1;
f[i]%=MOD-1;
}
all = 0;
memset(isp,0,sizeof isp);
isp[1]=1;
for(int i=2;i<MAX;i++)
{
if(!isp[i])pr[all++] = i;
for(int j=0;j<all;j++)
{
long long t = 1LL*pr[j]*i ;
if(t<MAX)
{
isp[t] = 1;
if(i%pr[j]==0)break;
}
else break;
}
}
}
vector<int>p;
void Div(ll x)//分解成素因子
{
for(ll i=2;i*i<=x;i++)
{
if(x%i==0)p.push_back(i);
while(x%i==0)x/=i;
}
if(x>1)p.push_back(x);
}
ll solve(ll n,ll k)//容斥计算
{
p.clear();
Div(n);
ll tot=c(n+k-1,k-1);
for(int i=(1<<p.size())-1;i>=1;i--)
{
int tmp=1;
for(int j=0;j<p.size();j++)if(i&(1<<j))tmp*=p[j];
if(__builtin_popcount(i)%2)(tot-=c(n/tmp+k-1,k-1))%=MOD;
else (tot+=c(n/tmp+k-1,k-1))%=MOD;
(tot+=MOD)%=MOD;
}
return tot;
}
int main()
{
init();
int T;
cin>>T;
while(T--)
{
ll n,k,ans=0;
scanf("%lld%lld",&n,&k);
for(int i=1;i*i<=n;i++)//枚举公约数
{
if(n%i)continue;
(ans+=solve(n/i,k)*cal(i)%MOD)%=MOD;
if(n!=i*i)(ans+=solve(i,k)*cal(n/i)%MOD)%=MOD;
(ans+=MOD)%=MOD;
}
ans=ans*POW(c(n+k-1,k-1),MOD-2,MOD)%MOD;
printf("%lld\n",ans);
}
return 0;
}