06.5 用 SymPy 操作特徵標與對稱函數#

import sympy as sp

x1, x2, x3 = sp.symbols('x1 x2 x3')
xs = [x1, x2, x3]

# SymPy 內建對稱多項式: sympy.polys.polytools.symmetric_poly
e1 = sp.symmetric_poly(1, *xs)
e2 = sp.symmetric_poly(2, *xs)
e3 = sp.symmetric_poly(3, *xs)
print("e1 =", e1)
print("e2 =", e2)
print("e3 =", e3)

# Viète: (t-x1)(t-x2)(t-x3) 的係數
t = sp.Symbol('t')
poly = sp.Poly(sp.expand((t - x1)*(t - x2)*(t - x3)), t)
print("展開:", poly.as_expr())
print("t^2 係數 =", sp.expand(-e1), "->", poly.coeff_monomial(t**2))
print("t^1 係數 =", sp.expand(e2),  "->", poly.coeff_monomial(t))
print("常數項   =", sp.expand(-e3), "->", poly.coeff_monomial(1))

e1 = x1 + x2 + x3
e2 = x1*x2 + x1*x3 + x2*x3
e3 = x1*x2*x3
展開: t**3 + t**2*(-x1 - x2 - x3) + t*(x1*x2 + x1*x3 + x2*x3) - x1*x2*x3
t^2 係數 = -x1 - x2 - x3 -> -x1 - x2 - x3
t^1 係數 = x1*x2 + x1*x3 + x2*x3 -> x1*x2 + x1*x3 + x2*x3
常數項   = -x1*x2*x3 -> -x1*x2*x3

import sympy as sp

xs = sp.symbols('x1:5')
n = 4
e = [sp.Integer(1)] + [sp.symmetric_poly(k, *xs) for k in range(1, n + 1)]
p = [sp.Integer(0)] + [sum(x**k for x in xs) for k in range(1, n + 1)]

for k in range(1, n + 1):
    expr = sp.expand(p[k] - e[1]*p[k-1] + (e[2]*p[k-2] if k >= 2 else 0)
                     - (e[3]*p[k-3] if k >= 3 else 0) + (e[4]*p[k-4] if k >= 4 else 0)
                     - (-1)**k * k * e[k])   # 符號整理: 統一形式
    # 直接驗證標準形式:
    s = sp.expand(p[k] + sum((-1)**j * e[j] * p[k-j] for j in range(1, k)) + (-1)**k * k * e[k])
    print(f"k = {k}: Newton 恆等式值 = {sp.simplify(s)}")

k = 1: Newton 恆等式值 = 0
k = 2: Newton 恆等式值 = 0
k = 3: Newton 恆等式值 = 0
k = 4: Newton 恆等式值 = 0

import sympy as sp
from math import factorial

def complete_h(k, xs):
    # h_k = sum over |mu|=k of monomial sym poly; 直接枚舉弱增指標
    m = len(xs)
    from itertools import combinations_with_replacement
    total = sp.Integer(0)
    for idx in combinations_with_replacement(range(m), k):
        term = sp.Integer(1)
        for i in idx:
            term *= xs[i]
        total += term
    return total

def schur(lam, xs):
    k = len(lam)
    H = sp.Matrix([[complete_h(lam[i] - i + j, xs) if lam[i] - i + j >= 0 else 0
                    for j in range(k)] for i in range(k)])
    return sp.expand(H.det())

def power_sum(k, xs):
    return sum(x**k for x in xs)

xs = sp.symbols('x1:4')
s21 = schur([2, 1], xs)
print("s_(2,1) =", s21)

# Jacobi–Trudi 行列式展開 = h_2*h_1 - h_3；另一個古典寫法 = e_1*e_2 - e_3
h1, h2, h3 = complete_h(1, xs), complete_h(2, xs), complete_h(3, xs)
e1 = sp.symmetric_poly(1, *xs)
e2 = sp.symmetric_poly(2, *xs)
e3 = sp.symmetric_poly(3, *xs)
print("h2*h1 - h3 =", sp.expand(h2*h1 - h3))
print("e1*e2 - e3 =", sp.expand(e1*e2 - e3))
print("相等:", sp.expand(s21 - (h2*h1 - h3)) == 0,
      sp.expand(s21 - (e1*e2 - e3)) == 0)

s_(2,1) = x1**2*x2 + x1**2*x3 + x1*x2**2 + 2*x1*x2*x3 + x1*x3**2 + x2**2*x3 + x2*x3**2
h2*h1 - h3 = x1**2*x2 + x1**2*x3 + x1*x2**2 + 2*x1*x2*x3 + x1*x3**2 + x2**2*x3 + x2*x3**2
e1*e2 - e3 = x1**2*x2 + x1**2*x3 + x1*x2**2 + 2*x1*x2*x3 + x1*x3**2 + x2**2*x3 + x2*x3**2
相等: True True

import sympy as sp
from math import factorial
from collections import Counter
from itertools import combinations_with_replacement

# ---- 共用工具（與範例 3 相同，此處列出使本段可獨立執行）----
def complete_h(k, xs):
    total = sp.Integer(0)
    for idx in combinations_with_replacement(range(len(xs)), k):
        term = sp.Integer(1)
        for i in idx:
            term *= xs[i]
        total += term
    return total

def schur(lam, xs):
    """Jacobi–Trudi: s_lam = det(h_{lam_i - i + j})。"""
    k = len(lam)
    H = sp.Matrix([[complete_h(lam[i] - i + j, xs) if lam[i] - i + j >= 0 else 0
                    for j in range(k)] for i in range(k)])
    return sp.expand(H.det())

def power_sum_of(k, xs):
    """冪和 p_k = sum_i x_i^k。"""
    return sum(x**k for x in xs)

def z_mu(mu):
    """z_mu = prod_i i^{m_i} m_i!（分割 mu 的 z 值）。"""
    c = Counter(mu)
    return sp.prod(i**cnt * factorial(cnt) for i, cnt in c.items())

def expand_schur_in_power(lam, xs):
    """解線性方程組: s_lam = sum_mu a_mu p_mu, 得 a_mu; 特徵標 chi(lam; mu) = a_mu * z_mu"""
    n = sum(lam)
    # 枚舉 n 的所有分割(循環型)
    mus = sorted({tuple(sorted([k for k in p for _ in range(p[k])], reverse=True))
                  for p in sp.utilities.iterables.partitions(n)}, reverse=True)

    # p_mu 作為多項式基
    p_exprs = {}
    for mu in mus:
        e = sp.Integer(1)
        for k in mu:
            e *= power_sum_of(k, xs)
        p_exprs[mu] = sp.expand(e)
    s_expr = schur(list(lam), xs)

    # 線性求解: s_expr = sum_mu a_mu p_exprs[mu]（係數比較，多餘方程但相容且唯一解）
    from sympy.polys.polytools import Poly
    polys = [Poly(p_exprs[mu], *xs) for mu in mus]
    target = Poly(s_expr, *xs)
    allm = sorted(set().union(*[set(p.monoms()) for p in polys + [target]]), reverse=True)
    A = sp.Matrix([[P.coeff_monomial(m) for m in allm] for P in polys]).T
    b = sp.Matrix([target.coeff_monomial(m) for m in allm])
    sol, params = A.gauss_jordan_solve(b)
    out = {}
    for j, mu in enumerate(sorted(mus, key=lambda t: (-len(t), t))):   # 顯示順序: 單位元類在前
        out[mu] = sp.simplify(sol[mus.index(mu)] * z_mu(mu))   # chi^lam(mu)
    return out

for lam in [(1,1,1), (3,), (2,1)]:
    n = sum(lam)
    # Frobenius 公式 s_lam = sum_mu chi^lam(mu) z_mu^{-1} p_mu 需要恰有 n 個變數
    xs = sp.symbols(f'x1:{n+1}')
    print(f"chi^{lam}:", expand_schur_in_power(lam, xs))

chi^(1, 1, 1): {(1, 1, 1): 1, (2, 1): -1, (3,): 1}
chi^(3,): {(1, 1, 1): 1, (2, 1): 1, (3,): 1}
chi^(2, 1): {(1, 1, 1): 2, (2, 1): 0, (3,): -1}