小球与盒子（矩阵乘法+快速幂）

Description of the Problem:

你有K个相同的盒子，N个互不相同的物品。你准备把这N个物品装入K个盒子，每个盒子最少要放入一个物品。问一共会有多少种分配方法。由于方案数很大，只需要输出方案总数除以10000的余数。

INPUT：

第一行有一个正整数 t ，表示数据组数(不多于50)。每组数据仅一行，两个整数， N 和K，其中1≤N ≤ 10^9，K≤min(50,N)。

OUTPUT：

每行输出一个整数，为方案总数除以10000的余数。

Analysis:

　　这是一个来自于同为HUSTer的高中童鞋的问题,从这个问题描述中,我们可以知道把N个物品放入K个相同的盒子的方法数就是 把N-1个物品放入K-1个盒子的方法数 加上 把N-1个物品放入K个盒子的方法数*K,为什么?我们只要考虑最后一个盒子和最后一个球即可,由于盒子是相同的,那么有两种情况
　　1.最后一个盒子只有一个球且为第N个球，即前面N-1个球没有填满K个盒子，最后一个球只能放在空盒子里。这种情况方法数就等于把N-1个物品放入K-1个盒子的方法数。
　　2.前面N-1个球填满的K个盒子，那么最后一个球可以任意放。这种情况方法数就等于 把N-1个物品放入K个盒子的方法数*K。

　　所以如果用f(N,K)来表示 N个物品K个盒子的方法数，可以得到以下递推式

　　f(N,K)=f(N-1,K-1)+f(N-1,K)*K;(即第二类斯特林数)

　　所以现在已经可以用程序通过递推算出结果（递归也可以但要慢一些，即使是记忆化递归）。但问题解决了吗？

　　答案是：NO。由于N很大所以要一项一项递推计算f(N,K)肯定会超时应为至少要算N*K次，所以显然要找更快速的算法让时间复杂度降低O(logn*K)。我们这时可以运用矩阵+快速幂的形式来解决这个问题。

　　首先来解释一下什么是快速幂：举个例子，对于计算q^k，我们其实可以将之看成q^[a(n)*2^n+a(n-1)*2^(n-1)+……a(1)*2^1+a(0)*2^0],其中a(n)=0或1。即把k看成许多{2^m}中某n+1项之和，那么怎么快速计算q^k呢？我们先把q^k看成n+1项之积，然后我们可以对k不断除2取余，如果第i次除2后余数为1，那么表示a(i-1)=1，即有n+1项乘积中有q^(i-1)这一项，具体描述步骤如下

　　首先我们用让一个数ans=1，用来存答案，用一个数x来存q^(i-1)，x开始为1，然后先对k进行第一次除以2，如果余数为1，证明a(0)=1，有q^0这一项，于是让ans*=x,如果不为1，证明不存在这一项(或者说着这一项为1)，就不用乘，然后进行完上述操作，x*=q;
　　然后我们对k进行第二次除以2，此时x=q^1，如果余数为1，证明a(1)=1，有q^1这一项，于是让ans*=x,如果不为1，证明不存在这一项(或者说着这一项为1)，就不用乘，然后进行完上述操作，x*=q;所以在对k进行第i次除以2操作时，x=q^(i-1)， 如果余数为1，证明a(i-1)=1，有q^(i-1)这一项，于是让ans*=x,如果不为1，证明不存在这一项(或者说着这一项为1)，就不用乘，然后进行完上述操作，x*=q;最后如果k=0证明已经除完了，ans已经算完那么就得到了结果，结束运算。
　　现在来分析一下算法效率，若按照原算法计算q^k要乘k次也就是计算k次，用现在的快速幂算法，则要除以[log2(k)]+1次，也就是x要乘[log2(k)]+1次，而如果a(n)均为1，那么要ans也要乘 [log2(k)]+1，所以最多算2{[log2(k)]+1}次算法效率为O(logk)。
现在来讲一讲矩阵乘法：
　　定义：设A为m×p的矩阵，B为p×n的矩阵，那么称m×n的矩阵C为矩阵A与B的乘积，记作C=A×B
　　计算方法：乘积C的第m行第n列的元素等于矩阵A的第m行的元素与矩阵B的第n列对应元素乘积之和。
性质：乘法结合律：(AB)C=A(BC)

　　所以我们可以把原问题的递推式子变成一个矩阵递推式

⎛⎝⎜⎜⎜⎜⎜⎜⎜f(n,k)f(n,k−1)⋮f(n,2)f(n,1)⎞⎠⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜k0⋮001k−1⋮0001⋮00⋯⋯⋱⋯⋯00⋮2000⋮11⎞⎠⎟⎟⎟⎟⎟⎟⎟⎛⎝⎜⎜⎜⎜⎜⎜⎜f(n−1,k)f(n−1,k−1)⋮f(n−1,2)f(n−1,1)⎞⎠⎟⎟⎟⎟⎟⎟⎟

　　所以可以用矩阵递推式把原递推式变成等比数列形式那么可以得到以下式子

⎛⎝⎜⎜⎜⎜⎜⎜⎜f(n,k)f(n,k−1)⋮f(n,2)f(n,1)⎞⎠⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜⎜⎜⎜⎜k0⋮001k−1⋮0001⋮00⋯⋯⋱⋯⋯00⋮2000⋮11⎞⎠⎟⎟⎟⎟⎟⎟⎟(n−k)⎛⎝⎜⎜⎜⎜⎜⎜⎜f(k,k)f(k,k−1)⋮f(k,2)f(k,1)⎞⎠⎟⎟⎟⎟⎟⎟⎟

　　然后我们将快速幂与矩阵乘法递推式结合起来，我们只需要先用递推计算出base矩阵，然后用快速幂算出per矩阵的(n-k)次方两者相乘得到ans 其中我们将ans矩阵不再初始化为1，而是初始化为矩阵中的“1”(其他矩阵相乘等于其他)。

⎛⎝⎜⎜⎜⎜⎜⎜⎜10⋮0001⋮0000⋮00⋯⋯⋱⋯⋯00⋮1000⋮01⎞⎠⎟⎟⎟⎟⎟⎟⎟

Code:

#include<iostream>#include<cstring>using namespace std;const int mod=10000;struct matrax{    long long m[66][66];};matrax per,ans,base,cal,unit; matrax mul(matrax a,matrax b,int x,int y,int z){    matrax c;    for (int i=0;i<x;i++)    for (int j=0;j<z;j++)    {        c.m[i][j]=0;        for (int k=0;k<y;k++)            c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j])%mod;        c.m[i][j]=c.m[i][j]%mod;    }    return c;}matrax qpow(int n,int k){    matrax p=per;    int ci=n-k;    while (ci)    {        if (ci&1)         {           unit=mul(unit,p,k,k,k);           ci--;        }        else         {            p=mul(p,p,k,k,k);            ci>>=1;        }    }    return unit;}void init(int k){    for (int i=0;i<50;i++)     for (int j=0;j<50;j++)    {        if (i==j&&i<k) per.m[i][j]=k-i;        else if (j-1==i&&j<k) per.m[i][j]=1;        else per.m[i][j]=0;        if(i==j) unit.m[i][j]=1;        else unit.m[i][j]=0;    }    for (int i=0;i<k;i++)    base.m[i][0]=cal.m[k-1][k-1-i];}void calculation(){    for (int i=0;i<50;i++) cal.m[i][i]=cal.m[i][0]=1;    for(int i=1;i<50;i++)    for(int j=1;j<i;j++)    cal.m[i][j]=(cal.m[i-1][j]*(j+1)+cal.m[i-1][j-1])%mod;}int main(){    int t,n,k;    cin>>t;    calculation();    while (t--)    {        cin>>n>>k;        init(k);        matrax xi=qpow(n,k);        ans=mul(xi,base,k,k,1);        cout<<ans.m[0][0]<<endl;    }    return 0;}

Addition：

　　下面给出一部分结果

n k 1 2 3 4 5 6 7 8 1 1 2 1 1 3 1 3 1 4 1 7 6 1 5 1 15 25 10 1 6 1 31 90 65 15 1 7 1 63 301 350 140 21 1 8 1 127 966 1701 1050 266 28 1