B - bookshelf
呜呜呜,几个特性不知道,结果式子都没推出来。。
题解链接搓这里
一个放x本书的层美观函数b(x)=2Fx?1b(x)=2Fx?1
整个书架的美观值为各层美观值的最大公约数,考察任意两项:
gcd(b(x),b(y))gcd(b(x),b(y))
=gcd(2Fx?1,2Fy?1)=gcd(2Fx?1,2Fy?1)
=2gcd(Fx,Fy)?1=2gcd(Fx,Fy)?1
=2Fgcd(x,y)?1=2Fgcd(x,y)?1
以上用了gcd的两个特性。
第一个是斐波那契数列的特性,即gcd(Fx,Fy)=Fgcd(x,y)gcd(Fx,Fy)=Fgcd(x,y)
第二个是gcd(ca+1,cb+1)gcd(ca+1,cb+1)也有特性,即=cgcd(a,b)+1=cgcd(a,b)+1
ok,式子推出来之后我们就可以枚举约数,然后对约数容斥即可。
我们可以计算约数至少为g的方案数:
1,每一层都是g的倍数,所以我们只要分配系数即可。
2,系数的和为n/g,共k个系数,每个非负,方案数是个经典组合数C(n/g+k?1,k?1)C(n/g+k?1,k?1)
然后容斥就行啦。
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+7;
const int phimod=1e9+6;
const int maxn=2e6+5;long long fmod(long long a)
{return a%phimod+((a>=phimod)?phimod:0);
}//***************************************************
//返回d=gcd(a,b);和对应于等式ax+by=d中的x,y
long long extgcd(long long a,long long b,long long &x,long long &y)
{if(a==0&&b==0)return -1;//无最大公约数if(b==0){x=1;y=0;return a;}long long d=extgcd(b,a%b,y,x);y-=a/b*x;return d;//返回gcd(a,b)
}
//****************求逆元******************************
//ax=1(mod n)
long long mod_reverse(long long a,long long n)
{long long x,y;long long d=extgcd(a,n,x,y);if(d==1)return (x%n+n)%n;return -1;
}long long f[maxn];
long long tab1[maxn],tab2[maxn];
void init()
{f[1]=1;for(int i=2;i<maxn;i++)f[i]=fmod(f[i-1]+f[i-2]);/*for(int i=1;i<=10;i++)cout<<i<<' '<<f[i]<<endl;*/tab1[0]=1;for(int i=1;i<maxn;i++)tab1[i]=tab1[i-1]*i%mod;tab2[maxn-1]=mod_reverse(tab1[maxn-1],mod);for(int i=maxn-2;i>=0;i--)tab2[i]=tab2[i+1]*(i+1)%mod;
}long long quick_mul(long long a,long long b)
{long long res=1;while(b){if(b%2)res=res*a%mod;b/=2;a=a*a%mod;}return res;
}
vector<int>V;
int n,k;long long C(long long a,long long b)
{return tab1[a]*tab2[b]%mod*tab2[a-b]%mod;
}
long long cnt[maxn];
void solve()
{int e=(int)sqrt(0.0+n);for(int i=1;i<=e;i++)if(n%i==0){if(i*i==n)V.push_back(i);else V.push_back(i),V.push_back(n/i);}sort(V.begin(),V.end());for(int i=0;i<V.size();i++){int now=V[i];cnt[i]=C(n/now+k-1,k-1);}long long sum=cnt[0];for(int i=V.size()-1;i>=0;i--)for(int j=i+1;j<V.size();j++)if(V[j]%V[i]==0)cnt[i]=((cnt[i]-cnt[j])%mod+mod)%mod;//cout<<quick_mul(8,6)<<endl;long long res=0;for(int i=0;i<V.size();i++){long long tmp=cnt[i];tmp=(tmp*(quick_mul(2,f[V[i]])%mod-1)%mod+mod)%mod;res=(res+tmp)%mod;}printf("%lld\n",res*mod_reverse(sum,mod)%mod);
}int main()
{init();int T;scanf("%d",&T);while(T--){scanf("%d %d",&n,&k);V.clear();solve();}return 0;}