HDU 6372 sacul（组合数学）

Description

定义$p^i\times p^i$矩阵$HMBB_i$，已知$HMBB_n[i][j]=(C(i,j)\ mod\ p>0)?1:0$，而$F[n][k]$表示$(HMBB_n)^k$元素之和，求$\sum\limits_{i=1}^n\sum\limits_{j=1}^kF[i][j]\ mod\ 10^9+7$

Input

第一行一整数$T$表示用例组数，每组用例输入三个整数$c,n,k$，其中$c$表示$p$为第$c$个素数

$(0<n\le 10^9,0<c,k\le 10^5)$

Output

输出结果

Sample Input

1
1 1 1

Sample Output

3

Solution

首先由$Lucas$定理$C(n,m)\%p=(C(n/p,m/p)\cdot C(n\%p,m\%p))\%p$，这意味着$C(n,m)\%p$的值等于将$n,m$拆成$p$进制后，每位做组合数的结果乘起来再模$p$，注意到$C(a,b)\%p>0,0\le b\le a<p$，故若$C(n,m)\%p>0$，则$n$的$p$进制每一位都不小于$m$的$p$进制这一位的取值

之后将$HMBB$矩阵看作一个$p^i$点的图的邻接矩阵，$HMBB[i][j]=1$说明$i$可一步到达$j$，而可达性等价于$i$的$p$进制表示每一位不小于$j$的$p$进制表示这一位的取值，那么$F[i][j]$即为找$j$个有$i$位$p$进制的数使得这$j$个数字依次可达的方案数，由于$i$位相互独立，每一位的问题变成找$j$个数$a_1,...,a_j$满足$0\le a_1\le a_2\le ...\le a_j<p$，由插板法知方案数为$C_{j+p}^{p-1}$，故答案即为$\sum\limits_{i=1}^n\sum\limits_{j=1}^k(C_{j+p}^{p-1})^i=\sum\limits_{j=1}^k\sum\limits_{i=1}^n(C_{j+p}^{p-1})^i$，内部求和用等比数列求和优化一下即可，时间复杂度$O(klogp)$

Code

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
#define maxn 1500005
int fact[maxn],inv[maxn],sum[maxn];
int p[2000005],res=0,mark[2000005];
void prime(int n=2e6)
{
    for(int i=2;i<=n;i++)
        if(!mark[i])
        {
            p[++res]=i;
            for(int j=2*i;j<=n;j+=i)mark[j]=1;
        }
}
#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;
}
void init(int n=1500000)
{
    fact[0]=1;
    for(int i=1;i<=n;i++)fact[i]=mul(i,fact[i-1]);
    inv[1]=1;
    for(int i=2;i<=n;i++)inv[i]=mul(mod-mod/i,inv[mod%i]);
    sum[0]=1;
    for(int i=1;i<=n;i++)sum[i]=mul(inv[i],sum[i-1]);
}
int C(int n,int m)
{
    if(m<0||m>n)return 0;
    return mul(fact[n],mul(sum[m],sum[n-m]));
}
int deal(int p,int n)
{
    if(p==1)return n;
    int ans=add(Pow(p,n),mod-1);
    ans=mul(ans,p);
    ans=mul(ans,Pow(p-1,mod-2));
    return ans;
}
int main()
{
    prime();
    init();
    int T,n,k,c;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d",&c,&n,&k);
        int ans=0;
        for(int i=1;i<=k;i++)
        {
            int t=C(i+p[c],p[c]-1);
            int temp=deal(t,n);
            ans=add(ans,temp);
        } 
        printf("%d\n",ans);
    }
    return 0;
}