<|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.

356. Extrasensory Perception
Time limit per test: 0.25 second(s)
Memory limit: 65536 kilobytes
input: standard
output: standard



His Royal Highness the King of Berland Berl XV has just decapitated his magician for being politically incorrect. It is a truth universally acknowledged that every king in the possession of a good fortune is in a want of a decent sorcerer. Heralds all over the country have spread the news about the coming Berland-wide wizard competition that is to end up with the announcement of the winner's name. The winner is to join the royal retinue as His Royal Highness' Counselor.

Magician Resool graduated from Bexford University ten years ago and since that time he has been chased by misfortunes and disasters. His wife, an accomplished, still pretty looking lady, believes her husband can do much better than just grumbling about the unfair life. She suggests that he should try his magic powers at the competition and reminds him of his famous trick. The gist of the trick is to see through people's souls.

After two weeks of training Resool makes himself enter the Great Hall of Berland Palace and demonstrate his exceptional abilities. His Royal Highness the King of Berland Berl XV gives him an assignment: N men and N women are presented, and the fact is known that they form N married couples; Resool is to identify a wife for every husband. After the magician casts a spell, burns some magic herbs, and performs ritual dance he still fails to answer correctly. He only manages to guess K married couples out of N. But The King is full of sympathy for the magician and can't just let the magician go away jobless and sorrowful. So, he hires him as a yard-keeper. Now Resool has a lot of free time on his new job and wants to calculate the chances he had at the competition. He wants to know the probability of guessing correctly exactly K married couples when given N men, N women, and a fact that they form N married couples. This will definitely help him during the competition next time.

Input
Input file contains two integers K and N (1 ≤ N ≤ 100; 0 ≤ K ≤ N) separated by a space.

Output
Print the required probability to the output file. If the answer is zero, simply print "0" (without quotes). Otherwise print the answer in a form of irreducible fraction "A/B" (without quotes), where A and B are positive integers without leading zeroes. See examples below for the format of output.

Example(s)
sample input
sample output
2 5
1/6

sample input
sample output
9 10
0

<|response|>
## 1) Abridged problem statement (concise)

There are \(N\) married couples. A magician randomly matches each of the \(N\) men to a distinct woman (all \(N!\) perfect matchings equally likely). Compute the probability that **exactly \(K\)** men are matched with their true wives.

Input: integers \(K, N\) \((1 \le N \le 100,\ 0 \le K \le N)\).  
Output: `0` if probability is zero; otherwise output an irreducible fraction `A/B`.

---

## 2) Key observations

1. **A matching is a permutation.**  
   Fix the “true” pairing as man \(i\) ↔ woman \(i\). Any guessed matching corresponds to a permutation \(\pi\) where man \(i\) gets woman \(\pi(i)\).

2. **Correct couples are fixed points.**  
   Man \(i\) is correctly matched iff \(\pi(i)=i\). So we need permutations of size \(N\) with exactly \(K\) fixed points.

3. **Count = choose fixed points × derange the rest.**  
   - Choose which \(K\) positions are fixed: \(\binom{N}{K}\).
   - Remaining \(m=N-K\) elements must have **no** fixed points: a **derangement** count \(D_m\).

   Hence number of valid permutations: \(\binom{N}{K} D_{N-K}\).

4. **Probability simplifies nicely.**
   \[
   P=\frac{\binom{N}{K} D_{N-K}}{N!}
   =\frac{D_{N-K}}{K!(N-K)!}
   \]
   So we only need derangements and factorials up to 100 (big integers required).

5. **Derangements recurrence (fast, exact):**
   - \(D_0=1\)
   - \(D_1=0\)
   - \(D_m=(m-1)(D_{m-1}+D_{m-2})\) for \(m\ge 2\)

---

## 3) Full solution approach

1. Read \(K, N\). Let \(m = N-K\).
2. Compute \(D_m\) using the recurrence above.
3. If \(D_m=0\), output `0` (this happens when \(m=1\), i.e. \(K=N-1\)).
4. Compute:
   - numerator \(= D_m\)
   - denominator \(= K! \cdot m!\)
5. Reduce the fraction by \(g=\gcd(\text{numerator},\text{denominator})\).
6. Output `numerator/g` + `/` + `denominator/g`.

Complexity: \(O(N)\) arithmetic operations; values are large, so use big integers.

---

## 4) C++ implementation

> Note: standard C++ has no built-in big integers, so this solution carries a self-contained arbitrary-precision `bigint` type. The derangement and factorial values for N up to 100 far exceed 64-bit range.

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

// base and base_digits must be consistent
const int base = 1000000000;
const int base_digits = 9;

struct bigint {
    vector<int> z;
    int sign;

    bigint() : sign(1) {}

    bigint(long long v) { *this = v; }

    bigint(const string& s) { read(s); }

    void operator=(const bigint& v) {
        sign = v.sign;
        z = v.z;
    }

    void operator=(long long v) {
        sign = 1;
        if(v < 0) {
            sign = -1, v = -v;
        }
        z.clear();
        for(; v > 0; v = v / base) {
            z.push_back(v % base);
        }
    }

    bigint operator+(const bigint& v) const {
        if(sign == v.sign) {
            bigint res = v;

            for(int i = 0, carry = 0;
                i < (int)max(z.size(), v.z.size()) || carry; ++i) {
                if(i == (int)res.z.size()) {
                    res.z.push_back(0);
                }
                res.z[i] += carry + (i < (int)z.size() ? z[i] : 0);
                carry = res.z[i] >= base;
                if(carry) {
                    res.z[i] -= base;
                }
            }
            return res;
        }
        return *this - (-v);
    }

    bigint operator-(const bigint& v) const {
        if(sign == v.sign) {
            if(abs() >= v.abs()) {
                bigint res = *this;
                for(int i = 0, carry = 0; i < (int)v.z.size() || carry; ++i) {
                    res.z[i] -= carry + (i < (int)v.z.size() ? v.z[i] : 0);
                    carry = res.z[i] < 0;
                    if(carry) {
                        res.z[i] += base;
                    }
                }
                res.trim();
                return res;
            }
            return -(v - *this);
        }
        return *this + (-v);
    }

    void operator*=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = 0, carry = 0; i < (int)z.size() || carry; ++i) {
            if(i == (int)z.size()) {
                z.push_back(0);
            }
            long long cur = z[i] * (long long)v + carry;
            carry = (int)(cur / base);
            z[i] = (int)(cur % base);
            // asm("divl %%ecx" : "=a"(carry), "=d"(a[i]) : "A"(cur),
            // "c"(base));
        }
        trim();
    }

    bigint operator*(int v) const {
        bigint res = *this;
        res *= v;
        return res;
    }

    friend pair<bigint, bigint> divmod(const bigint& a1, const bigint& b1) {
        int norm = base / (b1.z.back() + 1);
        bigint a = a1.abs() * norm;
        bigint b = b1.abs() * norm;
        bigint q, r;
        q.z.resize(a.z.size());

        for(int i = a.z.size() - 1; i >= 0; i--) {
            r *= base;
            r += a.z[i];
            int s1 = b.z.size() < r.z.size() ? r.z[b.z.size()] : 0;
            int s2 = b.z.size() - 1 < r.z.size() ? r.z[b.z.size() - 1] : 0;
            int d = ((long long)s1 * base + s2) / b.z.back();
            r -= b * d;
            while(r < 0) {
                r += b, --d;
            }
            q.z[i] = d;
        }

        q.sign = a1.sign * b1.sign;
        r.sign = a1.sign;
        q.trim();
        r.trim();
        return make_pair(q, r / norm);
    }

    friend bigint sqrt(const bigint& a1) {
        bigint a = a1;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        int n = a.z.size();

        int firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int norm = base / (firstDigit + 1);
        a *= norm;
        a *= norm;
        while(a.z.empty() || a.z.size() % 2 == 1) {
            a.z.push_back(0);
        }

        bigint r = (long long)a.z[n - 1] * base + a.z[n - 2];
        firstDigit = (int)sqrt((double)a.z[n - 1] * base + a.z[n - 2]);
        int q = firstDigit;
        bigint res;

        for(int j = n / 2 - 1; j >= 0; j--) {
            for(;; --q) {
                bigint r1 =
                    (r - (res * 2 * base + q) * q) * base * base +
                    (j > 0 ? (long long)a.z[2 * j - 1] * base + a.z[2 * j - 2]
                           : 0);
                if(r1 >= 0) {
                    r = r1;
                    break;
                }
            }
            res *= base;
            res += q;

            if(j > 0) {
                int d1 =
                    res.z.size() + 2 < r.z.size() ? r.z[res.z.size() + 2] : 0;
                int d2 =
                    res.z.size() + 1 < r.z.size() ? r.z[res.z.size() + 1] : 0;
                int d3 = res.z.size() < r.z.size() ? r.z[res.z.size()] : 0;
                q = ((long long)d1 * base * base + (long long)d2 * base + d3) /
                    (firstDigit * 2);
            }
        }

        res.trim();
        return res / norm;
    }

    bigint operator/(const bigint& v) const { return divmod(*this, v).first; }

    bigint operator%(const bigint& v) const { return divmod(*this, v).second; }

    void operator/=(int v) {
        if(v < 0) {
            sign = -sign, v = -v;
        }
        for(int i = (int)z.size() - 1, rem = 0; i >= 0; --i) {
            long long cur = z[i] + rem * (long long)base;
            z[i] = (int)(cur / v);
            rem = (int)(cur % v);
        }
        trim();
    }

    bigint operator/(int v) const {
        bigint res = *this;
        res /= v;
        return res;
    }

    int operator%(int v) const {
        if(v < 0) {
            v = -v;
        }
        int m = 0;
        for(int i = z.size() - 1; i >= 0; --i) {
            m = (z[i] + m * (long long)base) % v;
        }
        return m * sign;
    }

    void operator+=(const bigint& v) { *this = *this + v; }
    void operator-=(const bigint& v) { *this = *this - v; }
    void operator*=(const bigint& v) { *this = *this * v; }
    void operator/=(const bigint& v) { *this = *this / v; }

    bool operator<(const bigint& v) const {
        if(sign != v.sign) {
            return sign < v.sign;
        }
        if(z.size() != v.z.size()) {
            return z.size() * sign < v.z.size() * v.sign;
        }
        for(int i = z.size() - 1; i >= 0; i--) {
            if(z[i] != v.z[i]) {
                return z[i] * sign < v.z[i] * sign;
            }
        }
        return false;
    }

    bool operator>(const bigint& v) const { return v < *this; }
    bool operator<=(const bigint& v) const { return !(v < *this); }
    bool operator>=(const bigint& v) const { return !(*this < v); }
    bool operator==(const bigint& v) const {
        return !(*this < v) && !(v < *this);
    }
    bool operator!=(const bigint& v) const { return *this < v || v < *this; }

    void trim() {
        while(!z.empty() && z.back() == 0) {
            z.pop_back();
        }
        if(z.empty()) {
            sign = 1;
        }
    }

    bool isZero() const { return z.empty() || (z.size() == 1 && !z[0]); }

    bigint operator-() const {
        bigint res = *this;
        res.sign = -sign;
        return res;
    }

    bigint abs() const {
        bigint res = *this;
        res.sign *= res.sign;
        return res;
    }

    long long longValue() const {
        long long res = 0;
        for(int i = z.size() - 1; i >= 0; i--) {
            res = res * base + z[i];
        }
        return res * sign;
    }

    friend bigint gcd(const bigint& a, const bigint& b) {
        return b.isZero() ? a : gcd(b, a % b);
    }
    friend bigint lcm(const bigint& a, const bigint& b) {
        return a / gcd(a, b) * b;
    }

    void read(const string& s) {
        sign = 1;
        z.clear();
        int pos = 0;
        while(pos < (int)s.size() && (s[pos] == '-' || s[pos] == '+')) {
            if(s[pos] == '-') {
                sign = -sign;
            }
            ++pos;
        }
        for(int i = s.size() - 1; i >= pos; i -= base_digits) {
            int x = 0;
            for(int j = max(pos, i - base_digits + 1); j <= i; j++) {
                x = x * 10 + s[j] - '0';
            }
            z.push_back(x);
        }
        trim();
    }

    friend istream& operator>>(istream& stream, bigint& v) {
        string s;
        stream >> s;
        v.read(s);
        return stream;
    }

    friend ostream& operator<<(ostream& stream, const bigint& v) {
        if(v.sign == -1) {
            stream << '-';
        }
        stream << (v.z.empty() ? 0 : v.z.back());
        for(int i = (int)v.z.size() - 2; i >= 0; --i) {
            stream << setw(base_digits) << setfill('0') << v.z[i];
        }
        return stream;
    }

    static vector<int> convert_base(
        const vector<int>& a, int old_digits, int new_digits
    ) {
        vector<long long> p(max(old_digits, new_digits) + 1);
        p[0] = 1;
        for(int i = 1; i < (int)p.size(); i++) {
            p[i] = p[i - 1] * 10;
        }
        vector<int> res;
        long long cur = 0;
        int cur_digits = 0;
        for(int i = 0; i < (int)a.size(); i++) {
            cur += a[i] * p[cur_digits];
            cur_digits += old_digits;
            while(cur_digits >= new_digits) {
                res.push_back(int(cur % p[new_digits]));
                cur /= p[new_digits];
                cur_digits -= new_digits;
            }
        }
        res.push_back((int)cur);
        while(!res.empty() && res.back() == 0) {
            res.pop_back();
        }
        return res;
    }

    typedef vector<long long> vll;

    static vll karatsubaMultiply(const vll& a, const vll& b) {
        int n = a.size();
        vll res(n + n);
        if(n <= 32) {
            for(int i = 0; i < n; i++) {
                for(int j = 0; j < n; j++) {
                    res[i + j] += a[i] * b[j];
                }
            }
            return res;
        }

        int k = n >> 1;
        vll a1(a.begin(), a.begin() + k);
        vll a2(a.begin() + k, a.end());
        vll b1(b.begin(), b.begin() + k);
        vll b2(b.begin() + k, b.end());

        vll a1b1 = karatsubaMultiply(a1, b1);
        vll a2b2 = karatsubaMultiply(a2, b2);

        for(int i = 0; i < k; i++) {
            a2[i] += a1[i];
        }
        for(int i = 0; i < k; i++) {
            b2[i] += b1[i];
        }

        vll r = karatsubaMultiply(a2, b2);
        for(int i = 0; i < (int)a1b1.size(); i++) {
            r[i] -= a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            r[i] -= a2b2[i];
        }

        for(int i = 0; i < (int)r.size(); i++) {
            res[i + k] += r[i];
        }
        for(int i = 0; i < (int)a1b1.size(); i++) {
            res[i] += a1b1[i];
        }
        for(int i = 0; i < (int)a2b2.size(); i++) {
            res[i + n] += a2b2[i];
        }
        return res;
    }

    bigint operator*(const bigint& v) const {
        vector<int> a6 = convert_base(this->z, base_digits, 6);
        vector<int> b6 = convert_base(v.z, base_digits, 6);
        vll a(a6.begin(), a6.end());
        vll b(b6.begin(), b6.end());
        while(a.size() < b.size()) {
            a.push_back(0);
        }
        while(b.size() < a.size()) {
            b.push_back(0);
        }
        while(a.size() & (a.size() - 1)) {
            a.push_back(0), b.push_back(0);
        }
        vll c = karatsubaMultiply(a, b);
        bigint res;
        res.sign = sign * v.sign;
        for(int i = 0, carry = 0; i < (int)c.size(); i++) {
            long long cur = c[i] + carry;
            res.z.push_back((int)(cur % 1000000));
            carry = (int)(cur / 1000000);
        }
        res.z = convert_base(res.z, 6, base_digits);
        res.trim();
        return res;
    }
};

int k, n;

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

bigint derangements(int m) {
    if(m == 0) {
        return bigint(1);
    }
    if(m == 1) {
        return bigint(0);
    }

    bigint d_prev2 = 1, d_prev1 = 0;
    for(int i = 2; i <= m; i++) {
        bigint d_next = (d_prev1 + d_prev2) * (i - 1);
        d_prev2 = d_prev1;
        d_prev1 = d_next;
    }

    return d_prev1;
}

bigint factorial(int m) {
    bigint res = 1;
    for(int i = 2; i <= m; i++) {
        res *= i;
    }
    return res;
}

void solve() {
    // Classify each guessed matching by how it permutes the secret pairing.
    // Fixing exactly k of n pairs means choosing k correct positions and
    // deranging the remaining m = n - k (a permutation with no fixed point).
    // The favourable count is C(n, k) * D_m and the total is n!, so the
    // probability simplifies to D_m / (k! * m!).
    //
    // The derangement number is built with the recurrence
    // D_i = (i - 1) * (D_{i-1} + D_{i-2}), D_0 = 1, D_1 = 0. All quantities are
    // computed with the vendored arbitrary-precision bigint since n can reach
    // 100, then the fraction is reduced by its gcd.

    int m = n - k;
    bigint d = derangements(m);

    if(d.isZero()) {
        cout << 0 << '\n';
        return;
    }

    bigint numerator = d;
    bigint denominator = factorial(k) * factorial(m);

    bigint g = gcd(numerator, denominator);
    numerator /= g;
    denominator /= g;

    cout << numerator << '/' << denominator << '\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 (detailed comments)

Python has built-in big integers, so this becomes very compact and safe.

```python
from math import gcd, factorial

def derangements(m: int) -> int:
    """
    D_m = number of permutations of m elements with no fixed points.
    Recurrence:
      D_0 = 1
      D_1 = 0
      D_m = (m-1) * (D_{m-1} + D_{m-2})
    """
    if m == 0:
        return 1
    if m == 1:
        return 0

    d0, d1 = 1, 0  # D_0, D_1
    for i in range(2, m + 1):
        d0, d1 = d1, (i - 1) * (d1 + d0)
    return d1  # D_m

def main() -> None:
    # Input order: K N
    k, n = map(int, input().split())

    m = n - k                  # elements that must be non-fixed points
    num = derangements(m)      # D_{n-k}

    if num == 0:
        print("0")
        return

    # Probability = D_m / (k! * m!)
    den = factorial(k) * factorial(m)

    # Reduce fraction to irreducible form
    g = gcd(num, den)
    num //= g
    den //= g

    print(f"{num}/{den}")

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