第 9 章 進階主題與現代發展#

"""
minialg.py —— 用 SymPy 建構自己的小代數系統
配套：09.5 專題實作
"""
import itertools
import random as _random
import sympy as sp

# ============================================================
# 第 0 部分：公理檢查器（把數學定義寫成程式）
# ============================================================

def check_axioms(elements, op, identity, commutative=False, samples=None):
    """檢查 (elements, op) 是否構成群。
    samples=None 時做全枚舉（適合小集合），
    否則做 samples 次隨機抽樣測試（適合大集合）。"""
    E = list(elements)
    if samples is None or len(E) <= 40:
        test_pairs = list(itertools.product(E, repeat=2))
    else:
        test_pairs = [tuple(_random.sample(E, 2)) for _ in range(samples)]
    # 封閉性 + 結合律
    for a, b in test_pairs:
        if op(a, b) not in E:
            return f"封閉性失敗：{a} * {b} = {op(a,b)} 不在集合中"
    triples = (list(itertools.product(E, repeat=3))
               if samples is None or len(E) <= 25 else
               [tuple(_random.sample(E, 3)) for _ in range(samples)])
    for a, b, c in triples:
        if op(op(a, b), c) != op(a, op(b, c)):
            return f"結合律失敗於 ({a},{b},{c})"
    # 單位元
    if any(op(identity, a) != a or op(a, identity) != a for a in E):
        return "單位元失敗"
    # 反元素
    for a in E:
        if not any(op(a, b) == identity for b in E):
            return f"反元素失敗：{a}"
    if commutative:
        for a, b in test_pairs:
            if op(a, b) != op(b, a):
                return f"交換律失敗於 ({a},{b})"
    return "OK：滿足群公理" + ("（含交換律）" if commutative else "")

# ============================================================
# 第 1 部分：循環群 C_n 與二面體群 D_n（運算表風格）
# ============================================================

class CyclicGroup:
    """C_n = Z/nZ，元素為 0..n-1，運算為 mod n 加法"""
    def __init__(self, n):
        self.n = n
        self.elements = list(range(n))
        self.identity = 0
    def op(self, a, b):
        return (a + b) % self.n
    def order_of(self, g):
        x, k = g, 1
        while x != 0:
            x = self.op(x, g); k += 1
        return k

class DihedralGroup:
    """D_n = <r, s | r^n = s^2 = e, srs = r^{-1}>，
    元素以 (旋轉數, 鏡像) 表示"""
    def __init__(self, n):
        self.n = n
        self.elements = [(k, m) for m in (0, 1) for k in range(n)]
        self.identity = (0, 0)
    def op(self, a, b):
        k1, m1 = a; k2, m2 = b
        if m2 == 0:
            return ((k1 + k2) % self.n, m1)
        return ((k1 - k2) % self.n, 1 - m1)
    def order_of(self, g):
        x, k = self.identity, 0
        while True:
            x = self.op(x, g); k += 1
            if x == self.identity:
                return k

# --- 驗證 ---
C6 = CyclicGroup(6)
print("[C6]", check_axioms(C6.elements, C6.op, C6.identity,
                          commutative=True))
D4 = DihedralGroup(4)
print("[D4]", check_axioms(D4.elements, D4.op, D4.identity))
print("D4 元素階數：", sorted(D4.order_of(g) for g in D4.elements))
# Lagrange：所有階整除 |D4| = 8
assert all(D4.order_of(g) in (1, 2, 4) for g in D4.elements)

# ============================================================
# 第 2 部分：商環 Z[x]/(x^2 + 1) —— 高斯整數的化身
# ============================================================

class QuotientZxQuadratic:
    """元素 a + b·x，商規則 x^2 = -1；係數為整數"""
    def __init__(self, a, b):
        self.a, self.b = a, b
    def __add__(s, t):  return QuotientZxQuadratic(s.a + t.a, s.b + t.b)
    def __mul__(s, t):
        # (a + bx)(c + dx) = ac + (ad+bc)x + bd x^2 = (ac - bd) + (ad+bc)x
        return QuotientZxQuadratic(s.a*t.a - s.b*t.b,
                                   s.a*t.b + s.b*t.a)
    def __eq__(s, t):   return (s.a, s.b) == (t.a, t.b)
    def __repr__(s):    return f"({s.a} + {s.b}x)"

# 驗證：與複數對照（本體就是高斯整數 Z[i]）
u = QuotientZxQuadratic(2, 3)
v = QuotientZxQuadratic(1, -2)
print("\n[商環 Z[x]/(x^2+1)]")
print("  (2+3x)(1-2x) =", u * v)          # (8 - x) → 即 8 - i
z1, z2 = 2 + 3j, 1 - 2j
print("  複數對照 (2+3i)(1-2i) =", z1 * z2)   # (8-1j)

# ============================================================
# 第 3 部分：Galois 域 GF(2^8) —— AES 的數學基礎
# ============================================================

M = 0x11B   # m(x) = x^8 + x^4 + x^3 + x + 1 的整數編碼 100011011

def gf_add(a, b):
    """特徵 2 的加法：XOR"""
    return a ^ b

def gf_mul(a, b):
    """多項式乘法 + 對 m(x) 取餘（俄國農夫乘法）"""
    p = 0
    while b:
        if b & 1:
            p ^= a
        b >>= 1
        a <<= 1
        if a & 0x100:           # x^8 位出現 → 減去 m(x)
            a ^= M
    return p

def gf_pow(a, k):
    p = 1
    while k:
        if k & 1:
            p = gf_mul(p, a)
        a = gf_mul(a, a)
        k >>= 1
    return p

def gf_inv(a):
    """反元素：乘法群是 255 階循環群 → a^254 = a^{-1}
    （Fermat 小定理在有限體的版本，見 4.7 節）"""
    assert 0 < a < 256
    return gf_pow(a, 254)

def gf_is_field_check():
    """驗證 GF(2^8) 公理：加法群 + 非零元乘法群（抽樣 + 定理）"""
    # 加法交換群：XOR 顯然，抽驗結合律
    import random
    for _ in range(1000):
        a, b, c = (random.randrange(256) for _ in range(3))
        assert gf_add(gf_add(a, b), c) == gf_add(a, gf_add(b, c))
        assert gf_mul(gf_mul(a, b), c) == gf_mul(a, gf_mul(b, c))
        assert gf_mul(a, b) == gf_mul(b, a)
    # 非零元均可逆：由「乘法群循環」定理，驗證一個 255 階生成元即可
    g = 3
    assert gf_pow(g, 255) == 1 and any(gf_pow(g, k) == 1
                                       for k in range(1, 255)) is False
    order = 1
    x = g
    while x != 1:
        x = gf_mul(x, g); order += 1
    print("  生成元 3 的階 =", order, "（= 255，乘法群循環）")
    return "OK：GF(2^8) 體公理通過"

print("\n[GF(2^8)] m(x) = x^8+x^4+x^3+x+1 (0x11B)")
print(" 0x57 + 0x83 =", hex(gf_add(0x57, 0x83)))   # 0xd4（AES 論文原例）
print(" 0x57 * 0x83 =", hex(gf_mul(0x57, 0x83)))   # 0xc1
print(" (0x57)^-1 =", hex(gf_inv(0x57)))           # 0x8d
print(" 驗證 0x57·0x8d =", hex(gf_mul(0x57, 0x8d)))  # 0x1
print(gf_is_field_check())

# --- AES S-box：仿射變換層 ---
SBOX = []
for i in range(256):
    inv = gf_inv(i) if i else 0
    b = inv
    # 仿射：b_j' = b_j ^ b_{j+1} ^ b_{j+2} ^ b_{j+3} ^ b_{j+4} ^ c_j
    for _ in range(4):
        inv = ((inv << 1) | (inv >> 7)) & 0xFF    # 循環左移 xtime
    SBOX.append(b ^ inv ^ 0x63)
print(" SBOX[0x53] =", hex(SBOX[0x53]))            # 0xed（AES 標準）
print(" SBOX[0x00] =", hex(SBOX[0x00]))            # 0x63

# ============================================================
# 第 4 部分：用 SymPy 當參考答案，交叉驗證
# ============================================================

x = sp.Symbol('x')
# SymPy 驗證 m(x) 不可約（GF(2) 上）：
m_poly = sp.Poly(x**8 + x**4 + x**3 + x + 1, x, domain=sp.GF(2))
print("\n[SymPy 交叉驗證]")
print(" m(x) 在 GF(2) 上不可約？", m_poly.is_irreducible)   # True

# 用 sp.GF 直接做 GF(2^8) 乘法，比對 gf_mul
F = sp.GF(2**8, modulus=x**8 + x**4 + x**3 + x + 1)
a_ = F(0x57); b_ = F(0x83)
print(" SymPy GF(2^8): 0x57·0x83 =", hex(int(a_ * b_)))     # 0xc1
print(" SymPy GF(2^8): (0x57)^-1 =", hex(int(a_**-1)))      # 0x8d

# 全表比對：SymPy vs. 自製實作
ok = all(int(F(i) * F(j)) == gf_mul(i, j)
         for i in range(256) for j in range(256))
print(" 256×256 乘法表全比對：", "一致 ✓" if ok else "不一致 ✗")

# 乘法群循環性的 SymPy 驗證：找 3 的階
print(" SymPy n_order(3, 255) =",
      sp.n_order(3, 255))   # 255 —— 呼應 9.4 節

# 商環 Z[x]/(x^2+1) 也可用 SymPy 商環對照
R = sp.QQ[x]
I = R.old_poly_ring(x).quotient_ring(x**2 + 1)
print(" SymPy 商環範例：", I.one, "，此處以自製類別為主")

# ============================================================
# 第 5 部分（挑戰）：自製橢圓曲線點類別，與 9.2 節互通
# ============================================================

class ECPoint:
    """橢圓曲線 y^2 = x^3 + ax + b over GF(p) 的點；None 代表 O"""
    def __init__(self, pt, a, p):
        self.pt = pt; self.a = a; self.p = p
    def __add__(s, t):
        a, p = s.a, s.p
        if s.pt is None: return t
        if t.pt is None: return s
        (x1, y1), (x2, y2) = s.pt, t.pt
        if x1 == x2 and (y1 + y2) % p == 0:
            return ECPoint(None, a, p)
        m = ((3*x1**2 + a) * pow(2*y1, -1, p)) % p if s.pt == t.pt \
            else ((y2 - y1) * pow(x2 - x1, -1, p)) % p
        x3 = (m*m - x1 - x2) % p
        return ECPoint((x3, (m*(x1 - x3) - y1) % p), a, p)
    def __rmul__(s, k):          # 純量乘法
        R, B = ECPoint(None, s.a, s.p), s
        while k:
            if k & 1: R = R + B
            B = B + B; k >>= 1
        return R
    def __repr__(s):
        return "O" if s.pt is None else f"({s.pt[0]}, {s.pt[1]})"

a, b, p = 2, 3, 97
P = ECPoint((3, 6), a, p)
Q = 5 * P
print("\n[自製橢圓曲線] y^2 = x^3+2x+3 over GF(97)")
print(" 5P =", Q)                    # (85, 68)，與 9.2 節一致

[C6] OK：滿足群公理（含交換律）
[D4] OK：滿足群公理
D4 元素階數： [1, 2, 2, 2, 2, 4, 4, 4]
...
 0x57 * 0x83 = 0xc1
 (0x57)^-1 = 0x8d
 驗證 0x57·0x8d = 0x1
 生成元 3 的階 = 255（= 255，乘法群循環）
OK：GF(2^8) 體公理通過
 SBOX[0x53] = 0xed
 SBOX[0x00] = 0x63

[SymPy 交叉驗證]
 m(x) 在 GF(2) 上不可約？ True
 SymPy GF(2^8): 0x57·0x83 = 0xc1
 256×256 乘法表全比對： 一致 ✓
 SymPy n_order(3, 255) = 255

[自製橢圓曲線] y^2 = x^3+2x+3 over GF(97)
 5P = (85, 68)