## 1. Abridged problem statement

You are given `n` linearly independent integer vectors `v1, ..., vn` in `n`-dimensional space. They form the parallelepiped

\[
P=\left\{\sum_{i=1}^{n} t_i v_i \mid 0 \le t_i \le 1\right\}.
\]

There is one document at every integer lattice point inside or on the boundary of `P`.

For a document point, let `k` be the dimension of the smallest facet of the parallelepiped containing it. Equivalently, if

\[
x=\sum t_i v_i,
\]

then `k` is the number of parameters `t_i` strictly between `0` and `1`. The document price is `2^k`.

Compute the sum of all document prices modulo a given prime `p`.

---

## 2. Detailed editorial

Let the parallelepiped be generated by the independent integer vectors

\[
v_1,\dots,v_n.
\]

Every point inside the parallelepiped has a unique representation

\[
x = \sum_{i=1}^{n} t_i v_i,\qquad 0 \le t_i \le 1.
\]

Because the vectors are linearly independent, the parameters `t_i` are unique.

For such a point:

- if `t_i = 0` or `t_i = 1`, the point lies on a boundary wall in direction `i`;
- if `0 < t_i < 1`, the point is interior in direction `i`.

Therefore the dimension `k` of the smallest facet containing the point is

\[
k = \#\{i : 0 < t_i < 1\}.
\]

The price of the document is then

\[
2^k.
\]

So the direct problem is:

\[
\sum_{x \in P \cap \mathbb Z^n} 2^{\#\{i:0<t_i<1\}}.
\]

---

### Key idea: half-open fundamental parallelepiped

Define the half-open parallelepiped

\[
H = \left\{\sum_{i=1}^{n} t_i v_i \mid 0 \le t_i < 1\right\}.
\]

The integer vectors `v_i` generate a sublattice

\[
L = \left\{\sum_{i=1}^{n} z_i v_i \mid z_i \in \mathbb Z\right\}.
\]

Because the `v_i` are integer vectors, `L` is a sublattice of `\mathbb Z^n`.

The half-open parallelepiped `H` is a fundamental domain of this lattice. That means every point in space can be uniquely shifted by an element of `L` into `H`.

In particular, the integer lattice points in `H` represent the cosets of

\[
\mathbb Z^n / L.
\]

The number of such cosets is the lattice index

\[
[\mathbb Z^n : L] = |\det V|,
\]

where `V` is the matrix formed by the vectors `v_i`.

Thus:

\[
|H \cap \mathbb Z^n| = |\det V|.
\]

---

### Grouping points of the closed parallelepiped

Now consider the original closed parallelepiped

\[
P = \left\{\sum t_i v_i \mid 0 \le t_i \le 1\right\}.
\]

Every integer point in `P` can be reduced to a unique integer point in `H` by replacing every coordinate `t_i = 1` with `t_i = 0`, i.e. subtracting `v_i`.

So we group points of `P` by their representative in `H`.

Take one representative

\[
y = \sum s_i v_i,\qquad 0 \le s_i < 1.
\]

Suppose exactly `z` of the coordinates `s_i` are zero.

For each coordinate where `s_i = 0`, we have two possibilities in the closed parallelepiped:

- keep `t_i = 0`;
- change it to `t_i = 1`.

For coordinates with `0 < s_i < 1`, we cannot change them, because adding `1` would exceed the allowed range.

Therefore this group contains exactly

\[
2^z
\]

points.

For every point in this group:

- the `z` changed/free zero-coordinates are always on a wall, either `0` or `1`;
- the remaining `n - z` coordinates are strictly inside.

So every point in the group has price

\[
2^{n-z}.
\]

Hence the total contribution of this whole group is

\[
2^z \cdot 2^{n-z} = 2^n.
\]

This is the crucial observation: every representative in `H` contributes exactly `2^n`.

There are `|\det V|` representatives, so the total answer is

\[
2^n \cdot |\det V|.
\]

We only need this modulo `p`.

---

### Computing the answer

We need

\[
2^n \cdot |\det V| \pmod p.
\]

The determinant can be enormous, so the C++ solution computes:

1. `det(V) mod p` using Gaussian elimination modulo prime `p`;
2. the sign of the real determinant using floating-point Gaussian elimination;
3. converts `det mod p` into `|det| mod p`.

If

\[
\det V \equiv 0 \pmod p,
\]

then the final answer is `0`, regardless of the sign.

Otherwise:

- if `det(V) > 0`, then `|det(V)| mod p = det mod p`;
- if `det(V) < 0`, then `|det(V)| mod p = -det mod p = p - det`.

Finally multiply by

\[
2^n \bmod p.
\]

Complexity:

\[
O(n^3)
\]

which is easily enough for `n <= 50`.

---

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

int n;
int64_t p;
vector<vector<int64_t>> a;

void read() {
    cin >> n >> p;
    a.assign(n, vector<int64_t>(n));
    cin >> a;
}

int64_t mul_mod(int64_t x, int64_t y) { return (int64_t)((__int128)x * y % p); }

int64_t pow_mod(int64_t base, int64_t e) {
    int64_t r = 1 % p;
    base %= p;
    while(e > 0) {
        if(e & 1) {
            r = mul_mod(r, base);
        }

        base = mul_mod(base, base);
        e >>= 1;
    }

    return r;
}

void solve() {
    // Every point of the body can be written x = sum t_i v_i for one choice
    // of coefficients t_1..t_n; since the vectors are independent those
    // coefficients (call them the parameters of x) are unique. The point is
    // in the body exactly when every parameter is in [0, 1]: all-zero is the
    // origin corner, all-one is the opposite corner, and a parameter equal to
    // 0 or 1 means x sits on one of the flat walls of that direction.
    //
    // A facet is one of those flat boundary pieces: a corner (dimension 0),
    // an edge (dimension 1), a 2D face, and so on up to the solid body itself
    // (dimension n). To name a facet, pin some parameters to a wall (0 or 1)
    // and let the rest range over [0, 1]; the points you get form a facet
    // whose dimension is the number of parameters left free. The smallest
    // facet through a given x pins exactly the parameters that already sit at
    // 0 or 1 and frees the rest, so its dimension is the count of strictly
    // interior parameters
    // 
    //     k = #{ i : 0 < t_i < 1 }.
    //     
    // A corner has k = 0, a point in the solid interior has k = n, and the
    // document price 2^k is precisely the 2^0 .. 2^n scale in the statement.
    // The task is to add up 2^k over every integer point (every document)
    // whose parameters all lie in [0, 1].
    //
    // Adding 2^k over the whole body head-on is hard, because the points and
    // their prices both vary. So we compare against a region whose integer
    // points we can count exactly: the half-open box
    // 
    //     H = { sum t_i v_i : 0 <= t_i < 1 },
    //
    // i.e. the body with its far walls (the t_i = 1 ones) shaved off. Take a
    // copy of H and shift it by an integer combination of v1..vn. These
    // shifted copies cover space with no gaps or overlaps -- they tile it,
    // the way unit squares tile the plane -- and this works precisely because
    // v1..vn are independent: being a basis they give every point of space a
    // unique parameter vector, and peeling off the integer parts of those
    // parameters (which is an integer combination of the v_i) slides the
    // point into exactly one copy of H. Unique copy means total coverage with
    // no overlap, so it really is a tiling.
    //
    // Now count the grid points in one copy. Each copy is an integer
    // translate of H, and translating by an integer vector is a bijection of
    // the grid Z^n onto itself, so every copy holds the same number m of grid
    // points. Count the grid points of a large box two ways: directly it is
    // about (volume of the box), since the grid has one point per unit
    // volume; and it is also (number of copies inside) * m, with about
    // (volume of the box) / |det V| copies because each copy has volume
    // |det V|. Equating the two forces m = |det V|, so H holds exactly
    // |det V| grid points. (This is the geometry-of-numbers fact that a
    // sublattice's fundamental domain contains [Z^n : L] = |det V| grid
    // points; its two-dimensional area form is Pick's theorem.) The plan is
    // to charge the entire price sum to these |det V| points of H.
    //
    // Claim: group the body's integer points by the point of H they reduce
    // to, and every group sums to exactly 2^n in price. There are |det V|
    // groups, so the whole price sum is 2^n * |det V| -- the answer. The
    // per-group total is always the same because prices factor coordinate by
    // coordinate, and the only freedom inside a group is which wall a
    // boundary coordinate sits on.
    //
    // Reduce a body point x to its representative y in H like this: while some
    // parameter t_i equals 1, subtract v_i, which lowers that parameter to 0,
    // keeps x in the body, and ends at a y with every parameter in [0, 1).
    // Conversely, the body points that reduce to a fixed y are exactly those
    // built from y by taking each coordinate where y is 0 and setting it to 0
    // or 1 (setting it to 1 adds v_i); a coordinate of y already in (0, 1)
    // cannot move, since raising it would pass 1. So if y has z zero
    // coordinates, its group is the 2^z independent 0/1 choices on them.
    //
    // Now total the group. Write the price 2^k as one factor per coordinate:
    // 2 if that coordinate is strictly interior, 1 if it sits on a wall (0 or
    // 1). Summed over the group each coordinate contributes a factor of 2 -- a
    // coordinate with y in (0, 1) is interior in every member (factor 2), and
    // a coordinate with y = 0 runs over both walls 0 and 1 (factors 1 + 1 =
    // 2). The group therefore totals 2^n for every y. Concretely, the group
    // has 2^z members each of price 2^(n - z), and 2^z * 2^(n - z) = 2^n.
    //
    // With |det V| groups each worth 2^n, the answer is 2^n * |det V| mod p.
    //
    // Gaussian elimination over F_p gives det V mod p but only up to sign,
    // and the entries are too large to take an exact big-integer det, so the
    // true sign is read off a plain floating-point LU with partial pivoting
    // and used to turn the residue into the nonnegative |det V|. When det V
    // is a multiple of p the residue is already 0 and the sign does not
    // matter.

    vector<vector<int64_t>> mat = a;
    int64_t det = 1 % p;
    for(int col = 0; col < n; col++) {
        int piv = -1;
        for(int r = col; r < n; r++) {
            if(mat[r][col] % p != 0) {
                piv = r;
                break;
            }
        }

        if(piv == -1) {
            det = 0;
            break;
        }

        if(piv != col) {
            swap(mat[piv], mat[col]);
            det = (p - det) % p;
        }

        det = mul_mod(det, mat[col][col] % p);
        int64_t inv = pow_mod(mat[col][col], p - 2);
        for(int r = col + 1; r < n; r++) {
            int64_t factor = mul_mod(mat[r][col] % p, inv);
            for(int c = col; c < n; c++) {
                mat[r][c] =
                    (mat[r][c] - mul_mod(factor, mat[col][c]) % p + p) % p;
            }
        }
    }

    if(det != 0) {
        vector<vector<long double>> b(n, vector<long double>(n));
        for(int i = 0; i < n; i++) {
            for(int j = 0; j < n; j++) {
                b[i][j] = (long double)a[i][j];
            }
        }

        int sign = 1;
        for(int col = 0; col < n; col++) {
            int piv = col;
            for(int r = col + 1; r < n; r++) {
                if(fabsl(b[r][col]) > fabsl(b[piv][col])) {
                    piv = r;
                }
            }

            if(piv != col) {
                swap(b[piv], b[col]);
                sign = -sign;
            }

            if(b[col][col] < 0) {
                sign = -sign;
            }

            for(int r = col + 1; r < n; r++) {
                long double factor = b[r][col] / b[col][col];
                for(int c = col; c < n; c++) {
                    b[r][c] -= factor * b[col][c];
                }
            }
        }

        if(sign < 0) {
            det = (p - det) % p;
        }
    }

    cout << mul_mod(det, pow_mod(2, n)) << '\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();
        solve();
    }

    return 0;
}

```

---

## 4. Python solution with detailed comments

This Python version uses the same formula:

\[
\text{answer} = 2^n \cdot |\det V| \pmod p.
\]

Unlike the C++ version, Python can conveniently compute the determinant exactly using the Bareiss algorithm, a fraction-free Gaussian elimination method.

```python
import sys


def bareiss_det(matrix):
    """
    Computes the exact integer determinant of a square matrix.

    Bareiss algorithm is a fraction-free version of Gaussian elimination.
    It keeps all intermediate values as integers and has O(n^3) complexity.

    For n <= 50, this is easily practical in Python.
    """

    # Make a deep copy, because the algorithm modifies the matrix.
    a = [row[:] for row in matrix]

    # Matrix size.
    n = len(a)

    # The determinant sign changes whenever we swap two rows.
    sign = 1

    # Previous pivot used for exact division.
    # In Bareiss algorithm, this starts as 1.
    prev = 1

    # Process columns 0 through n-2.
    # After that, the last diagonal value contains the determinant up to sign.
    for k in range(n - 1):
        # Find a nonzero pivot in column k at or below row k.
        pivot_row = k
        while pivot_row < n and a[pivot_row][k] == 0:
            pivot_row += 1

        # If no pivot exists, determinant is zero.
        if pivot_row == n:
            return 0

        # Swap pivot row into position k if needed.
        if pivot_row != k:
            a[k], a[pivot_row] = a[pivot_row], a[k]
            sign *= -1

        # Current pivot.
        pivot = a[k][k]

        # Update the remaining submatrix.
        for i in range(k + 1, n):
            for j in range(k + 1, n):
                # Bareiss update:
                #
                # a[i][j] = (a[i][j] * pivot - a[i][k] * a[k][j]) / prev
                #
                # The division is exact for integer matrices.
                a[i][j] = (a[i][j] * pivot - a[i][k] * a[k][j]) // prev

        # Entries below the pivot are no longer needed.
        for i in range(k + 1, n):
            a[i][k] = 0

        # Store this pivot for the next Bareiss division.
        prev = pivot

    # For a 1x1 matrix, determinant is just its only element.
    # The problem has n >= 2, but this keeps the function general.
    if n == 1:
        return a[0][0]

    # Last diagonal entry is determinant up to row-swap sign.
    return sign * a[n - 1][n - 1]


def main():
    # Read all input tokens.
    data = sys.stdin.buffer.read().split()

    # First token is n, second is p.
    n = int(data[0])
    p = int(data[1])

    # Read the n by n matrix.
    values = list(map(int, data[2:]))

    matrix = []
    index = 0

    for _ in range(n):
        row = values[index:index + n]
        matrix.append(row)
        index += n

    # Compute exact determinant.
    det = bareiss_det(matrix)

    # We need |det| modulo p.
    abs_det_mod = abs(det) % p

    # Compute 2^n modulo p.
    pow2 = pow(2, n, p)

    # Final formula.
    answer = abs_det_mod * pow2 % p

    # Print answer.
    print(answer)


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

---

## 5. Compressed editorial

Every point of the parallelepiped has a unique representation

\[
x=\sum t_i v_i,\qquad 0\le t_i\le 1.
\]

The smallest facet dimension containing `x` is the number of coordinates with

\[
0<t_i<1.
\]

So the document price is `2^k`.

Consider the half-open parallelepiped

\[
H=\left\{\sum t_i v_i \mid 0\le t_i<1\right\}.
\]

It is a fundamental domain of the lattice generated by the integer vectors `v_i`. Therefore the number of integer points in `H` is

\[
|\det V|.
\]

Now group every integer point of the closed parallelepiped by reducing coordinates equal to `1` back to `0`, landing in `H`.

For a representative in `H`, suppose `z` of its coordinates are zero. Then it corresponds to `2^z` points in the closed parallelepiped, obtained by choosing each zero coordinate to be either `0` or `1`.

Each such point has exactly `n-z` strictly interior coordinates, hence price

\[
2^{n-z}.
\]

So the whole group contributes

\[
2^z \cdot 2^{n-z}=2^n.
\]

There are `|\det V|` groups, hence

\[
\text{answer}=2^n|\det V|\pmod p.
\]

Compute the determinant modulo `p`, adjust its sign to get the absolute determinant modulo `p`, and multiply by `2^n mod p`.
