n次方程式の判別式の具体的表示（1）

現代思想2011年4月号　特集=ガロアの思考　若き数学者の革命

この雑誌のどこだったかに、n次方程式の判別式は基本対称式で表される、と書いてありました。では、それはどのように表されるのでしょう。
判別式とは、根が重複しているかどうかが調べられる式です。2次方程式ax² + bx + cなら

D = b² - 4ac

が判別式ですね。ここからは、最高次数の係数は1として考えます。x² + ax + bの判別式は、

D = a² - 4b

です。さて、一般の場合はどう考えればよいでしょう。差積というものを考えます。要するに全ての差の積です。例えば、p, q, rの差積は

Δ(p, q, r) = (p - q)(p - r)(q - r)

です。このp, q, rを3次方程式の3つの根と考えると、根が重複すれば0、そうでなければそれ以外になります。ですからこれを判別式と呼んでもよいのですが、この式は少し性質が悪いのです。例えば、qとrを入れ替えてみましょう。そうすると差積の一つ目と2つ目の積は入れ替わるだけですが、3つ目の積は符号が変わり、結局全体として符号が変わってしまいます。好ましい性質とは入れ替えても式が変わらないことです。
入れ替えても式が変わらない式を対称式と呼びます。これは基本対称式で表されることが証明されています。基本対称式は、3変数なら3つあって、

p + q + r
pq + pr + qr
pqr

となります。これらは元の方程式の係数（あるいはその符号を変えたもの）になります。
ここで、差積の自乗を判別式としましょう。これは入れ替えても式が変わりません。そうすると、基本対称式で表される、すなわち元の方程式の係数で表されることになります。
まず、3次方程式でこの判別式を表してみましょう。

D = (p - q)²(p - r)²(q - r)²

を計算するのですが、ええっと、これは人間の計算するものではありませんね。オイラーならなんということもないのでしょうが。
しかし、我々にはプログラミングがあります。多項式クラスを作って計算しましょう。ここではp, q, rの代わりにa, b, cとします。c₁a^kb^lc^mという項は(k, l, m) => c₁という対応付けした辞書にします。この辞書を継承して、加減乗を定義します。これを使って計算すると、

a⁴b² + a⁴b² - 2a⁴bc + a⁴c² - 2a³b³ + 2a³b²c + 2a³bc² - 2a³c³ + a²b⁴ + 2a²b³c - 6a²b²c² + 2a²bc³ + a²c⁴ - 2ab⁴c + 2ab³c² + 2ab²c³- 2abc⁴ + b⁴c² - 2b³c³ + b²c⁴

3次でこれですか。大変ですね。19項あります。4次なら201項、5次なら2961項ありました。

from itertools import *

class cPolynomial(dict):
    def __init__(self, k = 0):
        if k != 0:
            self[(0,)*N] = k
    
    def __add__(self, a):
        return self.add_sub(True, a)
    
    def __sub__(self, a):
        return self.add_sub(False, a)
    
    def add_sub(self, isadd, a):
        b = self.copy()
        for term, coeff in a.iteritems():
            if isadd:
                b.add(term, coeff)
            else:
                b.add(term, -coeff)
        return b
    
    def __mul__(self, a):
        p = cPolynomial()
        for t1, c1 in self.iteritems():
            for t2, c2 in a.iteritems():
                t = tuple(e1 + e2 for e1, e2 in izip(t1, t2))
                p.add(t, c1 * c2)
        return p
    
    def add(self, term, coeff):
        if term in self:
            c = self.get(term)
            if c + coeff == 0:
                self.pop(term)
            else:
                self[term] += coeff
        else:
            self[term] = coeff
    
    def copy(self):
        b = cPolynomial()
        for key, value in self.iteritems():
            b[key] = value
        return b
    
    def __str__(self):
        items = self.items()
        items.sort()
        return self.toStr(list(reversed(items)))
    
    def toStr(self, terms):
        if len(terms) == 0:
            return "0"
        return self.strTerm(terms[0], True) \
                + "".join(self.strTerm(t, False) for t in terms)
    
    def strTerm(self, (term, coeff), ishead):
        return self.strCoeff(coeff, ishead) \
                + "".join(self.strTerm2(k, t) for k, t in enumerate(term))
    
    def strCoeff(self, c, ishead):
        if c == 1:
            return "" if ishead else " + "
        elif c == -1:
            return "-" if ishead else " - "
        elif c > 0:
            return str(c) if ishead else " + " + str(c)
        else:
            return str(c) if ishead else " - " + str(-c)
    
    def strTerm2(self, k, e):
        if e > 0:
            return chr(97 + k) + ("" if e == 1 else "^" + str(e))
        else:
            return ""

# x_k
def mono(k):
    p = cPolynomial()
    term = tuple(1 if k == i else 0 for i in range(N))
    p[term] = 1
    return p

def calc_diff_product():
    p = cPolynomial(1)
    for k, l in combinations(range(N), 2):
        p = p * (mono(k) - mono(l))
    return p

def determinant():
    p = calc_diff_product()
    return p * p

N = 3
print determinant()