<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

441. Set Division
Time limit per test: 0.25 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard



Ruslan has K friends. And all of them have birthday tomorrow. He has bought N different photo albums yesterday, and wants to present these photo albums to his friends. Of course, he can't give less than 1 photo album to anybody. Your task is to calculate how many possible ways is there to do it.

All photo albums are different. Two distributions are considered the same if they differ only by order of albums in the gifts or by the persons receiving gifts.

Input
In the only line of the input there are two numbers separated by a space — , .

Output
Output should contain one number — the number of possible ways to distribute the albums modulo 2007.

Example(s)
sample input
sample output
3 2
3

<|response|>
## 1. Abridged Problem Statement

You have **N distinct** photo albums and want to give them to **K friends** so that **each friend gets at least one album**.

Two distributions are considered the same if they differ only by:
- the order of albums inside each gift bundle, and/or
- which friend got which bundle (friends are **indistinguishable**).

Compute the number of such distributions **modulo 2007**.

---

## 2. Key Observations

1. Since albums are distinct but friends are indistinguishable, the result is the number of ways to partition an **N-element labeled set** into **K non-empty unlabeled subsets**.

2. This number is exactly the **Stirling number of the second kind**:
   S(N, K)

3. Classic recurrence:
   S(n, k) = S(n-1, k-1) + k * S(n-1, k)

   Intuition: place album n either as a new singleton group (S(n-1, k-1)) or into one of the existing k groups (k * S(n-1, k)).

4. If N < K, it is impossible to give everyone at least one album, so the answer is **0**.

5. N can be very large while K is relatively small, so O(NK) DP is too slow. The recurrence is **linear** in the vector [S(n,0), ..., S(n,K)], so we can use **matrix exponentiation** in O(K^3 log N).

---

## 3. Full Solution Approach

### Step A: Model as Stirling numbers

We need S(N, K) mod 2007.

### Step B: Turn recurrence into a linear transformation

Define a state vector v_n = [S(n,0), S(n,1), ..., S(n,K)]^T.

From S(n,p) = S(n-1,p-1) + p * S(n-1,p), we can write v_n = T * v_{n-1} where T is a (K+1) x (K+1) matrix:
- Diagonal term: T[p][p] = p (for p >= 1)
- Subdiagonal term: T[p][p-1] = 1 (equivalently T[p+1][p] = 1)

### Step C: Exponentiate

Base vector: v_0 = [1, 0, ..., 0]^T (since S(0,0) = 1).

Thus v_N = T^N * v_0 and the answer is S(N,K) = (T^N)[K][0].

### Step D: Edge case

If N < K: print 0.

---

## 4. C++ Implementation with Detailed Comments

```cpp
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

const int mod = 2007;

int64_t n;
int k;

using matrix = vector<vector<int64_t>>;

matrix multiply(const matrix& a, const matrix& b) {
    int sz = a.size();
    matrix c(sz, vector<int64_t>(sz, 0));
    for(int i = 0; i < sz; i++) {
        for(int j = 0; j < sz; j++) {
            for(int l = 0; l < sz; l++) {
                c[i][j] = (c[i][j] + a[i][l] * b[l][j]) % mod;
            }
        }
    }
    return c;
}

matrix mat_pow(matrix m, int64_t p) {
    int sz = m.size();
    matrix result(sz, vector<int64_t>(sz, 0));
    for(int i = 0; i < sz; i++) {
        result[i][i] = 1;
    }
    while(p > 0) {
        if(p & 1) {
            result = multiply(result, m);
        }
        m = multiply(m, m);
        p >>= 1;
    }
    return result;
}

void read() { cin >> n >> k; }

void solve() {
    // We can see that N is very large while K is small, so we should think of
    // either a closed form. Let's start with a simplification by making the
    // friends indistinguishable - we can think of this as sorting them by their
    // lowest album ID. After making this reduction, we can notice that each
    // sorted configuration is mapped to exactly K! of the original. The main
    // reason to make this change is that it let's us forget the ordering of the
    // friends and makes the transitions easier.
    //
    // We will start giving albums from 1 until N. Let's say we have P people
    // that have already received an album. We now have only two options -
    // give the current album to any of the people that already have albums in P
    // ways, or give it to someone that still hasn't received an album (bool)(K
    // - P) ways. This means we have a direct dp[N][P] solution and furthermore
    // each state depends only on the previous one. This is already enough to
    // create a O(K^3 log N) solution, because we can directly express dp[N][.]
    // as a linear combination of dp[N - 1][.] as the coefficients don't depend
    // on N. Then we can do binary exponentiation.
    //
    // There is also an alternative approach based on inclusion exclusion. Let's
    // drop the symmetry division by K!, and instead solve the original problem.
    // If we didn't have the "at least 1 album" constraint, the solution would
    // have been K^n. However, we can do inclusion-exclusion here. Particularly,
    // there are C(K, j) ways to have at least j people receive 0 albums, and
    // then (K - j)^n to distribute the rest. This leads us to the closed form
    // of sum_{j=0}^{K} C(K, j) * (-1)^j * (K - j)^n. With a precompute of the
    // factorials up to K, we can solve this in just O(K log N).
    //
    // People familiar with number theory might notice the problem here asks for
    // k! * S(n, k), or the Stirling numbers of the second kind. The two
    // approaches above are precisely the ways to compute it: due to a recursive
    // relation, and with inclusion-exclusion. Here we implement the first
    // solution.

    if(n < k) {
        cout << 0 << "\n";
        return;
    }

    matrix trans(k + 1, vector<int64_t>(k + 1, 0));
    for(int p = 0; p <= k; p++) {
        if(p > 0) {
            trans[p][p] = p % mod;
        }
        if(p < k) {
            trans[p + 1][p] = 1;
        }
    }

    matrix result = mat_pow(trans, n);
    cout << result[k][0] << "\n";
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    // cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}
```

---

## 5. Python Implementation with Detailed Comments

```python
import sys

MOD = 2007

def mat_mul(a, b):
    """Multiply two square matrices modulo MOD."""
    n = len(a)
    c = [[0] * n for _ in range(n)]
    for i in range(n):
        ai = a[i]
        for k in range(n):
            if ai[k] == 0:
                continue
            aik = ai[k]
            bk = b[k]
            for j in range(n):
                c[i][j] = (c[i][j] + aik * bk[j]) % MOD
    return c

def mat_pow(m, p):
    """Raise square matrix m to power p modulo MOD using binary exponentiation."""
    n = len(m)

    # Identity matrix
    res = [[0] * n for _ in range(n)]
    for i in range(n):
        res[i][i] = 1

    while p > 0:
        if p & 1:
            res = mat_mul(res, m)
        m = mat_mul(m, m)
        p >>= 1
    return res

def solve():
    data = sys.stdin.read().strip().split()
    if not data:
        return
    N = int(data[0])
    K = int(data[1])

    # Impossible to split N distinct albums into K non-empty groups if N < K
    if N < K:
        print(0)
        return

    # Transition matrix for Stirling numbers S(n,k):
    # S(n,p) = S(n-1,p-1) + p*S(n-1,p)
    # State vector is [S(n,0), S(n,1), ..., S(n,K)]^T
    size = K + 1
    T = [[0] * size for _ in range(size)]

    for p in range(size):
        if p > 0:
            T[p][p] = p % MOD      # coefficient of S(n-1,p)
        if p < K:
            T[p + 1][p] = 1        # coefficient of S(n-1,p) flowing to S(n,p+1)

    # Compute T^N. Since v0 = [1,0,0,...]^T, S(N,K) is (T^N)[K][0].
    TN = mat_pow(T, N)
    print(TN[K][0] % MOD)

if __name__ == "__main__":
    solve()
```

These solutions directly compute S(N,K) mod 2007 using matrix exponentiation, which is efficient when N is large and K is small.
