p407.ans1
======================
8

=================
p407.ans2
======================
21

=================
p407.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;
}

=================
p407.in1
======================
4 3

=================
p407.in2
======================
3 4

=================
p407.py
======================
# In C++ we would need long integers for this problem, so
# we will use Python's built-in int type which can handle large integers.
import sys


def matrix_mult_2x2(
    A_2x2: list[list[int]], B_2x2: list[list[int]]
) -> list[list[int]]:
    return [
        [
            A_2x2[0][0] * B_2x2[0][0] + A_2x2[0][1] * B_2x2[1][0],
            A_2x2[0][0] * B_2x2[0][1] + A_2x2[0][1] * B_2x2[1][1],
        ],
        [
            A_2x2[1][0] * B_2x2[0][0] + A_2x2[1][1] * B_2x2[1][0],
            A_2x2[1][0] * B_2x2[0][1] + A_2x2[1][1] * B_2x2[1][1],
        ],
    ]


def matrix_vector_mult_2x2(
    A_2x2: list[list[int]], v_2: list[int]
) -> list[int]:
    return [
        A_2x2[0][0] * v_2[0] + A_2x2[0][1] * v_2[1],
        A_2x2[1][0] * v_2[0] + A_2x2[1][1] * v_2[1],
    ]


def matrix_pow_2x2(A_2x2: list[list[int]], p: int) -> list[list[int]]:
    result = [[1, 0], [0, 1]]
    while p:
        if p % 2 == 1:
            result = matrix_mult_2x2(result, A_2x2)
        A_2x2 = matrix_mult_2x2(A_2x2, A_2x2)
        p //= 2

    return result


def solve(n: int, m: int) -> int:
    sys.set_int_max_str_digits(10000)

    # The main idea is to unite all of the non-start cells into a single
    # "super" cell. This is to compress the graph. You can notice that there
    # are only 4 options:
    # 1. You go from 0 to one of the n other cells. There are n options to do this,
    #    and this moves us from 0 to the super cell.
    # 2. You go from the super cell to 0. This is the opposite of option 1. Note that
    #    we have united the cells, so there is only one way to do this, not n (otherwise
    #    we overcount).
    # 3. We move from the super cell, to another cell in the super cell. Every cell has a 
    #    degree equal to 3, so there are exactly 2 ways to do this move, if we exclude the
    #    above option 2.
    # 4. We move from cell 0 to cell 0. This is impossible, since we are not allowed to
    #    stay in the same cell, so we do not count this option.

    A_2x2 = [[0, 1], [n, 2]]
    v_2 = [1, 0]
    Am_2x2 = matrix_pow_2x2(A_2x2, m)
    ans = matrix_vector_mult_2x2(Am_2x2, v_2)[0]

    return ans


def main():
    n, m = map(int, input().strip().split())
    result = solve(n, m)
    print(result)


if __name__ == "__main__":
    main()

=================
statement.txt
======================
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.

=================
