何かを書き留める何か

数学や読んだ本について書く何かです。最近は社会人として生き残りの術を学ぶ日々です。

Pythonによる体上の多項式環

整数剰余環の次は体上の多項式環である。 j2kun/finite-fields · GitHubを参考に…どころか数学的な部分はまる写しという体たらくである。 メタクラスを使ったり、コメントをつけたり、unittestを書いたのが私の頑張った範囲である。

プログラミングによる有限体は整数の素数による剰余環か、体上の多項式環の既約多項式による剰余環を使うのがほとんどであろうか。 どちらも、「可換環の極大イデアルによる剰余環は体である」という事実(誰が最初に証明したのだろう?)を利用している。

import fractions
import itertools


def polynomials_over(field=fractions.Fraction):
    """
    体上の多項式環を生成する。

    デフォルト引数は有理数体である。
    """
    class _Polynomial(type):

        """ 係数体を保持するためのメタクラス """
        @classmethod
        def __prepare__(metacls, name, bases, **kwags):
            return super().__prepare__(metacls, name, bases)

        def __new__(metacls, name, bases, ns, **kwags):
            return super().__new__(metacls, name, bases, ns)

        def __init__(cls, name, bases, ns, **kwags):
            type.__init__(cls, name, bases, ns)
            cls.fields = kwags['fields']

    class Polynomial(metaclass=_Polynomial, fields=field):

        """ 体 field を係数とする多項式環 """

        def __init__(self, c):
            if type(c) is Polynomial:
                self.coefficients = c.coefficients
            elif isinstance(c, Polynomial.fields):
                self.coefficients = [c]
            elif not hasattr(c, '__iter__') and not hasattr(c, 'iter'):
                self.coefficients = [Polynomial.fields(c)]
            else:
                self.coefficients = c

            if len(c) == 0:
                self.coefficients = c
            i = len(c) - 1
            while i >= 0 and c[i] == Polynomial.fields(0):
                i -= 1
            else:
                self.coefficients = c[:i + 1]

        @classmethod
        def factory(cls, coef):
            """ 係数リストの値を係数体の値に変換する """
            return Polynomial([cls.fields(c) for c in coef])

        def is_zero(self):
            """ 多項式が 0 であるかを判定する """
            return self.coefficients == []

        def __repr__(self):
            """ 多項式を出力する """
            if self.is_zero():
                return '0'

            return ' + '.join(['{} x^{}'.format(a, i)
                               if i > 0 else '{}'.format(a)
                               for i, a in enumerate(self.coefficients)])

        def __len__(self):
            """ 多項式の項数を返す """
            return len(self.coefficients)

        def __iter__(self):
            """ 多項式の各項の係数を辿るイテレータ """
            return iter(self.coefficients)

        def __neg__(self):
            """ 多項式の符号を反転させる """
            return Polynomial([-c for c in self])

        def __sub__(self, other):
            """ 多項式の減法 """
            return self + (-other)

        def leading_coef(self):
            """ 多項式の最高次数の係数 """
            return self.coefficients[-1]

        def degree(self):
            """ 多項式の次数を返す """
            return len(self) - 1

        def abs(self):
            """ 
            多項式の次数を返す

            euclid.pyとの整合性を取るために存在する
            """
            return self.degree()

        def __eq__(self, other):
            """ 多項式の等号 """
            return (self.degree() == other.degree()
                    and all([x == y for (x, y) in zip(self, other)]))

        def __add__(self, other):
            """ 多項式の加法 """
            new_coef = [sum(x, Polynomial.fields(0)) for x in itertools.zip_longest(
                self, other, fillvalue=Polynomial.fields(0))]
            return Polynomial(new_coef)

        def __mul__(self, other):
            """ 
            多項式の乗法

            素朴な実装であるので、多項式f, gの積のオーダーは
            O(def(f)*def(g))となる。
            """
            if self.is_zero() or other.is_zero():
                return Polynomial([])

            new_coef = [Polynomial.fields(0)
                        for _ in range(len(self) + len(other) - 1)]

            for i, a in enumerate(self):
                for j, b in enumerate(other):
                    new_coef[i + j] += a * b

            return Polynomial(new_coef)

        def __divmod__(self, divisor):
            """ 
            多項式の除法

            素朴な実装である。
            """
            quotient, remainder = Polynomial([]), self
            divisor_deg = divisor.degree()
            divisor_lc = divisor.leading_coef()

            while remainder.degree() >= divisor_deg:
                monomial_exponent = remainder.degree() - divisor_deg
                monomial_zeros = [Polynomial.fields(0)
                                  for _ in range(monomial_exponent)]
                monomial_divisor = Polynomial(
                    monomial_zeros + [remainder.leading_coef() / divisor_lc])

                quotient += monomial_divisor
                remainder -= monomial_divisor * divisor

            return quotient, remainder

        def __truediv__(self, divisor):
            """ 多項式の除法の商を返す """
            if divisor.is_zero():
                raise ZeroDivisionError
            return divmod(self, divisor)[0]

        def __mod__(self, divisor):
            """ 多項式の除法の剰余を返す """
            if divisor.is_zero():
                raise ZeroDivisionError
            return divmod(self, divisor)[1]

    return Polynomial

ユニットテストはこちら。

from fractions import Fraction
import unittest

import polynomials


class TestPolynomialsOverQ(unittest.TestCase):

    def setUp(self):
        self.p = polynomials.polynomials_over(Fraction).factory

    def test_factory(self):
        self.assertEqual(self.p([]), self.p([]))
        self.assertEqual(self.p([1, 2, 0]), self.p([1, 2, 0, 0]))

    def test_add(self):
        self.assertEqual(self.p([1, 2, 3]), self.p([1, 0, 3]) + self.p([0, 2]))
        self.assertEqual(self.p([1, 2, 3]), self.p([1, 2, 3]) + self.p([]))
        self.assertEqual(self.p([5, 2, 3]), self.p([4]) + self.p([1, 2, 3]))
        self.assertEqual(
            self.p([1, 2]), self.p([1, 2, 3]) + self.p([0, 0, -3]))

    def test_sub(self):
        self.assertEqual(
            self.p([1, -2, 3]), self.p([1, 0, 3]) - self.p([0, 2]))
        self.assertEqual(self.p([1, 2, 3]), self.p([1, 2, 3]) - self.p([]))
        self.assertEqual(self.p([-1, -2, -3]), self.p([]) - self.p([1, 2, 3]))

    def test_multiple(self):
        self.assertEqual(self.p([1, 2, 1]), self.p([1, 1]) * self.p([1, 1]))
        self.assertEqual(
            self.p([2, 5, 5, 3]), self.p([2, 3]) * self.p([1, 1, 1]))
        self.assertEqual(self.p([0, 7, 49]), self.p([0, 1, 7]) * self.p([7]))

    def test_division(self):
        self.assertEqual(self.p([1, 1, 1, 1, 1, 1]), self.p(
            [-1, 0, 0, 0, 0, 0, 1]) / self.p([-1, 1]))
        self.assertEqual(
            self.p([-1, 1, -1, 1, -1, 1]), self.p([1, 0, 0, 0, 0, 0, 1]) / self.p([1, 1]))
        self.assertEqual(self.p([]), self.p([]) / self.p([1, 1]))
        self.assertEqual(self.p([1, 1]), self.p([1, 1]) / self.p([1]))
        self.assertEqual(self.p([1, 1]), self.p([2, 2]) / self.p([2]))

    def test_modulus(self):
        self.assertEqual(self.p([]), self.p([1, 7, 49]) % self.p([7]))
        self.assertEqual(
            self.p([-7]), self.p([-3, 10, -5, 3]) % self.p([1, 3]))

if __name__ == '__main__':
    unittest.main()