HDU-6363：bookshelf（数论+容斥）

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 $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 $cnt_i(0≤cnt_i≤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 $stable_i=f[cnt_i]$.

   • Define the beauty value of i-th layers $beauty_i=2^{stable_i}−1$.

   • Define the whole beauty value of ZYG $score=gcd(beauty_1,beauty_2,...,beauty_k)(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≤10^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 $\frac pq$ (p and q are integer and coprime, q is positive), such that q is not divisible by $10^9+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

     思路：首先要知道2个公式：
     $1.gcd(2^a-1,2^b-1)=2^{gcd(a,b)}-1$
     $2.gcd(f[a],f[b])=f[gcd(a,b)]$
     其中$f$数组代表斐波那契数列

     那么所有的贡献是
     $ans=\sum_{gcd=1}^{gcd|n} (2^{f[gcd]}-1)*(n本书分到k个书架且最大公约数为gcd的方案数)$。

     其中$f[gcd]$可能很大，所以$2^{f[gcd]}-1$需要用欧拉降幂公式降幂。

     因为是求期望，所以$ans=\frac {ans}{总方案数}$ 。

     那么如何求n本书分到k个书架且大公约数为gcd的方案数呢？

     其中n本书分到k个书架的总方案数为$C_{n+k-1}^{k-1}$。

     那么n本书分到k个书架且$gcd$是$g$的倍数的方案数为$C_{\frac ng+k-1}^{k-1}$。

     要求出$gcd$只为 g g <script type="math/tex" id="MathJax-Element-95">g</script>的方案数就要用容斥了。

     #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;
     }