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

407. Number of Paths in the Empire
Time limit per test: 0.75 second(s)
Memory limit: 65536 kilobytes
input: standard
output: standard

During the period of Tsam dynasty ruling people rarely fought against each other and their neighbours, because they just could not afford themselves such a waste of time. The matter for it is that they were entirely absorbed with solving various problems which were connected with trade, craft, agriculture and other spheres of human activity. In a wide list of problems the ones of tax collection stand out. As one of such problems was posed by Emperor himself, it was of great importance. The problem was to count the number of different paths consisting of exactly m roads. Every path should have started and ended in the capital of Empire. Paths were supposed to cover the same towns and roads any times, moreover they could cover the capital several times. Now you are to solve this problem given information about Empire: there were n country towns situated at the foot of a hill, they formed a circle at the bottom, and the capital was on the top of the hill. The capital was connected to all other towns, and each country town was also connected to other two neighbouring country towns both to the left and to the right. Pic. 1 Empire comprising the capital (index 0) and four country towns (indices 1 — 4).
Input
The only line of input file contains two integers n and m (3 ≤ n ≤ 1000, 0 ≤ m ≤ 5000).
Output
Output the answer without leading zeros.
Example(s)
sample input
sample output
4 3
8

sample input
sample output
3 4
21

Commentary to the first sample test. There are 8 paths in the Empire. 0-1-2-0, 0-2-3-0, 0-3-4-0, 0-4-1-0, 0-2-1-0, 0-3-2-0, 0-4-3-0, 0-1-4-0.

<|response|>
1. Abridged Problem Statement
We have a "wheel" graph with one capital at node 0 and n peripheral towns (nodes 1…n) forming a cycle. The capital is connected to every town, and each town is also connected to its two cycle‐neighbors. Count the number of walks of exactly m edges that start and end at the capital (node 0), revisiting nodes and edges is allowed. Output the exact count (no modulo).

Constraints:
• 3 ≤ n ≤ 1000
• 0 ≤ m ≤ 5000

2. Key Observations
• All peripheral towns are symmetric, so instead of tracking each of the n towns separately, merge them into one "supernode."
• We only need to track two macro‐states at step k:
  – State 0: you are at the capital.
  – State 1: you are at any peripheral town (the supernode).
• Compute the number of ways to transition between these two states in one step:
  – From 0 → super: you can pick any of the n towns → n ways
  – From 0 → 0: impossible → 0 ways
  – From super → 0: aggregated as 1 way (choosing the spoke back)
  – From super → super: move along the cycle to one of two neighbors → 2 ways
• This yields a 2×2 transition matrix
   A = [ [0, 1],
       [n, 2] ]
  so that if v_k = [ways at 0 after k steps, ways at super after k steps], then
   v_{k+1} = A · v_k.
• We start with v₀ = [1, 0]. After m steps, v_m = A^m · v₀, and the answer is the first component of v_m.
• Compute A^m by fast exponentiation in O(log m) matrix multiplications. The counts grow huge and there is no modulus, so use big‐integer arithmetic (a vendored bigint struct in C++ or Python's built‐in int).

3. Full Solution Approach
1. Build the 2×2 matrix A = [[0,1],[n,2]].
2. Fast‐exponentiate A to the power m using binary exponentiation:
   • Initialize result = identity matrix I₂.
   • While p > 0:
     – If p odd, result = result × A
     – A = A × A
     – p >>= 1
3. Multiply the resulting A^m by the initial vector v₀ = [1,0]^T. The answer is the first entry of v_m = A^m·v₀, which equals (A^m)[0][0].
4. Print that big integer.

4. C++ Implementation
```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;
    }
};

using Mat = array<array<bigint, 2>, 2>;

Mat mat_mult(const Mat& a, const Mat& b) {
    Mat c;
    for(int i = 0; i < 2; i++) {
        for(int j = 0; j < 2; j++) {
            c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j];
        }
    }

    return c;
}

Mat mat_pow(Mat a, int64_t p) {
    Mat res;
    res[0][0] = 1;
    res[0][1] = 0;
    res[1][0] = 0;
    res[1][1] = 1;
    while(p > 0) {
        if(p & 1) {
            res = mat_mult(res, a);
        }
        a = mat_mult(a, a);
        p >>= 1;
    }

    return res;
}

int64_t n, m;

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

void solve() {
    // The capital (cell 0) connects to all n country towns, which form a cycle,
    // so every country town has degree 3. Collapse all country towns into one
    // "super" state and track only whether the current position is the capital
    // or a country town. From the capital there are n ways to step to a country
    // town; from a country town there is 1 way back to the capital and 2 ways
    // to an adjacent country town (its other two neighbours). This gives the
    // transition matrix [[0, 1], [n, 2]] acting on the (capital, town) state
    // vector. Raising it to the m-th power (binary exponentiation with the
    // vendored bigint, since the counts grow huge) and applying it to the
    // starting vector [1, 0] yields the number of length-m closed walks from
    // the capital as the first component.

    Mat a;
    a[0][0] = 0;
    a[0][1] = 1;
    a[1][0] = n;
    a[1][1] = 2;

    Mat am = mat_pow(a, m);
    cout << am[0][0] << '\n';
}

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

    read();
    solve();

    return 0;
}
```

5. Python Implementation with Detailed Comments
```python
import sys

def matrix_mult_2x2(A, B):
    """
    Multiply two 2×2 matrices A and B.
    A and B are lists of lists: [[a,b],[c,d]].
    Returns the 2×2 product.
    """
    return [
        [
            A[0][0]*B[0][0] + A[0][1]*B[1][0],
            A[0][0]*B[0][1] + A[0][1]*B[1][1]
        ],
        [
            A[1][0]*B[0][0] + A[1][1]*B[1][0],
            A[1][0]*B[0][1] + A[1][1]*B[1][1]
        ]
    ]

def matrix_pow_2x2(A, power):
    """
    Fast exponentiation of 2×2 matrix A to the integer power.
    Uses binary (square‐and‐multiply) method in O(log power) steps.
    """
    # Start with the identity matrix.
    result = [[1, 0], [0, 1]]
    while power > 0:
        if power & 1:
            result = matrix_mult_2x2(result, A)
        A = matrix_mult_2x2(A, A)
        power >>= 1
    return result

def solve(n, m):
    # Build transition matrix A = [[0,1],[n,2]]
    A = [[0, 1],
         [n, 2]]
    # Compute A^m
    Am = matrix_pow_2x2(A, m)
    # Initial vector v0 = [1, 0]^T → after m steps: v_m = Am · v0.
    # The number of ways to return to capital is v_m[0] = Am[0][0].
    return Am[0][0]

def main():
    data = sys.stdin.readline().split()
    if not data:
        return
    n, m = map(int, data)
    # Python ints are arbitrary precision by default.
    print(solve(n, m))

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