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

197. Nice Patterns Strike Back
time limit per test: 0.5 sec.
memory limit per test: 65536 KB
input: standard
output: standard



You might have noticed that there is the new fashion among rich people to have their yards tiled with black and white tiles, forming a pattern. The company Broken Tiles is well known as the best tiling company in our region. It provides the widest choices of nice patterns to tile your yard with. The pattern is nice if there is no square of size 2 × 2, such that all tiles in it have the same color. So patterns on the figure 1 are nice, while patterns on the figure 2 are not.





The president of the company wonders whether the variety of nice patterns he can provide to the clients is large enough. Thus he asks you to find out the number of nice patterns that can be used to tile the yard of size N × M. Now he is interested in the long term estimation, so he suggests N ≤ 10100. However, he does not like big numbers, so he asks you to find the answer modulo P.

Input

The input file contains three integer numbers: N (1 ≤ N ≤ 10100), M (1 ≤ M ≤ 5) and P (1 ≤ P ≤ 10000).

Output

Write the number of nice patterns of size N × M modulo P to the output file.

Sample test(s)

Input

Test #1

2 2 5

Test #2

3 3 23

Output

Test #1

4

Test #2

0


Author:	Andrew Stankevich
Resource:	Petrozavodsk Winter Trainings 2003
Date:	2003-02-06

<|response|>
1. Abridged Problem Statement
Given an N×M grid (1≤M≤5, N up to 10^100), color each cell black or white so that no 2×2 sub‐square is monochrome (all four tiles the same color). Compute the number of valid colorings modulo P (1≤P≤10 000).

2. Key Observations
- Any coloring of a row of width M can be encoded as an M‐bit mask (0=white, 1=black). There are S=2^M possible masks.
- Whether two consecutive rows (masks a and b) create a forbidden monochrome 2×2 square depends only on adjacent columns in these two rows. We can precompute a Boolean “valid(a,b)”.
- Let dp[i][mask] = number of ways to color rows 1..i with row i = mask. Then
  dp[i][cur] = Σ_{prev=0..S−1} dp[i−1][prev] × valid(prev,cur).
- This is a linear recurrence in vector form: V_i = T · V_{i−1}, where T is an S×S transition matrix.
- We need V_N = T^(N−1) · V_1. Since V_1 is all ones, the answer is simply the sum of all entries of T^(N−1).
- Since N can be up to 10^100, we keep N as a big integer and do fast exponentiation of T modulo P, reading the bits of the exponent N−1 via the big-integer mod-2 / divide-by-2 operations.

3. Full Solution Approach
a) State Encoding
   Each row is an integer mask in [0,2^M).
b) Build Transition Matrix T of size S×S, where
   T[a][b] = 1 if placing row b immediately below row a never forms a 2×2 monochrome block; otherwise 0.
   To check valid(a,b): for each column k=1..M−1, look at the four cells
     a₁ = (a>>(k−1))&1, a₂ = (a>>k)&1, b₁ = (b>>(k−1))&1, b₂ = (b>>k)&1
     and ensure they are not all 0 and not all 1.
c) Initial Vector V₁: for row 1 any mask is allowed, so V₁[mask]=1 for all mask.
d) Exponentiation
   We need T^(N−1) mod P. Store N as a big integer, form e = N−1, then run binary exponentiation:
   - R starts as the S×S identity matrix.
   - While e > 0:
       – If e is odd (e % 2 == 1), R = R × T mod P.
       – T = T × T mod P.
       – e = e / 2 (big-integer divide by 2).
e) The answer equals the sum of all entries of R = T^(N−1) (equivalently R·V₁ then summed), taken mod P.
f) Complexity
   - S = 2^M ≤ 32.
   - Matrix multiplication is O(S³) per multiply.
   - Exponentiation uses O(log N) ≈ O(330) squaring/multiplication steps.
   - Total ≈ 330 × 32³ ≈ 10⁷ basic ops, fits in 0.5 s.

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

bigint n;
int m, p;

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

vector<vector<int>> mat_mul(
    const vector<vector<int>>& a, const vector<vector<int>>& b
) {
    int sz = a.size();
    vector<vector<int>> c(sz, vector<int>(sz, 0));
    for(int i = 0; i < sz; i++) {
        for(int j = 0; j < sz; j++) {
            int64_t acc = 0;
            for(int k = 0; k < sz; k++) {
                acc += (int64_t)a[i][k] * b[k][j];
            }
            c[i][j] = acc % p;
        }
    }

    return c;
}

void solve() {
    // - A valid pattern has no monochromatic 2x2 square. Encode each row as a
    //   bitmask over M columns; a pair of consecutive rows is compatible iff for
    //   every adjacent column pair the four cells are not all equal, i.e. not all
    //   ones and not all zeros across the two rows.
    //
    // - Build the (2^M) x (2^M) transition matrix T where T[cur][next] = 1 if
    //   rows cur and next can be stacked. The number of valid N-row patterns is
    //   the sum of all entries of T^(N-1) times the all-ones start vector, all
    //   taken modulo P.
    //
    // - N is up to 10^100, so the exponent N-1 is handled as a big integer and
    //   binary exponentiation reads its bits via the modulo-2 / divide-by-2
    //   operations on the big integer type.

    int states = 1 << m;
    vector<vector<int>> trans(states, vector<int>(states, 0));
    for(int cur = 0; cur < states; cur++) {
        for(int nxt = 0; nxt < states; nxt++) {
            bool valid = true;
            for(int i = 1; i < m; i++) {
                bool all_one = ((cur >> (i - 1)) & 1) && ((cur >> i) & 1) &&
                               ((nxt >> (i - 1)) & 1) && ((nxt >> i) & 1);
                bool all_zero = !((cur >> (i - 1)) & 1) && !((cur >> i) & 1) &&
                                !((nxt >> (i - 1)) & 1) && !((nxt >> i) & 1);
                if(all_one || all_zero) {
                    valid = false;
                    break;
                }
            }

            if(valid) {
                trans[cur][nxt] = 1;
            }
        }
    }

    vector<vector<int>> result(states, vector<int>(states, 0));
    for(int i = 0; i < states; i++) {
        result[i][i] = 1 % p;
    }

    bigint e = n - bigint(1);
    while(e > bigint(0)) {
        if(e % 2) {
            result = mat_mul(result, trans);
        }

        trans = mat_mul(trans, trans);
        e /= 2;
    }

    int64_t ans = 0;
    for(int i = 0; i < states; i++) {
        for(int j = 0; j < states; j++) {
            ans += result[i][j];
        }
    }

    cout << ans % p << '\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

def matrix_multiply(A, B, mod):
    """Multiply two square matrices A and B under modulo mod."""
    n = len(A)
    C = [[0]*n for _ in range(n)]
    for i in range(n):
        for k in range(n):
            aik = A[i][k]
            if aik:
                for j in range(n):
                    C[i][j] = (C[i][j] + aik * B[k][j]) % mod
    return C

def matrix_vector_multiply(A, v, mod):
    """Multiply matrix A (n×n) by vector v (length n) under modulo mod."""
    n = len(A)
    res = [0]*n
    for i in range(n):
        s = 0
        for j in range(n):
            s += A[i][j] * v[j]
        res[i] = s % mod
    return res

def main():
    data = sys.stdin.read().split()
    N, M, P = int(data[0]), int(data[1]), int(data[2])

    S = 1 << M  # number of possible row masks

    # Build transition matrix T
    T = [[0]*S for _ in range(S)]
    for a in range(S):
        for b in range(S):
            ok = True
            for k in range(1, M):
                bitsum = ((a>>(k-1))&1) + ((a>>k)&1) + ((b>>(k-1))&1) + ((b>>k)&1)
                if bitsum == 0 or bitsum == 4:
                    ok = False
                    break
            if ok:
                T[a][b] = 1

    # Initial vector V1 (all ones)
    V = [1]*S

    # Fast exponentiation: R = T^(N-1) (Python ints are arbitrary precision)
    e = N - 1
    R = [[1 if i==j else 0 for j in range(S)] for i in range(S)]
    while e > 0:
        if e & 1:
            R = matrix_multiply(R, T, P)
        T = matrix_multiply(T, T, P)
        e >>= 1

    # Multiply R by V1 to get V_N
    VN = matrix_vector_multiply(R, V, P)

    # Sum up all entries modulo P
    answer = sum(VN) % P
    print(answer)

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