Pythonの組み込み函数sumについて

Pythonで有限体を実装しようと思い立ち、j2kun/finite-fields · GitHubを参考に体上の多項式環を実装しようとしていた。 多項式の和で組み込みのsumを利用しようとしたらTypeErrorが生じた。原因を調べたところ、sum内で組み込み型intと自作の体のクラスと__add__をしようとしてはまっていることがわかった。自作の体のクラスにはきちんと__add__を定義しているのに…と思ってsumの実装を調べてみた。

組み込みのsumのコードの本質的と思われる箇所の抜粋は以下の通りである。

static PyObject *
builtin_sum_impl(PyModuleDef *module, PyObject *iterable, PyObject *start)
{
    PyObject *result = start;
    PyObject *temp, *item, *iter;

    iter = PyObject_GetIter(iterable);
    if (iter == NULL)
        return NULL;

    if (result == NULL) {
        result = PyLong_FromLong(0);
        if (result == NULL) {
            Py_DECREF(iter);
            return NULL;
        }
    } else {
        /* reject string values for 'start' parameter */
        if (PyUnicode_Check(result)) {
            PyErr_SetString(PyExc_TypeError,
                "sum() can't sum strings [use ''.join(seq) instead]");
            Py_DECREF(iter);
            return NULL;
        }
        if (PyBytes_Check(result)) {
            PyErr_SetString(PyExc_TypeError,
                "sum() can't sum bytes [use b''.join(seq) instead]");
            Py_DECREF(iter);
            return NULL;
        }
        if (PyByteArray_Check(result)) {
            PyErr_SetString(PyExc_TypeError,
                "sum() can't sum bytearray [use b''.join(seq) instead]");
            Py_DECREF(iter);
            return NULL;
        }
        Py_INCREF(result);
    }

    for(;;) {
        item = PyIter_Next(iter);
        if (item == NULL) {
            /* error, or end-of-sequence */
            if (PyErr_Occurred()) {
                Py_DECREF(result);
                result = NULL;
            }
            break;
        }

        temp = PyNumber_Add(result, item);
        Py_DECREF(result);
        Py_DECREF(item);
        result = temp;
        if (result == NULL)
            break;
    }
    Py_DECREF(iter);
    return result;
}

わかったこととして、

    if (result == NULL) {
        result = PyLong_FromLong(0);
        if (result == NULL) {
            Py_DECREF(iter);
            return NULL;
        }

の部分でsum(iterable, start)のstartに何も入れないとint型の0がstartに入り、 自作の型と__add__をしようとしてエラーになった、と。 そのため、startに自作の有限体の零元を入れることで期待通りの動作をした。

また、sum関数は文字列を拒絶する旨がヘルプに書いてあるが、引用したコードを見るとバイト列も拒絶していることがわかる。 普段はPythonの内部ソースなんて気にしていないが、こう読んでみると面白い。

なお、ソースコードにはArgument Clinic How-To — Python 3.4.3 ドキュメントにあるツールが使われた痕跡があった。 Argument Clinicは"Monty Python's Flying Circus"の第3シーズンエピソード3の最後のスケッチが由来だと思われる。 "I'd like to have an argument please." という不思議な台詞から始まる。 案内された部屋に入ると罵声を浴びせられたり、謎の水掛け論の応酬が始まったり、仕舞いには突然ハンマーで殴られて「ここは頭を殴られる教室だ」と続き、そして警察が…とチャップマンとクリーズのケンブリッジ組が書いたと思われる会話主体の面白いスケッチである。

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()