SGU - 441 Set Division 斯特灵数

题目：N个物品分成K个非空子集的方法数。

思路：第二类斯特灵数。N很大，K很小，所以用矩阵做。

S(N,K)=K*S(N-1,K)+S(N-1.K-1)

代码：

#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<iostream>
#include<algorithm>
#include<ctime>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<string>
#include<vector>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<list>
#include<numeric>
using namespace std;
#define LL long long
#define ULL unsigned long long
#define INF 0x3f3f3f3f3f3f3f3f
#define mm(a,b) memset(a,b,sizeof(a))
#define PP puts("*********************");
template<class T> T f_abs(T a){ return a > 0 ? a : -a; }
template<class T> T gcd(T a, T b){ return b ? gcd(b, a%b) : a; }
template<class T> T lcm(T a,T b){return a/gcd(a,b)*b;}
// 0x3f3f3f3f3f3f3f3f

const int MOD=2007;
const int MAXN = 13;
int sz;
struct Matrix {
	int a[MAXN][MAXN];
	void toZero() {
		for(int i = 0; i < sz; ++i)
			for(int j = 0; j < sz; ++j)
				a[i][j] = 0;
	}
	void toOne() {
		for(int i = 0; i < sz; ++i) {
			for(int j = 0; j < sz; ++j)
				a[i][j] = (int)(i == j);
		}
	}
	Matrix operator + (const Matrix &B) const {
		Matrix A(*this);
		int i, j;
		for(i = 0; i < sz; ++i) {
			for(j = 0; j < sz; ++j)
				//A.a[i][j] += B.a[i][j];
				A.a[i][j] = (A.a[i][j] + B.a[i][j]) % MOD;
		} return A;
	}
	Matrix operator * (const Matrix &B) const {
		Matrix A(*this), res;
		res.toZero();
		int i, j, k;
		for(i = 0; i < sz; ++i) {
			for(j = 0; j < sz; ++j) {
				if(A.a[i][j] == 0) continue;
				for(k = 0; k < sz; ++k)
					//res.a[i][k] += A.a[i][j] * B.a[j][k];
					res.a[i][k] = (res.a[i][k] + ((LL)A.a[i][j] * B.a[j][k]) % MOD) % MOD;
			}
		} return res;
	}
	Matrix operator ^ (int p) const {
		Matrix res, X(*this);
		res.toOne();
		for(; p; p >>= 1) {
			if(p & 1) res = res * X;
			X = X * X;
		} return res;
	}
	void Disp() {
		printf("---------------------\n");
		for(int i = 0; i < sz; ++i) {
			for(int j = 0; j < sz; ++j)
				printf("%4d", a[i][j]);
			printf("\n");
		}
		printf("---------------------\n");
	}
};
int main(){

    int N,K;
    while(~scanf("%d%d",&N,&K)){
        sz=K+1;
        Matrix A;
        A.toZero();
        A.a[0][0]=1;
        for(int i=1;i<=K;i++){
            A.a[i][i]=i;
            A.a[i][i-1]=1;
        }
        A=A^(N-1);
        printf("%d\n",A.a[K][1]);
    }
    return 0;
}