Erasure Coding & Reed–Solomon — Practice Tasks¶

import os

def xor_blocks(blocks):
    """XOR a list of equal-length byte strings together."""
    out = bytearray(blocks[0])
    for b in blocks[1:]:
        for i in range(len(out)):
            out[i] ^= b[i]
    return bytes(out)

def raid5_encode(data_blocks):
    """data_blocks: list of k equal-length bytes. Returns k+1 blocks (data..., parity)."""
    parity = xor_blocks(data_blocks)
    return list(data_blocks) + [parity]

def raid5_reconstruct(stored, missing_index):
    """stored: list of k+1 blocks with stored[missing_index] = None. Returns the missing block."""
    present = [b for j, b in enumerate(stored) if j != missing_index]
    # The XOR of all present blocks equals the missing block, because the full
    # XOR of all k+1 blocks (data XOR ... XOR parity) is zero by construction.
    return xor_blocks(present)

def test_task1():
    for k in (2, 3, 4, 8):
        L = 64
        data = [os.urandom(L) for _ in range(k)]
        stored = raid5_encode(data)
        for i in range(k + 1):
            damaged = list(stored)
            lost = damaged[i]
            damaged[i] = None
            rec = raid5_reconstruct(damaged, i)
            assert rec == lost, f"k={k} index={i} mismatch"
            # Full byte-exact stripe check
            damaged[i] = rec
            assert damaged == stored
    print("Task 1 OK")

test_task1()

def replication(k, f):
    r = f + 1
    n = r * k
    return dict(scheme=f"{r}x-replication", n=n, overhead=n / k, usable=k / n)

def erasure(k, f):
    n = k + f
    return dict(scheme=f"RS({k+f},{k})", n=n, overhead=n / k, usable=k / n)

def gb_stored(overhead, user_gb=1.0):
    return overhead * user_gb

def test_task2():
    k = 10
    print(f"{'f':>2} {'scheme':>18} {'n':>3} {'overhead':>9} {'usable':>7} {'GB/GB':>7}")
    for f in (1, 2, 3):
        rep = replication(k, f)
        ec = erasure(k, f)
        for s in (rep, ec):
            print(f"{f:>2} {s['scheme']:>18} {s['n']:>3} {s['overhead']:>9.3f} "
                  f"{s['usable']:>7.3f} {gb_stored(s['overhead']):>7.3f}")
        assert ec['overhead'] < rep['overhead']
        assert (rep['n'] - k * (f + 1) == 0)  # sanity on replication count
    # Hand check
    assert abs(replication(10, 2)['overhead'] - 3.0) < 1e-9
    assert abs(erasure(10, 2)['overhead'] - 1.2) < 1e-9
    print("Task 2 OK")

test_task2()

PRIM = 0x11D  # x^8 + x^4 + x^3 + x^2 + 1
GEN = 3       # primitive element

def _russian_peasant_mul(a, b):
    """Reference GF(2^8) multiply without tables (for cross-checking)."""
    result = 0
    while b:
        if b & 1:
            result ^= a
        b >>= 1
        a <<= 1
        if a & 0x100:
            a ^= PRIM
    return result & 0xFF

# Build log / antilog tables.
EXP = [0] * 512   # antilog: EXP[i] = g^i  (doubled so we can index without mod)
LOG = [0] * 256
def _build_tables():
    x = 1
    for i in range(255):
        EXP[i] = x
        LOG[x] = i
        x = _russian_peasant_mul(x, GEN)
    for i in range(255, 512):
        EXP[i] = EXP[i - 255]
    LOG[0] = -1  # undefined; guard against use
_build_tables()

def gadd(a, b):
    return a ^ b

def gmul(a, b):
    if a == 0 or b == 0:
        return 0
    return EXP[LOG[a] + LOG[b]]          # log addition, no explicit mod (table doubled)

def gdiv(a, b):
    if b == 0:
        raise ZeroDivisionError("division by zero in GF(2^8)")
    if a == 0:
        return 0
    return EXP[(LOG[a] - LOG[b]) % 255]

def ginv(a):
    if a == 0:
        raise ZeroDivisionError("0 has no inverse in GF(2^8)")
    return EXP[(255 - LOG[a]) % 255]

def test_task4():
    # 1. addition
    for a in range(256):
        assert gadd(a, a) == 0
    # 5. mul matches reference
    import random
    for _ in range(5000):
        a, b = random.randrange(256), random.randrange(256)
        assert gmul(a, b) == _russian_peasant_mul(a, b)
    # 2. associativity, 3. distributivity
    for _ in range(5000):
        a, b, c = (random.randrange(256) for _ in range(3))
        assert gmul(gmul(a, b), c) == gmul(a, gmul(b, c))
        assert gmul(a, gadd(b, c)) == gadd(gmul(a, b), gmul(a, c))
    # 4. inverses
    for a in range(1, 256):
        assert gmul(a, ginv(a)) == 1
    try:
        ginv(0); assert False
    except ZeroDivisionError:
        pass
    print("Task 4 OK")

test_task4()

from itertools import combinations
import random

P = 257  # small prime field; symbols live in 0..P-1

def inv_mod(a, p=P):
    return pow(a, p - 2, p)   # Fermat inverse, p prime

def rs_prime_encode(data, n, p=P):
    """data: k symbols (poly coefficients). Evaluate at points 1..n => n codeword symbols.
    Returns list of (x, value). Systematic-ish: we keep x's so decode knows the points."""
    k = len(data)
    xs = list(range(1, n + 1))
    code = []
    for x in xs:
        # Horner evaluation of sum_{j} data[j] * x^j  mod p
        acc = 0
        for coef in reversed(data):
            acc = (acc * x + coef) % p
        code.append((x, acc))
    return code

def lagrange_interpolate_coeffs(points, p=P):
    """Given k points, return the k polynomial coefficients (constant-term first)."""
    k = len(points)
    # Build coefficients by summing y_i * L_i(x), each L_i is a degree k-1 poly.
    result = [0] * k
    for i in range(k):
        xi, yi = points[i]
        # numerator poly = prod_{j != i} (x - x_j); denominator = prod (x_i - x_j)
        num = [1]  # polynomial "1"
        den = 1
        for j in range(k):
            if j == i:
                continue
            xj = points[j][0]
            # multiply num by (x - xj): poly mult by [-xj, 1]
            new = [0] * (len(num) + 1)
            for d, c in enumerate(num):
                new[d] = (new[d] - c * xj) % p
                new[d + 1] = (new[d + 1] + c) % p
            num = new
            den = (den * (xi - xj)) % p
        scale = (yi * inv_mod(den % p, p)) % p
        for d in range(len(num)):
            result[d] = (result[d] + scale * num[d]) % p
    return result

def rs_prime_decode(survivors, k, p=P):
    """survivors: list of (x, value) of length >= k. Returns original k coefficients."""
    pts = survivors[:k]
    return lagrange_interpolate_coeffs(pts, p)

def test_task5():
    for k in (2, 3, 4):
        n = k + 3
        data = [random.randrange(P - 1) for _ in range(k)]
        code = rs_prime_encode(data, n)
        idxs = list(range(n))
        for erased in combinations(idxs, n - k):
            survivors = [code[i] for i in idxs if i not in erased]
            rec = rs_prime_decode(survivors, k)
            assert rec == data, f"k={k} erased={erased}: {rec} != {data}"
    print("Task 5 OK")

test_task5()

import os, random
from itertools import combinations
# Assumes gadd, gmul, gdiv, ginv, _build_tables from Task 4 are in scope.

def cauchy_matrix(rows, cols):
    """(rows x cols) Cauchy matrix over GF(2^8): C[i][j] = 1 / (x_i XOR y_j),
    with x_i and y_j disjoint sets of field elements (guarantees invertibility of
    every square submatrix)."""
    xs = list(range(rows))             # 0 .. rows-1
    ys = list(range(rows, rows + cols))  # rows .. rows+cols-1  (disjoint from xs)
    assert rows + cols <= 256, "field too small for these parameters"
    M = [[gdiv(1, gadd(xs[i], ys[j])) for j in range(cols)] for i in range(rows)]
    return M

def generator_matrix(n, k):
    """n x k systematic generator: identity on top, Cauchy parity below."""
    G = [[1 if i == j else 0 for j in range(k)] for i in range(k)]  # identity
    G += cauchy_matrix(n - k, k)                                     # parity rows
    return G

def mat_vec_symbols(M, data):
    """M: r x k matrix of field elements. data: k symbols (each a bytes of length L).
    Returns r symbols. Each output symbol o = sum_j M[i][j] * data[j] (GF byte-wise)."""
    L = len(data[0])
    out = []
    for row in M:
        acc = bytearray(L)
        for j, coef in enumerate(row):
            if coef == 0:
                continue
            dj = data[j]
            for b in range(L):
                acc[b] ^= gmul(coef, dj[b])
        out.append(bytes(acc))
    return out

def crs_encode(data, n, k):
    G = generator_matrix(n, k)
    return mat_vec_symbols(G, data)   # n symbols; first k == data (systematic)

def gf_solve(A, rhs):
    """Solve A x = rhs over GF(2^8). A: k x k field matrix; rhs: list of k symbols
    (each bytes of length L). Returns list of k symbols. Gaussian elimination."""
    k = len(A)
    L = len(rhs[0])
    A = [row[:] for row in A]
    rhs = [bytearray(s) for s in rhs]
    for col in range(k):
        # find pivot
        piv = next((r for r in range(col, k) if A[r][col] != 0), None)
        if piv is None:
            raise ValueError("singular matrix (should not happen for Cauchy submatrix)")
        if piv != col:
            A[col], A[piv] = A[piv], A[col]
            rhs[col], rhs[piv] = rhs[piv], rhs[col]
        # normalize pivot row
        inv = ginv(A[col][col])
        for c in range(col, k):
            A[col][c] = gmul(A[col][c], inv)
        for b in range(L):
            rhs[col][b] = gmul(rhs[col][b], inv)
        # eliminate this column from all other rows
        for r in range(k):
            if r == col or A[r][col] == 0:
                continue
            factor = A[r][col]
            for c in range(col, k):
                A[r][c] ^= gmul(factor, A[col][c])
            for b in range(L):
                rhs[r][b] ^= gmul(factor, rhs[col][b])
    return [bytes(s) for s in rhs]

def crs_decode(survivor_indices, survivor_symbols, n, k):
    """Recover original k data symbols from any k survivors."""
    G = generator_matrix(n, k)
    A = [G[i] for i in survivor_indices[:k]]
    rhs = survivor_symbols[:k]
    return gf_solve(A, rhs)

def test_task7():
    for (n, k) in [(6, 4), (9, 6), (14, 10)]:
        L = 37
        data = [os.urandom(L) for _ in range(k)]
        code = crs_encode(data, n, k)
        assert code[:k] == data  # systematic property
        idxs = list(range(n))

        def check(erased):
            survivors = [i for i in idxs if i not in erased]
            sym = [code[i] for i in survivors]
            rec = crs_decode(survivors, sym, n, k)
            assert rec == data, f"(n={n},k={k}) erased={erased}"

        if n <= 9:
            for erased in combinations(idxs, n - k):
                check(set(erased))
        else:
            for _ in range(400):
                check(set(random.sample(idxs, n - k)))
    print("Task 7 OK")

test_task7()