## 1. Abridged Problem Statement

Given **N distinct** photo albums and **K friends**, count the number of ways to distribute all albums so that **each friend gets at least one album**.

Two distributions are considered the same if they differ only by:
- the order of albums within a friend's gift, or
- the labeling/order of the friends (i.e., friends are indistinguishable).

Output the answer **modulo 2007**.

---

## 2. Detailed Editorial

### Key interpretation

A "distribution" is a partition of the **N labeled albums** into **K non-empty unlabeled groups** (since friends are indistinguishable). This is precisely the **Stirling number of the second kind** S(N, K), and we need S(N, K) mod 2007.

### DP recurrence

The standard recurrence for Stirling numbers is:

S(n, k) = S(n-1, k-1) + k * S(n-1, k)

Intuition by placing album n:
- It forms a new group by itself: S(n-1, k-1)
- It joins one of the existing k groups: k * S(n-1, k)

Equivalently, if P is the number of people already holding at least one album, giving album n to one of the P existing recipients keeps the count, while giving it to a new person increments P. Only the case where ultimately P = K counts.

Base conditions: S(0, 0) = 1; S(n, 0) = 0 for n > 0; S(n, k) = 0 for n < k.

Because N can be very large, the naive O(NK) DP is too slow.

### Matrix exponentiation

For fixed K, the transition from n-1 to n is **linear** with constant coefficients. Define the state vector:

v_n[p] = S(n, p), for p = 0..K

From the recurrence S(n, p) = p * S(n-1, p) + S(n-1, p-1), this is a linear transformation v_n = T * v_{n-1}, where T is a (K+1) x (K+1) matrix:

- T[p][p] = p  (for p >= 1)
- T[p+1][p] = 1  (for p <= K-1)

Then v_N = T^N * v_0, where v_0 = [1, 0, 0, ..., 0]^T (since S(0, 0) = 1). The desired answer is S(N, K) = (T^N)[K][0].

Compute T^N mod 2007 via fast matrix exponentiation in O(K^3 log N), which is efficient since K is small.

### Alternative: inclusion-exclusion

There is also a closed form via inclusion-exclusion. If friends were distinguishable, the count (without the "at least one" constraint) would be K^N. By inclusion-exclusion:

answer = sum_{j=0}^{K} C(K, j) * (-1)^j * (K - j)^N

This equals K! * S(N, K). The code implements the matrix exponentiation approach (which directly computes S(N, K)).

### Edge case

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

---

## 3. C++ Solution

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

---

## 4. Python Solution

```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()
```

---

## 5. Compressed Editorial

We need the number of partitions of N distinct albums into K non-empty unlabeled groups, i.e. S(N, K) (Stirling numbers of the 2nd kind), modulo 2007. If N < K, answer is 0.

Use recurrence S(n, k) = S(n-1, k-1) + k * S(n-1, k). Represent [S(n,0), ..., S(n,K)]^T as a vector and build a constant transition matrix T of size K+1:
- T[p][p] = p for p >= 1
- T[p+1][p] = 1 for p <= K-1

Then v_N = T^N * v_0 with v_0 = [1, 0, ..., 0]^T. Output (T^N)[K][0]. Compute T^N via fast exponentiation in O(K^3 log N), all mod 2007.
