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

探讨随机分配图书到多个书架时的美观度得分期望值计算方法。通过数学推导证明了两个关键结论,并利用组合数学及莫比乌斯反演进行高效求解。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

Description

nnn本书等概率随机放在kkk层书架上,如果第iii层书架上有cnticnt_icnti本书,那么该层书架的稳定值为stablei=f[cnti]stable_i=f[cnt_i]stablei=f[cnti],其中f[0]=0,f[1]=1,f[i]=f[i−1]+f[i−2]f[0]=0,f[1]=1,f[i]=f[i-1]+f[i-2]f[0]=0,f[1]=1,f[i]=f[i1]+f[i2]为斐波那契数列,美观度为beautyi=2stablei−1beauty_i=2^{stable_i}-1beautyi=2stablei1,此时得分为gcd(beauty1,...,beautyk)gcd(beauty_1,...,beauty_k)gcd(beauty1,...,beautyk),求得分期望值

Input

第一行一整数TTT表示用例组数,每组用例输入两个整数n,k(1≤T≤10,0&lt;n,k≤106)n,k(1\le T\le 10,0&lt;n,k\le 10^6)n,k(1T10,0<n,k106)

Output

输出得分期望值,结果模109+710^9+7109+7

Sample Input

1
6 8

Sample Output

797202805

Solution

首先证明两个结论:

1.gcd(2a−1,2b−1)=2gcd(a,b)−11.gcd(2^a-1,2^b-1)=2^{gcd(a,b)}-11.gcd(2a1,2b1)=2gcd(a,b)1

假设a≥ba\ge bab,则有

gcd(2a−1,2b−1)=gcd(2a−b(2b−1)+2a−b−1,2b−1)=gcd(2a−b−1,2b−1)...=gcd(2a%b−1,2b−1) \begin{array}{rcl} gcd(2^a-1,2^b-1)&amp;=&amp;gcd(2^{a-b}(2^b-1)+2^{a-b}-1,2^b-1)\\ &amp;=&amp;gcd(2^{a-b}-1,2^b-1)\\ &amp;...&amp;\\ &amp;=&amp;gcd(2^{a\%b}-1,2^b-1) \end{array} gcd(2a1,2b1)==...=gcd(2ab(2b1)+2ab1,2b1)gcd(2ab1,2b1)gcd(2a%b1,2b1)

辗转相除即得结论.

2.gcd(fa,fb)=fgcd(a,b)2.gcd(f_a,f_b)=f_{gcd(a,b)}2.gcd(fa,fb)=fgcd(a,b)

首先注意到
gcd(fx,fx−1)=gcd(fx−1+fx−2,fx−1)=gcd(fx−1,fx−2)=...=gcd(f1,f2)=1 gcd(f_x,f_{x-1})=gcd(f_{x-1}+f_{x-2},f_{x-1})=gcd(f_{x-1},f_{x-2})=...=gcd(f_1,f_2)=1 gcd(fx,fx1)=gcd(fx1+fx2,fx1)=gcd(fx1,fx2)=...=gcd(f1,f2)=1
假设a≥ba\ge bab,那么有
fx=fx−1+fx−2=2fx−2+fx−3=...=fy+1fx−y+fyfx−y−1 f_x=f_{x-1}+f_{x-2}=2f_{x-2}+f_{x-3}=...=f_{y+1}f_{x-y}+f_yf_{x-y-1} fx=fx1+fx2=2fx2+fx3=...=fy+1fxy+fyfxy1
则有
gcd(fx,fy)=gcd(fy+1fx−y+fyfx−y−1,fy)=gcd(fy+1fx−y,fy)=gcd(fx−y,fy)...=gcd(fx%y,fy) \begin{array}{rcl} gcd(f_x,f_y)&amp;=&amp;gcd(f_{y+1}f_{x-y}+f_yf_{x-y-1},f_y)\\ &amp;=&amp;gcd(f_{y+1}f_{x-y},f_y)\\ &amp;=&amp;gcd(f_{x-y},f_y)\\ &amp;...&amp;\\ &amp;=&amp;gcd(f_{x\%y},f_y) \end{array} gcd(fx,fy)===...=gcd(fy+1fxy+fyfxy1,fy)gcd(fy+1fxy,fy)gcd(fxy,fy)gcd(fx%y,fy)
辗转相除即得结论.

现在考虑原问题,由上面两个结论即得到当kkk个书架上书的数量为a1,...,aka_1,...,a_ka1,...,ak时,其得分即为2fgcd(a1,...,ak)−12^{f_{gcd(a_1,...,a_k)}}-12fgcd(a1,...,ak)1,总方案数即为将nnn本书放入kkk个书架,方案数Cn+k−1k−1C_{n+k-1}^{k-1}Cn+k1k1,假设gcd(a1,...,ak)=d,a1+...+ak=ngcd(a_1,...,a_k)=d,a_1+...+a_k=ngcd(a1,...,ak)=d,a1+...+ak=n的方案数有g(d)g(d)g(d)种,那么答案即为ans=∑d∣ng(d)(2fd−1)ans=\sum\limits_{d|n}g(d)(2^{f_d}-1)ans=dng(d)(2fd1),故只要求出g(1),...,g(d)g(1),...,g(d)g(1),...,g(d)即可

d∣gcd(a1,...,ak),a1+...+ak=nd|gcd(a_1,...,a_k),a_1+...+a_k=ndgcd(a1,...,ak),a1+...+ak=n的方案数为F(d)F(d)F(d),由插板法知F(d)=Cnd+k−1k−1F(d)=C_{\frac{n}{d}+k-1}^{k-1}F(d)=Cdn+k1k1

F(d)=∑d∣ig(i)F(d)=\sum\limits_{d|i}g(i)F(d)=dig(i),由莫比乌斯反演有g(d)=∑d∣iμ(id)F(i)g(d)=\sum\limits_{d|i}\mu(\frac{i}{d})F(i)g(d)=diμ(di)F(i),预处理莫比乌斯函数后直接求解即可

Code

#include<cstdio>
#include<algorithm>
#include<cstring>
using namespace std;
#define maxn 1000005
typedef long long ll;
int prime[maxn],mu[maxn],check[maxn],tot;
void Moblus(int n=1e6)
{
	memset(check,0,sizeof(check));
	mu[1]=1;
	tot=0;
	for(int i=2;i<=n;i++)
	{
		if(!check[i])
		{
			prime[tot++]=i;
			mu[i]=-1;
		}
		for(int j=0;j<tot&&i*prime[j]<=n;j++)
		{			
			check[i*prime[j]]=1;
			if(i%prime[j]==0)
			{
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
}
#define mod 1000000007
int mul(int x,int y)
{
	ll z=1ll*x*y;
	return z-z/mod*mod;
}
int add(int x,int y)
{
	x+=y;
	if(x>=mod)x-=mod;
	return x;
}
int Pow(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1)ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}
int a[maxn],fact[2*maxn],inv[2*maxn];
void init(int n=1e6)
{
	Moblus();
	a[0]=0,a[1]=1;
	for(int i=2;i<=1e6;i++)a[i]=(a[i-2]+a[i-1])%(mod-1);
	for(int i=1;i<=1e6;i++)a[i]=add(Pow(2,a[i]),mod-1);
	fact[0]=1;
	for(int i=1;i<=2*n;i++)fact[i]=mul(i,fact[i-1]);
	inv[1]=1;
	for(int i=2;i<=2*n;i++)inv[i]=mul(mod-mod/i,inv[mod%i]);
	inv[0]=1;
	for(int i=1;i<=2*n;i++)inv[i]=mul(inv[i],inv[i-1]);
}
int C(int n,int m)
{
	if(m<0||m>n)return 0;
	return mul(fact[n],mul(inv[m],inv[n-m]));
}
int F(int d,int k)
{
	return C(d+k-1,k-1);
}
int main()
{
	init();
	int T,n,k;
	scanf("%d",&T);
	while(T--)
	{
		scanf("%d%d",&n,&k);
		int ans=0;
		for(int d=1;d<=n;d++)
			if(n%d==0)
			{
				int f=0;
				for(int i=d;i<=n;i+=d)
					if(n%i==0)
					{
						if(mu[i/d]==1)f=add(f,F(n/i,k));
						else if(mu[i/d]==-1)f=add(f,mod-F(n/i,k));
					}
				ans=add(ans,mul(f,a[d]));
			}
		printf("%d\n",mul(ans,Pow(C(n+k-1,k-1),mod-2)));
	}
	return 0;
}
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值