Chapter 5 Block Matrices: Exercises

import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=100)

# ── Visual hierarchy helpers ──────────────────────────────
def print_header(title):
    print("=" * 60)
    print(f"  {title}")
    print("=" * 60)

def print_step(step, desc):
    print(f"\n▶ Step {step}: {desc}")
    print("-" * 40)

def compare_print(label, actual, expected_desc):
    print(f"\n[{label}]")
    print(f"  Computed value: {actual}")
    print(f"  Expected value: {expected_desc}")

def check(label, actual, expected, tol=1e-9):
    ok = np.allclose(np.asarray(actual, dtype=complex), np.asarray(expected, dtype=complex), atol=tol)
    print(f"  [{label}] {'✓' if ok else '✗'}")
    if not ok:
        print("    actual  =", actual)
        print("    expected=", expected)
    return ok

print("Helper functions loaded")

# ============================================================
print_header("Exercise 1 | Ways to Partition a 3×3 Matrix")
# ============================================================
from itertools import combinations
from math import comb

A = np.arange(1, 10).reshape(3, 3)

def all_partitions(M, s, t):
    """List all s×t partitions of M; return [(horizontal cut points, vertical cut points, table of sub-blocks)]"""
    m, n = M.shape
    result = []
    for rc in combinations(range(1, m), s - 1):        # horizontal partition lines: between which rows they are drawn
        for cc in combinations(range(1, n), t - 1):    # vertical partition lines: between which columns they are drawn
            rb, cb = [0, *rc, m], [0, *cc, n]
            blocks = [[M[rb[i]:rb[i+1], cb[j]:cb[j+1]] for j in range(t)]
                      for i in range(s)]
            result.append((rc, cc, blocks))
    return result

# The sub-blocks written out in the solution (in increasing order of r and c; flattened in the order A_00, A_01, A_10, A_11)
expected = {
    (2, 2): [[[[1]], [[2, 3]], [[4], [7]], [[5, 6], [8, 9]]],      # r=1, c=1
             [[[1, 2]], [[3]], [[4, 5], [7, 8]], [[6], [9]]],      # r=1, c=2
             [[[1], [4]], [[2, 3], [5, 6]], [[7]], [[8, 9]]],      # r=2, c=1
             [[[1, 2], [4, 5]], [[3], [6]], [[7, 8]], [[9]]]],     # r=2, c=2
    (1, 2): [[[[1], [4], [7]], [[2, 3], [5, 6], [8, 9]]],          # c=1
             [[[1, 2], [4, 5], [7, 8]], [[3], [6], [9]]]],         # c=2
    (2, 1): [[[[1, 2, 3]], [[4, 5, 6], [7, 8, 9]]],                # r=1
             [[[1, 2, 3], [4, 5, 6]], [[7, 8, 9]]]],               # r=2
}

for step, (s, t) in enumerate([(2, 2), (1, 2), (2, 1)], start=1):
    print_step(step, f"{s}×{t} partitions")
    parts = all_partitions(A, s, t)
    for (rc, cc, B), exp in zip(parts, expected[(s, t)]):
        flat = [blk for row in B for blk in row]
        print(f"horizontal cut r={list(rc)}, vertical cut c={list(cc)}:",
              [f"{b.shape[0]}×{b.shape[1]}" for b in flat])
        same = all(np.array_equal(b, np.array(e)) for b, e in zip(flat, exp))
        check("sub-blocks agree with the solution", same, True)
        check("np.block reassembles A", np.array_equal(np.block(B), A), True)
    check(f"{s}×{t} partition count = C(2,{s-1})·C(2,{t-1})", len(parts), comb(2, s-1) * comb(2, t-1))

print_step(4, "Special sub-blocks: whole columns and whole rows")
check("1×2 partition c=1: A_00 = col_0(A)", A[:, :1].ravel(), A[:, 0])
check("1×2 partition c=2: A_01 = col_2(A)", A[:, 2:].ravel(), A[:, 2])
check("2×1 partition r=1: A_00 = row_0(A)", A[:1, :].ravel(), A[0, :])
check("2×1 partition r=2: A_10 = row_2(A)", A[2:, :].ravel(), A[2, :])

# ============================================================
print_header("Exercise 2 | Compatibility for Block Addition")
# ============================================================
P = np.array([[1, 0, 3], [0, 1, 4], [2, 5, 6]])
Q = np.array([[-1, 0, 0], [1, -1, 0], [1, 1, -1]])

def split(M, rcuts, ccuts):
    """Cut out the table of sub-blocks at the horizontal cut points rcuts and the vertical cut points ccuts"""
    rb, cb = [0, *rcuts, M.shape[0]], [0, *ccuts, M.shape[1]]
    return [[M[rb[i]:rb[i+1], cb[j]:cb[j+1]] for j in range(len(cb) - 1)]
            for i in range(len(rb) - 1)]

def same_structure(XB, YB):
    """Whether two tables of sub-blocks have the same block structure (same numbers of blocks, corresponding sub-blocks of the same size)"""
    return (len(XB) == len(YB) and len(XB[0]) == len(YB[0]) and
            all(x.shape == y.shape for xr, yr in zip(XB, YB) for x, y in zip(xr, yr)))

PB = split(P, [1], [1])       # row partition {0}|{1,2}, column partition {0}|{1,2}
QB_a = split(Q, [1], [1])     # (a)
QB_b = split(Q, [2], [2])     # (b): row partition {0,1}|{2}, column partition {0,1}|{2}

print_step(1, "Decide whether the block structures are the same")
check("(a) same block structure", same_structure(PB, QB_a), True)
check("(b) different block structures", same_structure(PB, QB_b), False)
print("  (b) size of P_00", PB[0][0].shape, ", size of Q_00", QB_b[0][0].shape)

print_step(2, "(a) Add block by block")
SB = [[PB[i][j] + QB_a[i][j] for j in range(2)] for i in range(2)]
check("P_00 + Q_00 = 0", SB[0][0], [[0]])
check("P_01 + Q_01 = [0 3]", SB[0][1], [[0, 3]])
check("P_10 + Q_10 = [1;3]", SB[1][0], [[1], [3]])
check("P_11 + Q_11 = [[0,4],[6,5]]", SB[1][1], [[0, 4], [6, 5]])
S = np.block(SB)
print(S)
check("reassembled = P + Q", S, P + Q)
check("P + Q = [[0,0,3],[1,0,4],[3,6,5]]", P + Q, [[0, 0, 3], [1, 0, 4], [3, 6, 5]])

# ============================================================
print_header("Exercise 3 | Sub-blocks of a Block Transpose")
# ============================================================
A = np.arange(1, 10).reshape(3, 3)

print_step(1, "Partition of A: row partition {0,1}|{2}, column partition {0}|{1,2}")
A00, A01 = A[:2, :1], A[:2, 1:]
A10, A11 = A[2:, :1], A[2:, 1:]
for name, blk in [("A_00", A00), ("A_01", A01), ("A_10", A10), ("A_11", A11)]:
    print(f"  {name} ({blk.shape[0]}×{blk.shape[1]}) = {blk.tolist()}")

print_step(2, "Partition of A^T: row partition {0}|{1,2}, column partition {0,1}|{2}")
AT = A.T
B00, B01 = AT[:1, :2], AT[:1, 2:]
B10, B11 = AT[1:, :2], AT[1:, 2:]
check("B_00 = (A_00)^T = [1 4]", B00, A00.T); check("B_00 values", B00, [[1, 4]])
check("B_01 = (A_10)^T = [7]", B01, A10.T);   check("B_01 values", B01, [[7]])
check("B_10 = (A_01)^T", B10, A01.T);         check("B_10 values", B10, [[2, 5], [3, 6]])
check("B_11 = (A_11)^T", B11, A11.T);         check("B_11 values", B11, [[8], [9]])

print_step(3, "Reassemble A^T")
check("np.block(B) = A^T = [[1,4,7],[2,5,8],[3,6,9]]", np.block([[B00, B01], [B10, B11]]),
      [[1, 4, 7], [2, 5, 8], [3, 6, 9]])
check("B_01 ≠ (A_01)^T (different sizes)", B01.shape == A01.T.shape, False)

# ============================================================
print_header("Exercise 4 | 4D Tensor Indices of Block Matrices")
# ============================================================

def to_tensor(M, p, q):
    """Block matrix → 4D tensor 𝒜[i,j,k,l] = [A_ij]_kl (sub-blocks of size p×q)"""
    m, n = M.shape[0] // p, M.shape[1] // q
    return M.reshape(m, p, n, q).transpose(0, 2, 1, 3)

A_a = np.array([[0, 1, 2, 3], [5, 7, 11, 13], [17, 19, 23, 29], [31, 37, 41, 43]])
A_b = np.array([[0, 0, 1, 2, 3, 4], [0, 5, 6, 7, 8, 9]])
A_c = np.array([[1, 2, 6, 7, 15, 16], [3, 5, 8, 14, 17, 27], [4, 9, 13, 18, 26, 31],
                [10, 12, 19, 25, 32, 42], [11, 20, 24, 33, 41, 50], [21, 23, 34, 40, 51, 61]])

# (part, matrix, shape of 𝒜, 𝒜[0,1,1,0], entry to find, index given in the solution)
cases = [("(a)", A_a, (2, 2, 2, 2), 11, 19, (1, 0, 0, 1)),
         ("(b)", A_b, (1, 3, 2, 2), 6, 4, (0, 2, 0, 1)),
         ("(c)", A_c, (3, 3, 2, 2), 8, 26, (1, 2, 0, 0))]

for step, (name, M, shape, v0110, val, idx) in enumerate(cases, start=1):
    print_step(step, f"{name} 𝒜[0,1,1,0] and the index of {val} in 𝒜")
    T = to_tensor(M, 2, 2)
    check(f"shape of 𝒜 = {shape}", T.shape, shape)
    # Check the definition entry by entry: 𝒜[i,j,k,l] = A[2i+k, 2j+l]
    ok = all(T[i, j, k, l] == M[2*i + k, 2*j + l] for i, j, k, l in np.ndindex(T.shape))
    check("𝒜[i,j,k,l] = A[2i+k, 2j+l]", ok, True)
    print("  A_01 =", T[0, 1].tolist())
    check(f"𝒜[0,1,1,0] = {v0110}", T[0, 1, 1, 0], v0110)
    found = tuple(int(x) for x in np.argwhere(T == val)[0])
    compare_print(f"Entry {val}: 4D index", found, f"{idx}")
    check(f"{val} = 𝒜{list(idx)}", found, idx)

# ============================================================
print_header("Exercise 5 | Partition Compatibility for Block Multiplication")
# ============================================================
X = np.array([[1, 0, 3], [0, 1, 4], [2, 5, 6]])
Y = np.array([[-1, 0, 0], [1, -1, 0], [1, 1, -1]])
XY = X @ Y

def split(M, rcuts, ccuts):
    """Cut out the table of sub-blocks at the horizontal cut points rcuts and the vertical cut points ccuts"""
    rb, cb = [0, *rcuts, M.shape[0]], [0, *ccuts, M.shape[1]]
    return [[M[rb[i]:rb[i+1], cb[j]:cb[j+1]] for j in range(len(cb) - 1)]
            for i in range(len(rb) - 1)]

def block_matmul(XB, YB):
    """Block multiplication (XY)_ij = Σ_k X_ik Y_kj; raises ValueError if the block levels or the sub-block sizes are incompatible"""
    if len(XB[0]) != len(YB):
        raise ValueError("incompatible at the block level")
    return [[sum(XB[i][k] @ YB[k][j] for k in range(len(YB)))
             for j in range(len(YB[0]))] for i in range(len(XB))]

# (part, (horizontal, vertical) cut points of X, (horizontal, vertical) cut points of Y, verdict in the solution)
items = [("(a)", ([1], [1]), ([1], [1]), True),
         ("(b)", ([1], [1]), ([2], [2]), False),
         ("(c)", ([], [1]),  ([1], [2]), True),
         ("(d)", ([2], []),  ([], [1]),  True),
         ("(e)", ([], [1]),  ([1], []),  True),
         ("(f)", ([], [2]),  ([], [2]),  False)]

print_step(1, "Compatibility: do the vertical cuts of X (column partition) equal the horizontal cuts of Y (row partition)?")
results = {}
for name, xc, yc, ok_expected in items:
    rule = (xc[1] == yc[0])
    try:
        results[name] = block_matmul(split(X, *xc), split(Y, *yc))
        computed = True
    except ValueError:
        computed = False
    print(f"  {name} X column-partition cut points {xc[1]}, Y row-partition cut points {yc[0]}")
    check(f"{name} verdict = {'Yes' if ok_expected else 'No'} (rule)", rule, ok_expected)
    check(f"{name} can actually be multiplied block by block", computed, ok_expected)

print_step(2, "(a) sub-blocks and intermediate terms")
XB, YB = split(X, [1], [1]), split(Y, [1], [1])
check("X_00 Y_00 = -1", XB[0][0] @ YB[0][0], [[-1]])
check("X_01 Y_10 = 3", XB[0][1] @ YB[1][0], [[3]])
check("X_01 Y_11 = [3 -3]", XB[0][1] @ YB[1][1], [[3, -3]])
check("X_10 Y_00 = [0;-2]", XB[1][0] @ YB[0][0], [[0], [-2]])
check("X_11 Y_10 = [5;11]", XB[1][1] @ YB[1][0], [[5], [11]])
check("X_11 Y_11 = [[3,-4],[1,-6]]", XB[1][1] @ YB[1][1], [[3, -4], [1, -6]])
C = results["(a)"]
check("(XY)_00 = 2", C[0][0], [[2]]); check("(XY)_01 = [3 -3]", C[0][1], [[3, -3]])
check("(XY)_10 = [5;9]", C[1][0], [[5], [9]]); check("(XY)_11", C[1][1], [[3, -4], [1, -6]])

print_step(3, "(c) sub-blocks and intermediate terms")
XB, YB = split(X, [], [1]), split(Y, [1], [2])
check("X_00 Y_00", XB[0][0] @ YB[0][0], [[-1, 0], [0, 0], [-2, 0]])
check("X_01 Y_10", XB[0][1] @ YB[1][0], [[3, 3], [5, 3], [11, 1]])
check("X_01 Y_11 = [-3;-4;-6]", XB[0][1] @ YB[1][1], [[-3], [-4], [-6]])
C = results["(c)"]
check("(XY)_00 = [[2,3],[5,3],[9,1]]", C[0][0], [[2, 3], [5, 3], [9, 1]])
check("(XY)_01 = [-3;-4;-6]", C[0][1], [[-3], [-4], [-6]])

print_step(4, "(d) sub-blocks")
C = results["(d)"]
check("(XY)_00 = [2;5]", C[0][0], [[2], [5]]); check("(XY)_01", C[0][1], [[3, -3], [3, -4]])
check("(XY)_10 = 9", C[1][0], [[9]]);          check("(XY)_11 = [1 -6]", C[1][1], [[1, -6]])

print_step(5, "(e) outer product + block-product contribution")
x, w = X[:, 0], Y[0, :]                      # x = col_0(X), w^T = row_0(Y)
outer = np.outer(x, w)                        # |x⟩⟨w|
second = X[:, 1:] @ Y[1:, :]
check("x w^T = [[-1,0,0],[0,0,0],[-2,0,0]]", outer, [[-1, 0, 0], [0, 0, 0], [-2, 0, 0]])
check("rank of the outer product = 1", np.linalg.matrix_rank(outer), 1)
check("X_01 Y_10 = [[3,3,-3],[5,3,-4],[11,1,-6]]", second, [[3, 3, -3], [5, 3, -4], [11, 1, -6]])
check("(e) single-block result = sum of the two terms", results["(e)"][0][0], outer + second)

print_step(6, "All four feasible expressions reassemble to the same XY")
print(XY)
check("XY = [[2,3,-3],[5,3,-4],[9,1,-6]]", XY, [[2, 3, -3], [5, 3, -4], [9, 1, -6]])
for name in ["(a)", "(c)", "(d)", "(e)"]:
    check(f"{name} np.block = XY", np.block(results[name]), XY)

# ============================================================
print_header("Exercise 6 | Outer-Product Representation of Rank-One Matrices")
# ============================================================
rng = np.random.default_rng(2026)

def rank_one_factor(A, tol=1e-10):
    """Construction from the proof: u = the first nonzero column, v_j = the coefficient of col_j(A) along u"""
    j0 = next(j for j in range(A.shape[1]) if np.linalg.norm(A[:, j]) > tol)
    u = A[:, j0]
    v = A.T @ u / (u @ u)          # v_j = ⟨u|col_j(A)⟩ / ⟨u|u⟩
    return u, v, j0

print_step(1, "Randomly generate rank-1 matrices (including cases with zero columns, and constructions not written as outer products)")
tests = []
for t in range(4):
    m, n = rng.integers(2, 7, size=2)
    u0, v0 = rng.normal(size=m), rng.normal(size=n)
    v0[rng.random(n) < 0.3] = 0.0              # make some columns zero
    if not v0.any():
        v0[-1] = 1.0
    tests.append(np.outer(u0, v0))
# Another construction: P E_00 Q (P, Q invertible; E_00 has a 1 only in entry (0,0)), not written directly as an outer product
for t in range(2):
    m, n = rng.integers(2, 7, size=2)
    P, Q = rng.normal(size=(m, m)), rng.normal(size=(n, n))
    E = np.zeros((m, n)); E[0, 0] = 1.0
    tests.append(P @ E @ Q)

print_step(2, "Verify each conclusion of the proof in turn")
for idx, A in enumerate(tests):
    print(f"\n  Instance {idx}: shape {A.shape}")
    check("rank(A) = 1", np.linalg.matrix_rank(A), 1)
    u, v, j0 = rank_one_factor(A)
    check("every column col_j(A) = v_j u", all(np.allclose(A[:, j], v[j] * u) for j in range(A.shape[1])), True)
    check(f"v_{{j0}} = v_{j0} = 1", v[j0], 1.0)
    check("u ≠ 0 and v ≠ 0", np.linalg.norm(u) > 0 and np.linalg.norm(v) > 0, True)
    check("A = u v^T = |u⟩⟨v|", np.outer(u, v), A)
    check("a_ij = u_i v_j", all(np.isclose(A[i, j], u[i] * v[j]) for i, j in np.ndindex(A.shape)), True)
    check("not unique: (2u)(v/2)^T = A", np.outer(2 * u, v / 2), A)

# ============================================================
print_header("Exercise 7 | Computations with Block Diagonal Matrices")
# ============================================================
from sympy import Matrix, Rational, diag as sdiag

print_step(1, "(a) M^2 = diag(M_0^2, M_1^2)")
M0 = np.array([[3, 2], [1, -1]]); M1 = np.array([[0, 1], [1, 0]])
Z = np.zeros((2, 2), dtype=int)
M = np.block([[M0, Z], [Z, M1]])
check("M_0^2 = [[11,4],[2,3]]", M0 @ M0, [[11, 4], [2, 3]])
check("M_1^2 = I", M1 @ M1, np.eye(2))
M2 = M @ M
print(M2)
check("M^2 = diag(M_0^2, M_1^2)", M2, np.block([[M0 @ M0, Z], [Z, M1 @ M1]]))
check("M^2 values", M2, [[11, 4, 0, 0], [2, 3, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])

print_step(2, "(b) Multiply the diagonal blocks block by block (block sizes 1,2,1)")
X = np.array([[2, 0, 0, 0], [0, -1, 2, 0], [0, 3, 2, 0], [0, 0, 0, -2]])
Y = np.array([[4, 0, 0, 0], [0, 2, 2, 0], [0, -1, -3, 0], [0, 0, 0, 7]])
X1, Y1 = X[1:3, 1:3], Y[1:3, 1:3]
check("X_1 Y_1 = [[-4,-8],[4,0]]", X1 @ Y1, [[-4, -8], [4, 0]])
check("corners: 2·4 = 8, (-2)·7 = -14", [X[0, 0] * Y[0, 0], X[3, 3] * Y[3, 3]], [8, -14])
XY = X @ Y
print(XY)
check("XY values", XY, [[8, 0, 0, 0], [0, -4, -8, 0], [0, 4, 0, 0], [0, 0, 0, -14]])

print_step(3, "(c) Invert block by block (sympy for exact fractions)")
N0 = Matrix([[3, 2], [5, 3]]); N1 = Matrix([[0, 2], [2, 0]])
check("det N_0 = -1", N0.det(), -1)
check("N_0^{-1} = [[-3,2],[5,-3]]", np.array(N0.inv(), dtype=float), [[-3, 2], [5, -3]])
check("N_1^{-1} = [[0,1/2],[1/2,0]]", np.array(N1.inv(), dtype=float), [[0, 0.5], [0.5, 0]])
N = sdiag(N0, N1)
Ninv = N.inv()
print(Ninv)
expected = Matrix([[-3, 2, 0, 0], [5, -3, 0, 0], [0, 0, 0, Rational(1, 2)], [0, 0, Rational(1, 2), 0]])
check("N^{-1} = diag(N_0^{-1}, N_1^{-1}) (exact)", Ninv == expected, True)
check("N N^{-1} = I_4", np.array(N * Ninv, dtype=float), np.eye(4))

# ============================================================
print_header("Exercise 8 | Inverse of a Block Upper Triangular Matrix")
# ============================================================
rng = np.random.default_rng(8)

def upper_block_inverse(A, B, D):
    """The formula of Exercise 8: [[A, B], [0, D]]^{-1}"""
    Ai, Di = np.linalg.inv(A), np.linalg.inv(D)
    Z = np.zeros((D.shape[0], A.shape[0]))
    return np.block([[Ai, -Ai @ B @ Di], [Z, Di]])

def schur_block_inverse(A, B, C, D):
    """Block inversion formula (thm-block-inverse-schur); requires A and S = D - C A^{-1} B to be invertible"""
    Ai = np.linalg.inv(A)
    Si = np.linalg.inv(D - C @ Ai @ B)
    return np.block([[Ai + Ai @ B @ Si @ C @ Ai, -Ai @ B @ Si],
                     [-Si @ C @ Ai, Si]])

# (b) deliberately include instances with p ≠ q
for k, (p, q) in enumerate([(2, 2), (2, 3), (3, 1), (1, 4)], start=1):
    print_step(k, f"Random instance: A is {p}×{p}, D is {q}×{q}, B is {p}×{q}")
    A = rng.standard_normal((p, p))
    B = rng.standard_normal((p, q))
    D = rng.standard_normal((q, q))
    Zqp = np.zeros((q, p))
    M = np.block([[A, B], [Zqp, D]])
    N = upper_block_inverse(A, B, D)
    check("M N = I", M @ N, np.eye(p + q))
    check("N M = I", N @ M, np.eye(p + q))
    check("N equals np.linalg.inv(M)", N, np.linalg.inv(M))
    # (c) When C = 0 the Schur complement is S = D, and the block inversion formula gives the same result
    check("when C = 0, S = D - C A^{-1} B = D", D - Zqp @ np.linalg.inv(A) @ B, D)
    check("agrees with the block inversion formula", schur_block_inverse(A, B, Zqp, D), N)

# ============================================================
print_header("Exercise 9 | Inverse of a Block Lower Triangular Matrix")
# ============================================================
rng = np.random.default_rng(9)

def lower_block_inverse(A, C, D):
    """The formula of Exercise 9: [[A, 0], [C, D]]^{-1}"""
    Ai, Di = np.linalg.inv(A), np.linalg.inv(D)
    Z = np.zeros((A.shape[0], D.shape[0]))
    return np.block([[Ai, Z], [-Di @ C @ Ai, Di]])

for k, (p, q) in enumerate([(2, 2), (3, 2), (1, 3)], start=1):
    print_step(k, f"Random instance: A is {p}×{p}, D is {q}×{q}, C is {q}×{p}")
    A = rng.standard_normal((p, p))
    C = rng.standard_normal((q, p))
    D = rng.standard_normal((q, q))
    Zpq = np.zeros((p, q))
    M = np.block([[A, Zpq], [C, D]])
    N = lower_block_inverse(A, C, D)
    check("N equals np.linalg.inv(M)", N, np.linalg.inv(M))
    # The transpose route: ((M^T)^{-1})^T, with (M^T)^{-1} from the formula of Exercise 8
    via_T = upper_block_inverse(A.T, C.T, D.T).T
    check("the transpose route ((M^T)^{-1})^T gives the same result", via_T, N)
    check("agrees with the block inversion formula (B = 0)", schur_block_inverse(A, Zpq, C, D), N)

# ============================================================
print_header("Exercise 10 | Numerical Block Inversion")
# ============================================================
from fractions import Fraction

def frac_matrix(rows):
    """Convert integers or 'a/b' strings into an array of Fraction objects (exact rational arithmetic)"""
    return np.array([[Fraction(x) for x in r] for r in rows], dtype=object)

def frac_inverse(M):
    """Exact inverse by Gauss–Jordan elimination (Fraction arithmetic)"""
    n = M.shape[0]
    aug = np.concatenate([M.copy(), frac_matrix(np.eye(n, dtype=int))], axis=1)
    for c in range(n):
        p = next(r for r in range(c, n) if aug[r, c] != 0)  # the row containing the pivot of column c
        aug[[c, p]] = aug[[p, c]]                             # swap rows c and p
        aug[c] = aug[c] / aug[c, c]
        for r in range(n):
            if r != c:
                aug[r] = aug[r] - aug[r, c] * aug[c]
    return aug[:, n:]

def frac_show(name, M):
    """Print a matrix laid out with fractions"""
    cells = [[str(x) for x in r] for r in M]
    w = max(len(s) for r in cells for s in r)
    print(f"{name} =")
    for r in cells:
        print("  [ " + "  ".join(s.rjust(w) for s in r) + " ]")

def check_exact(label, actual, expected):
    """Exact (Fraction) entry-by-entry comparison"""
    ok = np.array_equal(np.asarray(actual, dtype=object), np.asarray(expected, dtype=object))
    print(f"  [{label}] {'✓' if ok else '✗'} (exact)")
    return ok

# ---------- First matrix (1+2 partition, block upper triangular) ----------
print_step(1, "3×3 matrix: A = (2), B = (1 5), D = [[1,2],[2,3]]")
M1 = frac_matrix([[2, 1, 5], [0, 1, 2], [0, 2, 3]])
A, B, Z, D = M1[:1, :1], M1[:1, 1:], M1[1:, :1], M1[1:, 1:]
check_exact("lower-left block is zero", Z, frac_matrix([[0], [0]]))
Ai, Di = frac_inverse(A), frac_inverse(D)
check_exact("A^{-1} = (1/2)", Ai, frac_matrix([['1/2']]))
check_exact("D^{-1} = [[-3,2],[2,-1]]", Di, frac_matrix([[-3, 2], [2, -1]]))
Y = -Ai @ B @ Di
check_exact("upper-right block -A^{-1} B D^{-1} = (-7/2, 3/2)", Y, frac_matrix([['-7/2', '3/2']]))
M1_inv = np.block([[Ai, Y], [frac_matrix([[0], [0]]), Di]])
frac_show("M1^{-1}", M1_inv)
check_exact("agrees with the solution", M1_inv, frac_matrix([['1/2', '-7/2', '3/2'], [0, -3, 2], [0, 2, -1]]))
check_exact("M1 M1^{-1} = I", M1 @ M1_inv, frac_matrix(np.eye(3, dtype=int)))
check_exact("M1^{-1} M1 = I", M1_inv @ M1, frac_matrix(np.eye(3, dtype=int)))

# ---------- Second matrix (2+2 partition, block inversion formula) ----------
print_step(2, "4×4 matrix: 2+2 partition and the Schur complement")
M2 = frac_matrix([[3, 2, 0, 1], [5, 3, 1, 2], [0, 1, 0, 2], [1, 2, 2, 0]])
A, B, C, D = M2[:2, :2], M2[:2, 2:], M2[2:, :2], M2[2:, 2:]
Ai = frac_inverse(A)
check_exact("A^{-1} = [[-3,2],[5,-3]]", Ai, frac_matrix([[-3, 2], [5, -3]]))
check_exact("C A^{-1}", C @ Ai, frac_matrix([[5, -3], [7, -4]]))
check_exact("C A^{-1} B", C @ Ai @ B, frac_matrix([[-3, -1], [-4, -1]]))
S = D - C @ Ai @ B
check_exact("S = [[3,3],[6,1]]", S, frac_matrix([[3, 3], [6, 1]]))
Si = frac_inverse(S)
check_exact("S^{-1}", Si, frac_matrix([['-1/15', '1/5'], ['2/5', '-1/5']]))

print_step(3, "Substitute block by block into the block inversion formula")
check_exact("A^{-1} B", Ai @ B, frac_matrix([[2, 1], [-3, -1]]))
TR = -Ai @ B @ Si
BL = -Si @ C @ Ai
corr = Ai @ B @ Si @ C @ Ai
TL = Ai + corr
check_exact("upper-right block -A^{-1} B S^{-1}", TR, frac_matrix([['-4/15', '-1/5'], ['1/5', '2/5']]))
check_exact("lower-left block -S^{-1} C A^{-1}", BL, frac_matrix([['-16/15', '3/5'], ['-3/5', '2/5']]))
check_exact("A^{-1} B S^{-1} C A^{-1}", corr, frac_matrix([['41/15', '-8/5'], ['-19/5', '11/5']]))
check_exact("upper-left block", TL, frac_matrix([['-4/15', '2/5'], ['6/5', '-4/5']]))
M2_inv = np.block([[TL, TR], [BL, Si]])
frac_show("M2^{-1}", M2_inv)
check_exact("agrees with the exact Gauss–Jordan inverse", M2_inv, frac_inverse(M2))
check_exact("M2 M2^{-1} = I", M2 @ M2_inv, frac_matrix(np.eye(4, dtype=int)))
check_exact("M2^{-1} M2 = I", M2_inv @ M2, frac_matrix(np.eye(4, dtype=int)))
check("agrees with np.linalg.inv (floating point)", M2_inv.astype(float), np.linalg.inv(M2.astype(float)))

# ============================================================
print_header("Exercise 11 | Deriving the Block Inversion Formula by Block Gaussian Elimination")
# ============================================================
rng = np.random.default_rng(11)

for k, (p, q) in enumerate([(2, 2), (3, 2), (2, 4)], start=1):
    print_step(k, f"Random instance: A is {p}×{p}, D is {q}×{q}")
    A = rng.standard_normal((p, p)); B = rng.standard_normal((p, q))
    C = rng.standard_normal((q, p)); D = rng.standard_normal((q, q))
    M = np.block([[A, B], [C, D]])
    Ai = np.linalg.inv(A)
    S = D - C @ Ai @ B
    Ip, Iq = np.eye(p), np.eye(q)
    E = np.block([[Ip, np.zeros((p, q))], [-C @ Ai, Iq]])
    U = np.block([[A, B], [np.zeros((q, p)), S]])
    # (a) the block Gaussian elimination identity
    check("(a) E M = U", E @ M, U)
    # (b) block forms of E^{-1} and U^{-1}
    E_inv = np.block([[Ip, np.zeros((p, q))], [C @ Ai, Iq]])
    check("E^{-1} = [[I, 0], [C A^{-1}, I]]", E_inv, np.linalg.inv(E))
    check("U^{-1} (formula of Exercise 8)", upper_block_inverse(A, B, S), np.linalg.inv(U))
    # M^{-1} = U^{-1} E and the block inversion formula
    check("M^{-1} = U^{-1} E", upper_block_inverse(A, B, S) @ E, np.linalg.inv(M))
    check("U^{-1} E equals the block inversion formula", upper_block_inverse(A, B, S) @ E, schur_block_inverse(A, B, C, D))

# ============================================================
print_header("Exercise 12 | Inverting a 4×4 Matrix with the Block Inversion Formula")
# ============================================================
M = frac_matrix([[1, 1, 0, 1], [1, 2, 1, 0], [0, 1, 0, 0], [1, 0, 0, 0]])
A, B, C, D = M[:2, :2], M[:2, 2:], M[2:, :2], M[2:, 2:]
P = frac_matrix([[0, 1], [1, 0]])
I2 = frac_matrix(np.eye(2, dtype=int))
O2 = frac_matrix(np.zeros((2, 2), dtype=int))

print_step(1, "Partition: B = C = P, D = 0")
check_exact("B = P", B, P)
check_exact("C = P", C, P)
check_exact("D = 0", D, O2)

print_step(2, "A^{-1} and the Schur complement")
Ai = frac_inverse(A)
check_exact("A^{-1} = [[2,-1],[-1,1]]", Ai, frac_matrix([[2, -1], [-1, 1]]))
check_exact("C A^{-1} (swap the two rows)", C @ Ai, frac_matrix([[-1, 1], [2, -1]]))
check_exact("C A^{-1} B (then swap the two columns)", C @ Ai @ B, frac_matrix([[1, -1], [-1, 2]]))
S = D - C @ Ai @ B
check_exact("S = [[-1,1],[1,-2]]", S, frac_matrix([[-1, 1], [1, -2]]))
Si = frac_inverse(S)
check_exact("S^{-1} = [[-2,-1],[-1,-1]]", Si, frac_matrix([[-2, -1], [-1, -1]]))

print_step(3, "Substitute block by block into the formula")
check_exact("A^{-1} B", Ai @ B, frac_matrix([[-1, 2], [1, -1]]))
check_exact("A^{-1} B S^{-1}", Ai @ B @ Si, frac_matrix([[0, -1], [-1, 0]]))
check_exact("S^{-1} C A^{-1}", Si @ C @ Ai, frac_matrix([[0, -1], [-1, 0]]))
check_exact("A^{-1} B S^{-1} C A^{-1}", Ai @ B @ Si @ C @ Ai, frac_matrix([[-2, 1], [1, -1]]))
TL = Ai + Ai @ B @ Si @ C @ Ai
TR = -Ai @ B @ Si
BL = -Si @ C @ Ai
check_exact("upper-left block = 0", TL, O2)
check_exact("upper-right block = P", TR, P)
check_exact("lower-left block = P", BL, P)

print_step(4, "Assemble and verify")
M_inv = np.block([[TL, TR], [BL, Si]])
frac_show("M^{-1}", M_inv)
check_exact("agrees with the solution", M_inv, frac_matrix([[0, 0, 0, 1], [0, 0, 1, 0], [0, 1, -2, -1], [1, 0, -1, -1]]))
check_exact("upper-right block check A P + P S^{-1} = 0", A @ P + P @ Si, O2)
check_exact("M M^{-1} = I_4", M @ M_inv, frac_matrix(np.eye(4, dtype=int)))
check_exact("M^{-1} M = I_4", M_inv @ M, frac_matrix(np.eye(4, dtype=int)))

# ============================================================
print_header("Exercise 13 | The Block Inversion Formula with D as Pivot")
# ============================================================
rng = np.random.default_rng(13)

def schur_D_block_inverse(A, B, C, D):
    """The formula of Exercise 13: requires D and T = A - B D^{-1} C to be invertible (A may or may not be invertible)"""
    Di = np.linalg.inv(D)
    Ti = np.linalg.inv(A - B @ Di @ C)
    return np.block([[Ti, -Ti @ B @ Di],
                     [-Di @ C @ Ti, Di + Di @ C @ Ti @ B @ Di]])

for k, (p, q) in enumerate([(2, 2), (3, 2), (2, 4), (4, 3)], start=1):
    print_step(k, f"A is a non-invertible {p}×{p} matrix, D is {q}×{q}")
    # construct A of rank p-1 (not invertible)
    A = rng.standard_normal((p, p - 1)) @ rng.standard_normal((p - 1, p))
    B = rng.standard_normal((p, q)); C = rng.standard_normal((q, p)); D = rng.standard_normal((q, q))
    M = np.block([[A, B], [C, D]])
    print(f"  rank(A) = {np.linalg.matrix_rank(A)} < {p}, so A is not invertible")
    check("A is not invertible (rank(A) = p-1)", np.linalg.matrix_rank(A), p - 1)
    N = schur_D_block_inverse(A, B, C, D)
    check("formula = np.linalg.inv(M)", N, np.linalg.inv(M))
    # "Reading it off directly": Π M Π^T = [[D, C], [B, A]]
    Pi = np.block([[np.zeros((q, p)), np.eye(q)], [np.eye(p), np.zeros((p, q))]])
    Mp = Pi @ M @ Pi.T
    check("Π M Π^T = [[D, C], [B, A]]", Mp, np.block([[D, C], [B, A]]))
    check("M^{-1} = Π^T (M')^{-1} Π, with (M')^{-1} from the block inversion formula", Pi.T @ schur_block_inverse(D, C, B, A) @ Pi, N)

print_step(5, "Example: A = 0, B = C = D = I_n (n = 3)")
n = 3
I, O = np.eye(n), np.zeros((n, n))
M = np.block([[O, I], [I, I]])
check("M^{-1} = [[-I, I], [I, 0]]", schur_D_block_inverse(O, I, I, I), np.block([[-I, I], [I, O]]))
check("M M^{-1} = I", M @ np.block([[-I, I], [I, O]]), np.eye(2 * n))

print_step(6, "Note: when A, D, S, T are all invertible the two formulas agree (Woodbury form)")
p, q = 3, 2
A = rng.standard_normal((p, p)); B = rng.standard_normal((p, q))
C = rng.standard_normal((q, p)); D = rng.standard_normal((q, q))
Ai = np.linalg.inv(A)
lhs = np.linalg.inv(A - B @ np.linalg.inv(D) @ C)
rhs = Ai + Ai @ B @ np.linalg.inv(D - C @ Ai @ B) @ C @ Ai
check("(A - B D^{-1} C)^{-1} = A^{-1} + A^{-1} B S^{-1} C A^{-1}", lhs, rhs)
check("the two block formulas give the same M^{-1}", schur_D_block_inverse(A, B, C, D), schur_block_inverse(A, B, C, D))

# ============================================================
print_header("Exercise 14 | Must M Be Non-invertible When Neither Diagonal Block Is Invertible?")
# ============================================================
print_step(1, "Smallest counterexample: A = D = (0), B = C = (1)")
M = np.array([[0., 1.], [1., 0.]])
check("A = (0) is not invertible (rank 0)", np.linalg.matrix_rank(M[:1, :1]), 0)
check("M M = I_2, so M^{-1} = M", M @ M, np.eye(2))
for n in [2, 3]:
    I, O = np.eye(n), np.zeros((n, n))
    Mn = np.block([[O, I], [I, O]])
    check(f"n = {n}: [[0, I], [I, 0]]^2 = I", Mn @ Mn, np.eye(2 * n))

print_step(2, "Counterexample with nonzero diagonal blocks: M = [[J_2, I_2], [I_2, J_2]] = J_4 - P")
J2, I2 = np.ones((2, 2)), np.eye(2)
M = np.block([[J2, I2], [I2, J2]])
J4, P = np.ones((4, 4)), np.fliplr(np.eye(4))   # P: the antidiagonal permutation matrix
print("M =\n", M)
check("M = J_4 - P", M, J4 - P)
check("rank(J_2) = 1, so the diagonal blocks are not invertible", np.linalg.matrix_rank(J2), 1)
check("J_4^2 = 4 J_4", J4 @ J4, 4 * J4)
check("J_4 P = P J_4 = J_4", np.stack([J4 @ P, P @ J4]), np.stack([J4, J4]))
check("P^2 = I_4", P @ P, np.eye(4))
M_inv = J4 / 3 - P
check("M^{-1} = J_4/3 - P", M_inv, np.linalg.inv(M))
check("agrees with the solution matrix (1/3)[[1,1,1,-2],...]", M_inv,
      np.array([[1, 1, 1, -2], [1, 1, -2, 1], [1, -2, 1, 1], [-2, 1, 1, 1]]) / 3)
check("M M^{-1} = I_4", M @ M_inv, np.eye(4))

print_step(3, "Comparison: block upper triangular with A not invertible forces M to be non-invertible; A, D invertible does not make M invertible")
rng = np.random.default_rng(14)
A = np.array([[1., 2.], [2., 4.]])            # not invertible
B, D = rng.standard_normal((2, 2)), rng.standard_normal((2, 2))
Mt = np.block([[A, B], [np.zeros((2, 2)), D]])
v = np.array([2., -1.])                         # A v = 0
check("A v = 0", A @ v, np.zeros(2))
check("M [v; 0] = 0, so M is not invertible", Mt @ np.concatenate([v, np.zeros(2)]), np.zeros(4))
check("[[1,1],[1,1]]: A = D = (1) invertible but rank(M) = 1", np.linalg.matrix_rank(np.ones((2, 2))), 1)

# ============================================================
print_header("Exercise 15 | Choosing the Pivot to Invert a 4×4 Matrix")
# ============================================================
I4 = frac_matrix(np.eye(4, dtype=int))

# ---------- M1: A is not invertible, so use the Schur complement T of D ----------
print_step(1, "M1: A is not invertible, D is invertible")
M1 = frac_matrix([[1, 1, 0, 1], [1, 1, 0, 0], [0, 1, 1, 1], [0, 0, 1, 2]])
A, B, C, D = M1[:2, :2], M1[:2, 2:], M1[2:, :2], M1[2:, 2:]
check("rank(A) = 1, so A is not invertible", np.linalg.matrix_rank(A.astype(float)), 1)
Di = frac_inverse(D)
check_exact("D^{-1} = [[2,-1],[-1,1]]", Di, frac_matrix([[2, -1], [-1, 1]]))

print_step(2, "M1: Schur complement T = A - B D^{-1} C")
check_exact("B D^{-1}", B @ Di, frac_matrix([[-1, 1], [0, 0]]))
check_exact("B D^{-1} C", B @ Di @ C, frac_matrix([[0, -1], [0, 0]]))
T = A - B @ Di @ C
check_exact("T = [[1,2],[1,1]]", T, frac_matrix([[1, 2], [1, 1]]))
Ti = frac_inverse(T)
check_exact("T^{-1} = [[-1,2],[1,-1]]", Ti, frac_matrix([[-1, 2], [1, -1]]))

print_step(3, "M1: substitute block by block into the formula of Exercise 13")
TR = -Ti @ B @ Di
check_exact("upper-right block -T^{-1} B D^{-1}", TR, frac_matrix([[-1, 1], [1, -1]]))
check_exact("D^{-1} C", Di @ C, frac_matrix([[0, 2], [0, -1]]))
BL = -Di @ C @ Ti
check_exact("lower-left block -D^{-1} C T^{-1}", BL, frac_matrix([[-2, 2], [1, -1]]))
check_exact("D^{-1} C T^{-1} B D^{-1}", Di @ C @ Ti @ B @ Di, frac_matrix([[-2, 2], [1, -1]]))
BR = Di + Di @ C @ Ti @ B @ Di
check_exact("lower-right block = [[0,1],[0,0]]", BR, frac_matrix([[0, 1], [0, 0]]))
M1_inv = np.block([[Ti, TR], [BL, BR]])
frac_show("M1^{-1}", M1_inv)
check_exact("agrees with the solution", M1_inv, frac_matrix([[-1, 2, -1, 1], [1, -1, 1, -1], [-2, 2, 0, 1], [1, -1, 0, 0]]))
print("  inner products of row 0 with columns 0 and 1:", M1[0] @ M1_inv[:, 0], M1[0] @ M1_inv[:, 1])
check_exact("M1 M1^{-1} = I_4", M1 @ M1_inv, I4)
check_exact("M1^{-1} M1 = I_4", M1_inv @ M1, I4)

# ---------- M2: A is invertible, so use the block inversion formula ----------
print_step(4, "M2: A^{-1} and the Schur complement S")
M2 = frac_matrix([[1, 0, 0, 1], [0, 2, 1, 0], [0, 1, 1, 0], [1, 0, 0, 2]])
A, B, C, D = M2[:2, :2], M2[:2, 2:], M2[2:, :2], M2[2:, 2:]
Ai = frac_inverse(A)
check_exact("A^{-1} = diag(1, 1/2)", Ai, frac_matrix([[1, 0], [0, '1/2']]))
check_exact("C A^{-1}", C @ Ai, frac_matrix([[0, '1/2'], [1, 0]]))
check_exact("C A^{-1} B", C @ Ai @ B, frac_matrix([['1/2', 0], [0, 1]]))
S = D - C @ Ai @ B
check_exact("S = diag(1/2, 1)", S, frac_matrix([['1/2', 0], [0, 1]]))
Si = frac_inverse(S)
check_exact("S^{-1} = diag(2, 1)", Si, frac_matrix([[2, 0], [0, 1]]))

print_step(5, "M2: substitute block by block into the block inversion formula")
P = frac_matrix([[0, 1], [1, 0]])
check_exact("A^{-1} B", Ai @ B, frac_matrix([[0, 1], ['1/2', 0]]))
check_exact("A^{-1} B S^{-1} = P", Ai @ B @ Si, P)
check_exact("S^{-1} C A^{-1} = P", Si @ C @ Ai, P)
check_exact("A^{-1} B S^{-1} C A^{-1} = diag(1, 1/2)", Ai @ B @ Si @ C @ Ai, frac_matrix([[1, 0], [0, '1/2']]))
TL = Ai + Ai @ B @ Si @ C @ Ai
check_exact("upper-left block = diag(2, 1)", TL, frac_matrix([[2, 0], [0, 1]]))
M2_inv = np.block([[TL, -Ai @ B @ Si], [-Si @ C @ Ai, Si]])
frac_show("M2^{-1}", M2_inv)
check_exact("agrees with the solution", M2_inv, frac_matrix([[2, 0, 0, -1], [0, 1, -1, 0], [0, -1, 2, 0], [-1, 0, 0, 1]]))
check_exact("M2 M2^{-1} = I_4", M2 @ M2_inv, I4)
check_exact("M2^{-1} M2 = I_4", M2_inv @ M2, I4)

print_step(6, "M2: indices {0,3} and {1,2} are not coupled (block diagonal after reordering)")
i03, i12 = [0, 3], [1, 2]
check_exact("entries of M2 in {0,3}×{1,2} are all 0", M2[np.ix_(i03, i12)], frac_matrix([[0, 0], [0, 0]]))
check_exact("[[1,1],[1,2]]^{-1} sits at the {0,3} positions of M2^{-1}",
            frac_inverse(M2[np.ix_(i03, i03)]), M2_inv[np.ix_(i03, i03)])
check_exact("[[2,1],[1,1]]^{-1} sits at the {1,2} positions of M2^{-1}",
            frac_inverse(M2[np.ix_(i12, i12)]), M2_inv[np.ix_(i12, i12)])

# ============================================================
print_header("Exercise 16 | Inverse of the Block Matrix [[I, X], [Y, 0]]")
# ============================================================
rng = np.random.default_rng(16)

def ixy0_inverse(X, Y):
    """The formula of Exercise 16 (a); X is n×k, Y is k×n, and YX is invertible"""
    n, k = X.shape
    R = np.linalg.inv(Y @ X)
    return np.block([[np.eye(n) - X @ R @ Y, X @ R],
                     [R @ Y, -R]])

for step, n in enumerate([2, 3, 4], start=1):
    print_step(step, f"n = {n}: (a) formula, (b) two-sided check, (c) simplification")
    X = rng.standard_normal((n, n)); Y = rng.standard_normal((n, n))
    M = np.block([[np.eye(n), X], [Y, np.zeros((n, n))]])
    N = ixy0_inverse(X, Y)
    check("(a) formula = np.linalg.inv(M)", N, np.linalg.inv(M))
    check("(a) formula = block inversion formula (A = I, D = 0)", N, schur_block_inverse(np.eye(n), X, Y, np.zeros((n, n))))
    check("(b) M M^{-1} = I_{2n}", M @ N, np.eye(2 * n))
    check("(b) M^{-1} M = I_{2n}", N @ M, np.eye(2 * n))
    Xi, Yi = np.linalg.inv(X), np.linalg.inv(Y)
    simple = np.block([[np.zeros((n, n)), Yi], [Xi, -Xi @ Yi]])
    check("(c) M^{-1} = [[0, Y^{-1}], [X^{-1}, -X^{-1} Y^{-1}]]", simple, N)
    check("(c) upper-left block I - X (YX)^{-1} Y = 0", N[:n, :n], np.zeros((n, n)))

print_step(4, "Note: in the square case, YX invertible ⟹ X, Y invertible (checked via rank)")
n = 3
X = rng.standard_normal((n, n)); Y = rng.standard_normal((n, n))
check("when rank(YX) = n, rank(X) = rank(Y) = n",
      [np.linalg.matrix_rank(Y @ X), np.linalg.matrix_rank(X), np.linalg.matrix_rank(Y)], [n, n, n])

print_step(5, "Note: the rectangular case, X is n×k and Y is k×n (k < n)")
for n, k in [(4, 2), (5, 3)]:
    X = rng.standard_normal((n, k)); Y = rng.standard_normal((k, n))
    M = np.block([[np.eye(n), X], [Y, np.zeros((k, k))]])
    N = ixy0_inverse(X, Y)
    check(f"n = {n}, k = {k}: M M^{{-1}} = I_{{n+k}}", M @ N, np.eye(n + k))
    check(f"n = {n}, k = {k}: M^{{-1}} M = I_{{n+k}}", N @ M, np.eye(n + k))
    Pr = X @ np.linalg.inv(Y @ X) @ Y
    check(f"n = {n}, k = {k}: X (YX)^{{-1}} Y is a projection (P^2 = P)", Pr @ Pr, Pr)

# ============================================================
print_header("Exercise 17 | Inverse of the Block Matrix [[A, B], [-B, A]]")
# ============================================================
rng = np.random.default_rng(17)

for step, n in enumerate([2, 3, 4], start=1):
    print_step(step, f"(a) n = {n}: general (noncommuting) A, B")
    A = rng.standard_normal((n, n)); B = rng.standard_normal((n, n))
    print(f"  ||AB - BA|| = {np.linalg.norm(A @ B - B @ A):.4f} (noncommuting)")
    M = np.block([[A, B], [-B, A]])
    Mi = np.linalg.inv(M)
    Ai = np.linalg.inv(A)
    K = Ai @ B
    I = np.eye(n)
    S = A + B @ Ai @ B
    check("S = A (I + K^2)", S, A @ (I + K @ K))
    G = np.linalg.inv(I + K @ K)
    check("G K = K G", G @ K, K @ G)
    X, Y = G @ Ai, -K @ G @ Ai
    check("X = (A + B A^{-1} B)^{-1}", X, np.linalg.inv(S))
    check("M^{-1} = [[X, Y], [-Y, X]]", np.block([[X, Y], [-Y, X]]), Mi)
    check("M^{-1} upper-left block = lower-right block", Mi[:n, :n], Mi[n:, n:])
    check("M^{-1} lower-left block = -upper-right block", Mi[n:, :n], -Mi[:n, n:])
    J = np.block([[np.zeros((n, n)), I], [-I, np.zeros((n, n))]])
    check("Note: M J = J M and M^{-1} J = J M^{-1}", np.stack([M @ J, Mi @ J]), np.stack([J @ M, J @ Mi]))

print_step(4, "(b) commuting A, B: take B to be a polynomial in A")
for n in [2, 3, 4]:
    A = rng.standard_normal((n, n))
    B = A @ A - 2 * A + 0.5 * np.eye(n)        # B = A^2 - 2A + I/2, which commutes with A
    M = np.block([[A, B], [-B, A]])
    R = np.linalg.inv(A @ A + B @ B)
    formula = np.block([[A @ R, -B @ R], [B @ R, A @ R]])
    check(f"n = {n}: AB = BA", A @ B, B @ A)
    check(f"n = {n}: M^{{-1}} = [[A R, -B R], [B R, A R]]", formula, np.linalg.inv(M))
    check(f"n = {n}: R A = A R", R @ A, A @ R)

print_step(5, "The formula of (b) does not need A to be invertible: A = diag(0, 1), B = diag(1, 0)")
A = np.diag([0., 1.]); B = np.diag([1., 0.])
M = np.block([[A, B], [-B, A]])
R = np.linalg.inv(A @ A + B @ B)
check("A^2 + B^2 = I", A @ A + B @ B, np.eye(2))
check("M^{-1} = [[A R, -B R], [B R, A R]]", np.block([[A @ R, -B @ R], [B @ R, A @ R]]), np.linalg.inv(M))

print_step(6, "Note: n = 1 corresponds to the complex number a + bi (a = 3, b = 4)")
a, b = 3., 4.
Mc = np.array([[a, b], [-b, a]])
check("[[a,b],[-b,a]]^{-1} = (1/(a^2+b^2)) [[a,-b],[b,a]]", np.linalg.inv(Mc), np.array([[a, -b], [b, a]]) / (a**2 + b**2))
z_inv = 1 / complex(a, b)
check("corresponds to 1/(a+bi) = (a - bi)/(a^2+b^2)", [z_inv.real, z_inv.imag], [a / 25, -b / 25])

# ============================================================
print_header("Exercise 18 | Computing Tensor Products and Their Noncommutativity")
# ============================================================
I2 = np.eye(2, dtype=int)
M = np.array([[1, 2], [3, 4]])
U = np.array([[1, 2, 4], [0, 1, 8], [0, 0, 1]])
V = np.array([[1, -1], [0, 1]])
a = np.array([[1], [0], [0]])   # a column vector in R^3
b = np.array([[1], [0]])        # a column vector in R^2

# --- Part 1 ---
print_step(1, "Part 1: I_2 ⊗ M and M ⊗ I_2")
K1, K2 = np.kron(I2, M), np.kron(M, I2)
print("I_2 ⊗ M =\n", K1)
print("M ⊗ I_2 =\n", K2)
Z2 = np.zeros((2, 2), dtype=int)
check("I_2 ⊗ M = diag(M, M)", K1, np.block([[M, Z2], [Z2, M]]))
check("M ⊗ I_2 agrees with the hand computation", K2, [[1, 0, 2, 0], [0, 1, 0, 2], [3, 0, 4, 0], [0, 3, 0, 4]])
check("I_2 ⊗ M ≠ M ⊗ I_2", np.array_equal(K1, K2), False)
check("row 0: (1,2,0,0) and (1,0,2,0)", [K1[0, :], K2[0, :]], [[1, 2, 0, 0], [1, 0, 2, 0]])

# --- Part 2 ---
print_step(2, "Part 2: U ⊗ V and V ⊗ U")
K3, K4 = np.kron(U, V), np.kron(V, U)
print("U ⊗ V =\n", K3)
print("V ⊗ U =\n", K4)
check("U ⊗ V agrees with the hand computation", K3, [[1, -1, 2, -2, 4, -4],
                             [0, 1, 0, 2, 0, 4],
                             [0, 0, 1, -1, 8, -8],
                             [0, 0, 0, 1, 0, 8],
                             [0, 0, 0, 0, 1, -1],
                             [0, 0, 0, 0, 0, 1]])
check("V ⊗ U = [[U, -U], [0, U]]", K4, np.block([[U, -U], [np.zeros_like(U), U]]))
check("U ⊗ V ≠ V ⊗ U", np.array_equal(K3, K4), False)
check("(0,1) entries: -1 and 2", [K3[0, 1], K4[0, 1]], [-1, 2])
for name, K in [("U ⊗ V", K3), ("V ⊗ U", K4)]:
    check(f"{name} is unit upper triangular", [np.allclose(K, np.triu(K)), np.allclose(np.diag(K), 1)], [True, True])

# --- Part 3 ---
print_step(3, "Part 3: a ⊗ b and b ⊗ a")
ab, ba = np.kron(a, b), np.kron(b, a)
print("(a ⊗ b)^T =", ab.ravel(), "  (b ⊗ a)^T =", ba.ravel())
check("a ⊗ b = (1,0,0,0,0,0)^T", ab.ravel(), [1, 0, 0, 0, 0, 0])
check("a ⊗ b = b ⊗ a", ab, ba)
a2 = np.array([[0], [1], [0]])
print("Taking a' = (0,1,0)^T instead: (a'⊗b)^T =", np.kron(a2, b).ravel(), " (b⊗a')^T =", np.kron(b, a2).ravel())
check("a' ⊗ b = (0,0,1,0,0,0)^T", np.kron(a2, b).ravel(), [0, 0, 1, 0, 0, 0])
check("b ⊗ a' = (0,1,0,0,0,0)^T", np.kron(b, a2).ravel(), [0, 1, 0, 0, 0, 0])

# --- Note: the permutation matrix P ---
print_step(4, "Note: B ⊗ A = P (A ⊗ B) P^T")
def swap_perm(m, n):
    # P(x ⊗ y) = y ⊗ x: send index i*n + j to j*m + i
    P = np.zeros((m * n, m * n), dtype=int)
    for i in range(m):
        for j in range(n):
            P[j * m + i, i * n + j] = 1
    return P
P22, P32 = swap_perm(2, 2), swap_perm(3, 2)
check("M ⊗ I_2 = P (I_2 ⊗ M) P^T", K2, P22 @ K1 @ P22.T)
check("V ⊗ U = P (U ⊗ V) P^T", K4, P32 @ K3 @ P32.T)
rng = np.random.default_rng(18)
x, y = rng.standard_normal(3), rng.standard_normal(2)
check("P (x ⊗ y) = y ⊗ x (random vectors)", P32 @ np.kron(x, y), np.kron(y, x))

# ============================================================
print_header("Exercise 19 | Verifying the Transpose and Trace Properties")
# ============================================================
A = np.array([[1, 3], [0, -1]])
B = np.array([[2, 0], [-2, 3]])

print_step(1, "Compute A ⊗ B and (A ⊗ B)^T")
K = np.kron(A, B)
print("A ⊗ B =\n", K)
print("(A ⊗ B)^T =\n", K.T)
check("A ⊗ B agrees with the hand computation", K, [[2, 0, 6, 0], [-2, 3, -6, 9], [0, 0, -2, 0], [0, 0, 2, -3]])

print_step(2, "Transpose property (A ⊗ B)^T = A^T ⊗ B^T")
KT = np.kron(A.T, B.T)
print("A^T ⊗ B^T =\n", KT)
check("(A ⊗ B)^T = A^T ⊗ B^T", K.T, KT)

print_step(3, "Test the reversed version (A ⊗ B)^T =? B^T ⊗ A^T")
KR = np.kron(B.T, A.T)
print("B^T ⊗ A^T =\n", KR)
check("B^T ⊗ A^T agrees with the hand computation", KR, [[2, 0, -2, 0], [6, -2, -6, 2], [0, 0, 3, 0], [0, 0, 9, -3]])
check("(A ⊗ B)^T ≠ B^T ⊗ A^T", np.array_equal(K.T, KR), False)
check("(0,1) entries: -2 and 0", [K.T[0, 1], KR[0, 1]], [-2, 0])
check("B^T ⊗ A^T = (B ⊗ A)^T", KR, np.kron(B, A).T)

print_step(4, "Trace property tr(A ⊗ B) = tr(A)·tr(B)")
print("diagonal entries of A ⊗ B:", np.diag(K))
check("diagonal entries = (2, 3, -2, -3)", np.diag(K), [2, 3, -2, -3])
compare_print("tr(A ⊗ B)", np.trace(K), f"tr(A)·tr(B) = {np.trace(A)}·{np.trace(B)} = {np.trace(A) * np.trace(B)}")
check("tr(A ⊗ B) = tr(A)·tr(B) = 0", [np.trace(K), np.trace(A) * np.trace(B)], [0, 0])
check("tr(A) = 0, tr(B) = 5", [np.trace(A), np.trace(B)], [0, 5])

# ============================================================
print_header("Exercise 20 | Verifying the Mixed-Product Property")
# ============================================================
A = np.array([[1, 2], [2, 5]])
B = np.array([[1, 0], [-1, 1]])
C = np.array([[5, -2], [-2, 1]])
D = np.array([[2, 0], [2, 2]])

print_step(1, "Compute A ⊗ B and C ⊗ D")
AB, CD = np.kron(A, B), np.kron(C, D)
print("A ⊗ B =\n", AB)
print("C ⊗ D =\n", CD)
check("A ⊗ B agrees with the hand computation", AB, [[1, 0, 2, 0], [-1, 1, -2, 2], [2, 0, 5, 0], [-2, 2, -5, 5]])
check("C ⊗ D agrees with the hand computation", CD, [[10, 0, -4, 0], [10, 10, -4, -4], [-4, 0, 2, 0], [-4, -4, 2, 2]])

print_step(2, "Left-hand side (A ⊗ B)(C ⊗ D)")
lhs = AB @ CD
print("(A ⊗ B)(C ⊗ D) =\n", lhs)
check("left-hand side = 2 I_4", lhs, 2 * np.eye(4))
# Block multiplication: block (0,0) = B(5D) + (2B)(-2D) = BD
check("block (0,0) = BD", lhs[:2, :2], B @ (5 * D) + (2 * B) @ (-2 * D))

print_step(3, "Right-hand side (AC) ⊗ (BD)")
print("AC =\n", A @ C)
print("BD =\n", B @ D)
rhs = np.kron(A @ C, B @ D)
print("(AC) ⊗ (BD) =\n", rhs)
check("AC = I_2", A @ C, np.eye(2))
check("BD = 2 I_2", B @ D, 2 * np.eye(2))
check("(A ⊗ B)(C ⊗ D) = (AC) ⊗ (BD)", lhs, rhs)

print_step(4, "Observation: C = A^{-1}, D = 2 B^{-1}")
check("C = A^{-1}", C, np.linalg.inv(A))
check("D = 2 B^{-1}", D, 2 * np.linalg.inv(B))
check("(A ⊗ B)^{-1} = A^{-1} ⊗ B^{-1} = (C ⊗ D)/2",
      [np.linalg.inv(AB), np.kron(np.linalg.inv(A), np.linalg.inv(B))], [CD / 2, CD / 2])

# ============================================================
print_header("Exercise 21 | Tensor Products of Triangular Matrices and Identity Matrices")
# ============================================================
rng = np.random.default_rng(21)
sizes = [(1, 4), (2, 3), (3, 2), (4, 4), (5, 2)]   # (n, m): A is n×n, B is m×m

def is_upper(K):
    return np.allclose(K, np.triu(K))

print_step(1, "Entry formula (A ⊗ B)[i*m + k, j*m + l] = a_ij * b_kl")
A = rng.standard_normal((3, 3)); B = rng.standard_normal((2, 2)); n, m = 3, 2
K = np.kron(A, B)
ok = all(np.isclose(K[i * m + k, j * m + l], A[i, j] * B[k, l])
         for i in range(n) for j in range(n) for k in range(m) for l in range(m))
check("entry formula (n=3, m=2, random matrices)", ok, True)

print_step(2, "Part 1: upper triangular ⊗ upper triangular is upper triangular")
for n, m in sizes:
    A = np.triu(rng.standard_normal((n, n)))
    B = np.triu(rng.standard_normal((m, m)))
    check(f"n={n}, m={m}", is_upper(np.kron(A, B)), True)

print_step(3, "Part 2: unit upper triangular ⊗ unit upper triangular is unit upper triangular")
for n, m in sizes:
    A = np.triu(rng.standard_normal((n, n)), 1) + np.eye(n)
    B = np.triu(rng.standard_normal((m, m)), 1) + np.eye(m)
    K = np.kron(A, B)
    check(f"n={n}, m={m}", [is_upper(K), np.allclose(np.diag(K), 1)], [True, True])

print_step(4, "Part 3: I_n ⊗ I_m = I_{nm}")
for n, m in sizes:
    check(f"I_{n} ⊗ I_{m} = I_{n * m}", np.kron(np.eye(n), np.eye(m)), np.eye(n * m))

# ============================================================
print_header("Exercise 22 | Constructing the Natural Basis of R⁶ from Tensor Products")
# ============================================================
E2, E3, E6 = np.eye(2, dtype=int), np.eye(3, dtype=int), np.eye(6, dtype=int)
e = [E2[:, i] for i in range(2)]   # e_0, e_1
f = [E3[:, j] for j in range(3)]   # f_0, f_1, f_2
g = [E6[:, r] for r in range(6)]   # g_0, ..., g_5

print_step(1, "Compute each e_i ⊗ f_j and its index r = 3i + j")
indices = []
for i in range(2):
    for j in range(3):
        v = np.kron(e[i], f[j])
        r = 3 * i + j
        indices.append(r)
        print(f"  e_{i} ⊗ f_{j} = {v}  → g_{r}")
        check(f"e_{i} ⊗ f_{j} = g_{r}", v, g[r])

print_step(2, "The six tensor products run through exactly g_0, ..., g_5")
check("in lexicographic order the indices are 0,1,...,5 (bijection)", indices, list(range(6)))
check("divmod(r, 3) recovers (i, j)", [divmod(r, 3) for r in range(6)],
      [(i, j) for i in range(2) for j in range(3)])

print_step(3, "Rearranged in C order into a 2×3 matrix, this is E_{i,j} = e_i f_j^T")
for i in range(2):
    for j in range(3):
        check(f"reshape(e_{i} ⊗ f_{j}) = E_{i},{j}", np.kron(e[i], f[j]).reshape(2, 3), np.outer(e[i], f[j]))

# ============================================================
print_header("Exercise 23 | Tensor Products of Natural Bases")
# ============================================================
print_step(1, "Verify e_i ⊗ f_j = g_{i*n + j} for several (m, n)")
for m, n in [(1, 3), (2, 2), (2, 3), (3, 2), (3, 4), (4, 5)]:
    Em, En, Emn = np.eye(m), np.eye(n), np.eye(m * n)
    ok = all(np.array_equal(np.kron(Em[:, i], En[:, j]), Emn[:, i * n + j])
             for i in range(m) for j in range(n))
    idx = sorted(i * n + j for i in range(m) for j in range(n))
    check(f"m={m}, n={n}: e_i ⊗ f_j = g_(i*n+j)", ok, True)
    check(f"m={m}, n={n}: indices i*n+j run through exactly 0..{m * n - 1}", idx, list(range(m * n)))

print_step(2, "Division with remainder recovers (i, j): divmod(r, n) = (i, j)")
m, n = 3, 4
check("m=3, n=4: divmod(i*n + j, n) = (i, j)",
      [divmod(i * n + j, n) for i in range(m) for j in range(n)],
      [(i, j) for i in range(m) for j in range(n)])

# ============================================================
print_header("Exercise 24 | Tensor Products of the Natural Bases of Matrix Spaces")
# ============================================================
def unit_matrix(rows, cols, i, j):
    # the matrix with 1 in row i, column j and 0 elsewhere
    X = np.zeros((rows, cols), dtype=int)
    X[i, j] = 1
    return X

print_step(1, "Method 1: E_{i,j} ⊗ F_{p,q} = G_{ik+p, jl+q}")
for m, n, k, l in [(2, 2, 2, 2), (2, 3, 3, 2), (1, 2, 3, 1), (3, 2, 2, 4)]:
    ok, positions = True, []
    for i in range(m):
        for j in range(n):
            for p in range(k):
                for q in range(l):
                    K = np.kron(unit_matrix(m, n, i, j), unit_matrix(k, l, p, q))
                    r, s = i * k + p, j * l + q
                    ok &= np.array_equal(K, unit_matrix(m * k, n * l, r, s))
                    positions.append((r, s))
    check(f"(m,n,k,l)=({m},{n},{k},{l}): E ⊗ F = G_(ik+p, jl+q)", ok, True)
    check(f"(m,n,k,l)=({m},{n},{k},{l}): positions pairwise distinct and covering all {m * k}×{n * l} of them",
          sorted(positions), [(r, s) for r in range(m * k) for s in range(n * l)])

print_step(2, "Method 2: (e_i f_j^T) ⊗ (x_p y_q^T) = (e_i ⊗ x_p)(f_j ⊗ y_q)^T")
m, n, k, l = 2, 3, 3, 2
Im, In, Ik, Il = np.eye(m), np.eye(n), np.eye(k), np.eye(l)
ok = True
for i in range(m):
    for j in range(n):
        for p in range(k):
            for q in range(l):
                lhs = np.kron(np.outer(Im[:, i], In[:, j]), np.outer(Ik[:, p], Il[:, q]))
                rhs = np.outer(np.kron(Im[:, i], Ik[:, p]), np.kron(In[:, j], Il[:, q]))
                ok &= np.array_equal(lhs, rhs)
check("(m,n,k,l)=(2,3,3,2): the mixed-product form holds", ok, True)

# ============================================================
print_header("Exercise 25 | The Outer Product Is a Kronecker Product")
# ============================================================
rng = np.random.default_rng(25)

# --- Step 1 ---
print_step(1, "Concrete example: u = (1, 2)^T, v = (3, 4, 5)^T")
u = np.array([1, 2])
v = np.array([3, 4, 5])
U = u.reshape(-1, 1)      # m×1 matrix (a column vector)
Vt = v.reshape(1, -1)     # 1×n matrix (a row vector)
outer = U @ Vt            # outer product u v^T = |u⟩⟨v|
kron = np.kron(U, Vt)     # Kronecker product u ⊗ v^T
print("u v^T =\n", outer)
print("u ⊗ v^T =\n", kron)
check("u v^T = u ⊗ v^T", outer, kron)
check("u v^T = [[3,4,5],[6,8,10]]", outer, [[3, 4, 5], [6, 8, 10]])
check("v^T ⊗ u also equals u v^T", np.kron(Vt, U), outer)

# --- Step 2 ---
print_step(2, "Random vectors of various dimensions")
for m, n in [(1, 1), (2, 3), (4, 2), (5, 5)]:
    U = rng.standard_normal((m, 1))
    Vt = rng.standard_normal((1, n))
    check(f"m={m}, n={n}: u v^T = u ⊗ v^T", U @ Vt, np.kron(U, Vt))

# --- Step 3 ---
print_step(3, "Shape analysis: AB (requires n = p) and A⊗B have the same shape ⟺ n = p = 1")
ok_all = True
for m in range(1, 5):
    for n in range(1, 5):
        for q in range(1, 5):
            p = n                                   # AB is defined only when n = p
            same_shape = (m, q) == (m * p, n * q)   # AB is m×q, A⊗B is mp×nq
            ok_all &= (same_shape == (n == 1))
check("m, n, q ∈ {1,…,4}: same shape ⟺ n = p = 1", ok_all, True)
A2 = rng.standard_normal((2, 2))
B2 = rng.standard_normal((2, 2))
compare_print("two 2×2 matrices", f"shape of AB {(A2 @ B2).shape}, shape of A⊗B {np.kron(A2, B2).shape}",
              "(2, 2) and (4, 4), cannot be equal")

# ============================================================
print_header("Exercise 26 | Block Form of the Kronecker Product and the Commutation Matrix")
# ============================================================
rng = np.random.default_rng(26)

# --- Step 1 ---
print_step(1, "Symbolic form: represent the entries of A⊗B as string products")
a_sym = np.array([["a", "b"], ["c", "d"]])
b_sym = np.array([["x", "y"], ["z", "w"]])

def kron_sym(S, T):
    # (S⊗T)[2i+k, 2j+l] = S[i,j]·T[k,l], with string concatenation standing for the product
    return np.array([[S[r // 2, s // 2] + T[r % 2, s % 2] for s in range(4)]
                     for r in range(4)])

def normalize(X):
    # scalar multiplication is commutative: sort the letters of each product before comparing (ax and xa count as equal)
    return np.vectorize(lambda t: "".join(sorted(t)))(X)

AkB = kron_sym(a_sym, b_sym)
print("A⊗B =\n", AkB)

# --- Step 2 ---
print_step(2, "Permutation matrix P: row r has a 1 in column perm[r]")
perm = [0, 2, 1, 3]
P = np.eye(4, dtype=int)[perm]
print("P =\n", P)
check("P^T = P", P.T, P)
check("P P^T = I_4", P @ P.T, np.eye(4))

# --- Steps 3, 4 ---
print_step("3–4", "First swap rows 1 and 2, then swap columns 1 and 2")
PX = AkB[perm, :]           # left-multiply by P: swap rows 1 and 2
PXPt = PX[:, perm]          # right-multiply by P^T: swap columns 1 and 2
print("P(A⊗B) =\n", PX)
print("P(A⊗B)P^T =\n", PXPt)
expected = np.array([["ax", "bx", "ay", "by"],
                     ["cx", "dx", "cy", "dy"],
                     ["az", "bz", "aw", "bw"],
                     ["cz", "dz", "cw", "dw"]])
check("P(A⊗B)P^T matches the matrix in Step 4", (PXPt == expected).all(), True)

# --- Step 5 ---
print_step(5, "Compare with B⊗A (symbolic)")
BkA = kron_sym(b_sym, a_sym)
print("B⊗A =\n", BkA)
check("P(A⊗B)P^T = B⊗A (symbolic, products commute)", (normalize(PXPt) == normalize(BkA)).all(), True)

# --- Step 6 ---
print_step(6, "Numerical check: an integer example and random matrices; the index formula")
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
check("integer example: P(A⊗B)P^T = B⊗A", P @ np.kron(A, B) @ P.T, np.kron(B, A))
for t in range(5):
    A = rng.standard_normal((2, 2))
    B = rng.standard_normal((2, 2))
    check(f"random example {t}: P(A⊗B)P^T = B⊗A", P @ np.kron(A, B) @ P.T, np.kron(B, A))
ok = all(perm[2 * k + i] == 2 * i + k for i in range(2) for k in range(2))
check("P sends row 2k+i to 2i+k (digit swap)", ok, True)

# ============================================================
print_header("Exercise 27 | Matrices Satisfying A⊗X = X⊗A")
# ============================================================
import itertools
rng = np.random.default_rng(27)

A_cases = {
    "(1) A = I_2": np.eye(2),
    "(2) A = N":   np.array([[0., 1.], [0., 0.]]),
    "(3) A = [[0,1],[1,1]]": np.array([[0., 1.], [1., 1.]]),
}
# order of the unknowns vec(X) = (x00, x01, x10, x11); E[t] is the corresponding standard basis matrix
E = [np.eye(4)[t].reshape(2, 2) for t in range(4)]

for step, (name, A) in enumerate(A_cases.items(), start=1):
    print_step(step, name)
    # the 16×4 coefficient matrix of the linear mapping X ↦ A⊗X − X⊗A
    L = np.column_stack([(np.kron(A, Et) - np.kron(Et, A)).ravel() for Et in E])
    nullity = 4 - np.linalg.matrix_rank(L)
    null_vec = np.linalg.svd(L)[2][-1]          # a unit basis vector of the null space
    print("dimension of the null space =", nullity, "; basis vector (reshaped to 2×2) =\n", null_vec.reshape(2, 2))
    check("dimension of the null space = 1", nullity, 1)
    check("the null space is spanned by A (rank of [null, vec A] is 1)",
          np.linalg.matrix_rank(np.column_stack([null_vec, A.ravel()])), 1)
    c = rng.standard_normal()
    check("X = cA is a solution", np.kron(A, c * A), np.kron(c * A, A))
    # exhaustive search of the integer grid x_ij ∈ {−2,…,2}: the solutions are exactly cA, c ∈ {−2,…,2}
    sols = [np.array(x, float).reshape(2, 2) for x in itertools.product(range(-2, 3), repeat=4)
            if np.allclose(np.kron(A, np.array(x).reshape(2, 2)),
                           np.kron(np.array(x).reshape(2, 2), A))]
    print("solutions on the integer grid:", [S.ravel().tolist() for S in sols])
    check("number of integer-grid solutions = 5", len(sols), 5)
    check("every integer solution is a multiple of A",
          all(np.linalg.matrix_rank(np.column_stack([S.ravel(), A.ravel()])) <= 1 for S in sols), True)

print_step(4, "In general: for a random nonzero A (including 3×3), only X = cA")
for n in [2, 3]:
    A = rng.standard_normal((n, n))
    Es = [np.eye(n * n)[t].reshape(n, n) for t in range(n * n)]
    L = np.column_stack([(np.kron(A, Et) - np.kron(Et, A)).ravel() for Et in Es])
    check(f"n={n}: dimension of the null space = 1", n * n - np.linalg.matrix_rank(L), 1)

# ============================================================
print_header("Exercise 28 | Partial Traces of a Given 4×4 Block Matrix")
# ============================================================

# Partial trace tools (reused in Exercises 29 and 30)
def to_4d(M, dA, dB):
    # 4D tensor: 𝓜[i,j,k,l] = M[i·dB + k, j·dB + l] (block indices i,j; within-block indices k,l)
    return M.reshape(dA, dB, dA, dB).transpose(0, 2, 1, 3)

def ptrace_A(M, dA, dB):
    # [Tr_A(M)]_{kl} = Σ_i 𝓜[i,i,k,l]
    return np.einsum('iikl->kl', to_4d(M, dA, dB))

def ptrace_B(M, dA, dB):
    # [Tr_B(M)]_{ij} = Σ_k 𝓜[i,j,k,k]
    return np.einsum('ijkk->ij', to_4d(M, dA, dB))

M = np.array([[ 0,  1,  2,  3],
              [ 5,  7, 11, 13],
              [17, 19, 23, 29],
              [31, 37, 41, 43]])

# --- Step 1 ---
print_step(1, "The 4D tensor representation and the four blocks")
T = to_4d(M, 2, 2)
ok = all(T[i, j, k, l] == M[2 * i + k, 2 * j + l]
         for i in range(2) for j in range(2) for k in range(2) for l in range(2))
check("𝓜[i,j,k,l] = M[2i+k, 2j+l]", ok, True)
M00, M01, M10, M11 = M[:2, :2], M[:2, 2:], M[2:, :2], M[2:, 2:]
check("𝓜[0,0] = M_00 = [[0,1],[5,7]]", T[0, 0], [[0, 1], [5, 7]])
check("𝓜[0,1] = M_01 = [[2,3],[11,13]]", T[0, 1], [[2, 3], [11, 13]])
check("𝓜[1,0] = M_10 = [[17,19],[31,37]]", T[1, 0], [[17, 19], [31, 37]])
check("𝓜[1,1] = M_11 = [[23,29],[41,43]]", T[1, 1], [[23, 29], [41, 43]])

# --- Step 2 ---
print_step(2, "Tr_A(M): set i = j and sum")
TrA = ptrace_A(M, 2, 2)
print("Tr_A(M) =\n", TrA)
check("Tr_A(M) = [[23,30],[46,50]]", TrA, [[23, 30], [46, 50]])
check("Tr_A(M) = M_00 + M_11", TrA, M00 + M11)

# --- Step 3 ---
print_step(3, "Tr_B(M): set k = l and sum")
TrB = ptrace_B(M, 2, 2)
print("Tr_B(M) =\n", TrB)
check("Tr_B(M) = [[7,15],[54,66]]", TrB, [[7, 15], [54, 66]])
check("Tr_B(M) = [[tr M_00, tr M_01],[tr M_10, tr M_11]]", TrB,
      [[np.trace(M00), np.trace(M01)], [np.trace(M10), np.trace(M11)]])

# --- Step 4 ---
print_step(4, "The three traces")
compare_print("tr(Tr_A M), tr(Tr_B M), tr(M)",
              (int(np.trace(TrA)), int(np.trace(TrB)), int(np.trace(M))), "(73, 73, 73)")
check("tr(Tr_A M) = 73", np.trace(TrA), 73)
check("tr(Tr_B M) = 73", np.trace(TrB), 73)
check("tr(M) = 73", np.trace(M), 73)
check("Σ_{i,k} 𝓜[i,i,k,k] = 73", np.einsum('iikk->', T), 73)

# ============================================================
print_header("Exercise 29 | Block Formulas for the Partial Traces of a 2×2 Block Matrix")
# ============================================================
rng = np.random.default_rng(29)

for step, n in enumerate([1, 2, 3, 5], start=1):
    print_step(step, f"Random 2×2 block matrix, each block {n}×{n}")
    M = rng.standard_normal((2 * n, 2 * n))
    # extract the four blocks directly by slicing (independently of the einsum-based 4D tensor definition)
    M00, M01 = M[:n, :n], M[:n, n:]
    M10, M11 = M[n:, :n], M[n:, n:]
    check("Tr_A(M) = M_00 + M_11", ptrace_A(M, 2, n), M00 + M11)
    check("Tr_B(M) = [[tr M_ij]]", ptrace_B(M, 2, n),
          [[np.trace(M00), np.trace(M01)], [np.trace(M10), np.trace(M11)]])
    check("Tr_A(M) has shape n×n, Tr_B(M) has shape 2×2",
          (ptrace_A(M, 2, n).shape, ptrace_B(M, 2, n).shape) == ((n, n), (2, 2)), True)

# ============================================================
print_header("Exercise 30 | Partial Traces of a Kronecker Product and the Multiplicativity of the Trace")
# ============================================================
rng = np.random.default_rng(30)

# A is n×n (system A, block indices), B is m×m (system B, within-block indices)
for step, (n, m) in enumerate([(2, 3), (3, 2), (1, 4), (4, 1), (3, 5)], start=1):
    print_step(step, f"A is {n}×{n}, B is {m}×{m}")
    A = rng.standard_normal((n, n))
    B = rng.standard_normal((m, m))
    K = np.kron(A, B)
    check("𝓜[i,j,k,l] = a_ij · b_kl", to_4d(K, n, m), np.einsum('ij,kl->ijkl', A, B))
    check("Tr_A(A⊗B) = tr(A) B", ptrace_A(K, n, m), np.trace(A) * B)
    check("Tr_B(A⊗B) = tr(B) A", ptrace_B(K, n, m), np.trace(B) * A)
    check("tr(A⊗B) = tr(A) tr(B)", np.trace(K), np.trace(A) * np.trace(B))
    # a general M (not of product form): the partial trace preserves the trace
    M = rng.standard_normal((n * m, n * m))
    check("general M: tr(Tr_A M) = tr(M)", np.trace(ptrace_A(M, n, m)), np.trace(M))
    check("general M: tr(Tr_B M) = tr(M)", np.trace(ptrace_B(M, n, m)), np.trace(M))

print_step(6, "Product state: once a factor of trace 1 is traced out, the other factor is recovered in full")
rhoA = np.diag([0.25, 0.75])
rhoB = np.array([[0.5, 0.5], [0.5, 0.5]])
rho = np.kron(rhoA, rhoB)
check("Tr_A(ρ_A⊗ρ_B) = ρ_B", ptrace_A(rho, 2, 2), rhoB)
check("Tr_B(ρ_A⊗ρ_B) = ρ_A", ptrace_B(rho, 2, 2), rhoA)