判別式を根で表せたので、今度は基本対称式で表すことを考えましょう。3次方程式の場合基本対称式は根をp, q, rで表したとき、
s1 = p + q + r
s2 = pq + pr + qr
s3 = pqr
の3つがあります。判別式は6次なので、これらの基本対称式の積で6次になる式の線形結合になります。そのような6次の積は以下の7つがあります。
t1 = s16
t2 = s14s2
t3 = s13s3
t4 = s12s2
t5 = s1s2s3
t6 = s23
t7 = s32
これらを根で表します。
t1 = p6 + 6p5q + 6p5r + 15p4q2 + 30p4qr + 15p4r2 + 20p3q3 + 60p3q2r + 60p3qr2 + 20p3r3 + 15p2q4 + 60p2q3r + 90p2q2r2 + 60p2qr3 + 15p2r4 + 6pq5 + 30pq4r + 60pq3r2 + 60pq2r3 + 30pqr4 + 6pr5 + q6 + 6q5r + 15q4r2 + 20q3r3 + 15q2r4 + 6qr5 + r6
t2 = p5q + p5r + 4p4q2 + 9p4qr + 4p4r2 + 6p3q3 + 22p3q2r + 22p3qr2 + 6p3r3 + 4p2q4 + 22p2q3r + 36p2q2r2 + 22p2qr3 + 4p2r4 + pq5 + 9pq4r + 22pq3r2 + 22pq2r3 + 9pqr4 + pr5 + q5r + 4q4r2 + 6q3r3 + 4q2r4 + qr5
t3 = p4qr + 3p3q2r + 3p3qr2 + 3p2q3r + 6p2q2r2 + 3p2qr3 + pq4r + 3pq3r2 + 3pq2r3 + pqr4
t4 = p4q2 + 2p4qr + p4r2 + 2p3q3 + 8p3q2r + 8p3qr2 + 2p3r3 + p2q4 + 8p2q3r + 15p2q2r2 + 8p2qr3 + p2r4 + 2pq4r + 8pq3r2 + 8pq2r3 + 2pqr4 + q4r2 + 2q3r3 + q2r4
t5 = p3q2r + p3qr2 + p2q3r + 3p2q2r2 + p2qr3 + pq3r2 + pq2r3
t6 = p3q3 + 3p3q2r + 3p3qr2 + p3r3 + 3p2q3r + 6p2q2r2 + 3p2qr3 + 3pq3r2 + 3pq2r3 + q3r3
t7 = p2q2r2
これらで判別式を表せばよいわけです。なかなか大変そうですね。しかし、例えばp4q2が対称式の中にあれば、
p4r2 p2q4 q4r2 p2r4 q2r4
もあるはずです。だから、一つに代表させればよいわけです。p, q, rの指数が降順になる項で代表させましょう。代表する項は、
p6 p5q p4q2 p4qr p3q3 p3q2r p2q2r2
の7つあります。基本対称式の積も7つありましたね。これは偶然でしょうか。
ちょっと考えると、
s1 -> p
s2 -> pq
s3 -> pqr
と置き換えると対応がつくことがわかります。また、基本対称式の積は線形独立なのは明らかなので、これで対称式は基本対称式で表されることが証明されました。
基本対称式の積を代表項r1, ..., r7で表すと、
t1 = r1 + 6r2 + 15r3 + 30r4 + 20r5 + 60r6 + 90r7
t2 = r2 + 4r3 + 9r4 + 6r5 + 22r6 + 36r7
t3 = r4 + 3r6 + 6r7
t4 = r3 + 2r4 + 2r5 + 8r6 + 15r7
t5 = r6 + 3r7
t6 = r5 + 3r6 + 6r7
t7 = r7
これを解いて、
r1 = t1 + -6t2 + 6t3 + 9t4 + -12t5 + -2t6 + 3t7
r2 = t2 + -1t3 + -4t4 + 7t5 + 2t6 + -3t7
r3 = -2t3 + t4 + 4t5 + -2t6 + -3t7
r4 = t3 + -3t5 + 3t7
r5 = -3t5 + t6 + 3t7
r6 = t5 + -3t7
r7 = t7
これを判別式に代入すれば、判別式を3次方程式の係数で表すことができます。
-4a3c + a2b2 + 18abc - 4b3 - 27c2
4次は
-27a4d2 + 18a3bcd - 4a3c3 - 4a2b3d + a2b2c2 + 144a2bd2 - 6a2c2d - 80ab2cd + 18abc3 - 192acd2 + 16b4d - 4b3c2 - 128b2d2 + 144bc2d - 27c4 + 256d3
5次は
256a5e3 - 192a4bde2 - 128a4c2e2 + 144a4cd2e - 27a4d4 + 144a3b2ce2 - 6a3b2d2e - 80a3bc2de + 18a3bcd3 - 1600a3be3 + 16a3c4e - 4a3c3d2 + 160a3cde2 - 36a3d3e - 27a2b4e2 + 18a2b3cde - 4a2b3d3 - 4a2b2c3e + a2b2c2d2 + 1020a2b2de2 + 560a2bc2e2 - 746a2bcd2e + 144a2bd4 + 24a2c3de - 6a2c2d3 + 2000a2ce3 - 50a2d2e2 - 630ab3ce2 + 24ab3d2e + 356ab2c2de - 80ab2cd3 + 2250ab2e3 - 72abc4e + 18abc3d2 - 2050abcde2 + 160abd3e - 900ac3e2 + 1020ac2d2e - 192acd4 - 2500ade3 + 108b5e2 - 72b4cde + 16b4d3 + 16b3c3e - 4b3c2d2 - 900b3de2 + 825b2c2e2 + 560b2cd2e - 128b2d4 - 630bc3de + 144bc2d3 - 3750bce3 + 2000bd2e2 + 108c5e - 27c4d2 + 2250c2de2 - 1600cd3e + 256d5 + 3125e4
となりました。
from itertools import * from fractions import Fraction def inverse(A): def find_not_zero(m): for m1 in range(m + 1, n): if A[m1][m] != 0: return m1 def swap(m): m1 = find_not_zero(m) A[m], A[m1] = A[m1], A[m] B[m], B[m1] = B[m1], B[m] def div(row, d): for k in range(n): row[k] /= d def sub(row, row0, m): for k in range(n): row[k] -= row0[k] * m def is_all_int(B): return all(e.denominator == 1 for v in B for e in v) A = [ map(Fraction, v) for v in A ] n = len(A) B = [ [ 1 if i == j else 0 for j in range(n) ] for i in range(n) ] for m in range(n): if A[m][m] == 0: swap(m) a = A[m][m] div(A[m], a) div(B[m], a) for k in range(m + 1, n): a = A[k][m] sub(A[k], A[m], a) sub(B[k], B[m], a) for m in range(n - 2, -1, -1): for k in range(m + 1, n): sub(B[m], B[k], A[m][k]) if not is_all_int(B): print "error : B includes fractions." exit(1) return [ map(int, v) for v in B ] def mul(A, v): N = len(A) return [ sum(A[j][i] * v[j] for j in range(N)) for i in range(N) ] 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 __pow__(self, e): if e == 0: return cPolynomial(1) return self * self ** (e - 1) def add(self, term, coeff): if coeff == 0: return 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[1:]) 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 # x_k1..x_kn def mono2(ks): p = cPolynomial() a = list(ks) term = tuple(1 if a.count(i) == 1 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 def collect_reprs(n, m, upper = -1): if upper == -1: upper = m if n == 0: if m == 0: yield () return for k in xrange(min(m, upper), -1, -1): for t in collect_reprs(n - 1, m - k, k): yield (k,) + t def bases(N, M): def add(s1, s2): return tuple(e1 + e2 for e1, e2 in izip(s1, s2)) def mono(k): return tuple(int(l == k) for l in xrange(1, N + 1)) bs = [ [] ] for k in xrange(1, M + 1): b = set() if k <= N: b.add(mono(k)) for l in (l for l in xrange(1, N + 1) if l < k): b.update(add(s1, mono(l)) for s1 in bs[k-l]) bs.append(list(b)) return list(reversed(sorted(bs[-1]))) def sym2poly(b): return reduce(lambda x, (k, e): x * ss[k] ** e, enumerate(b), cPolynomial(1)) def encode(p): def is_descending(t): return all(k >= l for k, l in izip(t, t[1:])) a = [ 0 ] * L for t, c in p.iteritems(): if is_descending(t): a[rs[t]] = c return a def convert(p): for t, c in p.iteritems(): s = sum((k + 1) % 2 * e for k, e in enumerate(t)) if s % 2 == 1: p[t] = -c N = 5 M = N * (N - 1) # dimension ss = [ reduce(lambda p, ks: p + mono2(ks), combinations(range(N), n), cPolynomial()) for n in xrange(1, N + 1) ] rs = dict((r, k) for k, r in enumerate(collect_reprs(N, M))) L = len(rs) bs = bases(N, M) bs2 = map(sym2poly, bs) A = map(encode, bs2) D = determinant() v = encode(D) w = mul(inverse(A), v) p = cPolynomial() p.update((bs[k], w[k]) for k in range(L) if w[k] != 0) convert(p) print p
