about 5 years ago

Jason Liszka さんの Infinite lazy polynomials という記事を見て感動しました。
あ、遅延評価使えば無限次の多項式もプログラムで扱えるのか!とw
元記事ではScalaで実装してあったのですが浅学でよくわからなかったのでPythonで実装しなおして見ました。

無限次多項式クラス Poly

これから
$$
f(x)=\sum_{n=0}^{\infty}a_nx^n
$$
の形で書かれる無限次の多項式を考えます。詳しくは wiki:[冪級数]
例えば
$$
1+x+\frac{1}{2!}x^2+\frac{1}{3!}x^3+\frac{1}{4!}x^4+\cdots
$$
のようなものです。
$$
1+2x+x^2
$$
も\(x^3\)以降の項の係数が0だと考えれば無限次多項式に入ります。

Infinite lazy polynomialsではこのような無限次多項式をクラスとしてScala上で定義し種々の演算をできるようにしています。
面白かったので自分もこれに倣ってPython上でモジュールとして使ってみました。
以下では元記事に倣ってPythonのコードを元に解説してみます。

基本的な定義

まず無限次多項式を作るには係数\(a_n\)を与える必要があります。
この係数は無限個必要なので少し工夫が必要です。
そこで\(a_n\)を自然数から数への写像と考えて定義します。
例えば

x = Poly(lambda n: 1 if n==1 else 0)

こんなふうに定義をします。
Polyクラスはラムダ式を渡して初期化しており、このラムダ式が係数に当たるわけです。
この方法で
$$
\exp(x) = 1+x+\frac{1}{2!}x^2+\frac{1}{3!}x^3+\frac{1}{4!}x^4+\cdots
$$
を実装しようとすると以下のようになります。

from math import factorial
exp = Poly(lambda n: 1./factorial(n))

では実際にPolyクラスを定義していきましょう。

class Poly:
    def __init__(self,coeffs):
        self.coeffs = coeffs

    def __getitem__(self,n):
        return self.coeffs(n)

    def __call__(self,x,precision=10):
        return sum([ self[i] * x**i for i in xrange(precision) ])

まず初期化時にself.coeffsに係数を記憶しておき、x[n]のようなアクセスが来たら\(x^n\)番目の係数を返し、x(t)のようにアクセスされたら\(t\)における値を近似的に計算して返しています。
さらに

    def __repr__(self):
        return ", ".join([ str(self[i]) for i in xrange(10) ] + ["..."])

とすることでprint文で呼び出された時に\(x^9\)までの係数を","区切りで並べて表示するようにしています。

加算減算

無限次多項式の足し算引き算は以下のように定義されます。
$$
\sum_{n=0}^{\infty}a_nx^n + \sum_{n=0}^{\infty}b_nx^n = \sum_{n=0}^{\infty}(a_n+b_n)x^n \\
\sum_{n=0}^{\infty}a_nx^n - \sum_{n=0}^{\infty}b_nx^n = \sum_{n=0}^{\infty}(a_n-b_n)x^n
$$
これは以下のように実装しました。
Pythonは演算子のオーバーロードが使えるのでPolyクラスを直感的に使えるようになります。

    def __add__(self,other):
        return Poly(lambda n: self.coeffs(n) + other.coeffs(n))

    def __sub__(self,other):
        return Poly(lambda n: self.coeffs(n) - other.coeffs(n))

実際の動作です

プログラムで

$$
1+x
$$

などを計算するときにいちいち

Poly(lambda n: 1 if n==0 else 0) + x

と書くのは面倒なので 数値+Poly と計算された時は自動的に数値をPolyクラスに変換してくれるようにしました。

def isNumber(n):
    return isinstance(n,int) or isinstance(n,float) or isinstance(n,long) or isinstance(n,complex)

def numberToPoly(number):
    return Poly(lambda n: number if n==0 else 0)
...
    #In Poly Class Definition
    def __add__(self,other):
        if isNumber(other): other = numberToPoly(other)
        return Poly(lambda n: self.coeffs(n) + other.coeffs(n))


__radd____rsub__に関しても同様の実装をしています。

乗算

無限次多項式の乗算は以下のように定義されます
$$
\sum_{n=0}^{\infty}a_nx^n \times \sum_{n=0}^{\infty}b_nx^n = \sum_{n=0}^{\infty}\left(\sum_{k=0}^na_kb_{n-k}\right)x^n
$$
各項の係数は畳み込みになります。

これをプログラムで実装すると次のようになります

    def __mul__(self,other):
        return Poly(lambda n: sum([ self.coeffs(k) * other.coeffs(n-k) for k in xrange(n+1) ]))

実際の動作です

除算

無限次多項式による除算は無限次多項式の中だけで考えるのはある意味無意味なことかもしれません。
ですがここまで来たら実装したいということで無理矢理にでも実装してみましょう。
除算を行うには逆元を求めて掛ければいいのでまずは逆元について考えてみます
x が良い値をとれば
$$
\frac{1}{1-x} = 1+x+x^2+\cdots
$$
が成り立ちます。
一般の多項式\(q(x)\)に対しても
$$
\frac{1}{1-q(x)} = 1+q(x)+q(x)^2+\cdots
$$
が成り立つことが期待できます。
これを用いて一般の多項式\(p(x)\)に対して
$$
a = p(0), q(x) = 1 - \frac{p(x)}{a}
$$
とおけば
$$
\frac{1}{p(x)} = \frac{1}{a}\frac{1}{1-q(x)}=\frac{1}{a}\left(1+q(x)+q(x)^2+\cdots\right)
$$
と計算できます。
更に今\(q(0)=0\)であるので\(q(x)^n\)は\(x^n\)以上の項しか含んでいません。
これを利用して実装すると

    def inv(self):
        a = self[0]
        if a == 0 : return
        else:
            q = 1 - self / a
            return Poly(lambda n: sum([ (q**i)[n] for i in xrange(n+1) ]) / a)


となります。
これを用いれば除算は

    def __div__(self,other):
        return self * other.inv()

とかけますね。

べき乗

べき乗については掛け算を組み合わせることで実装できます。
実際のコードは以下のようになります。

    def __pow__(self,other):
        if other==0: return numberToPoly(1)
        p2 = self ** (other/2)
        if other%2==0: return p2 * p2
        else: return p2 * p2 * self

指数対数

指数関数、対数関数についてもpExppLogという形で実装しています。
実装方法は元記事に詳しく書いてあるので参照ください。

微分

微分に関しては元に記事には直接は書かれていないのですが実装してみました。
無限次多項式の微分は次のように定義されます
$$
\left(\sum_{n=0}^{\infty}a_nx^n\right)' = \sum_{n=0}^{\infty}a_{n+1}(n+1)x^n
$$
これを実装すると

    def derivative(self):
        return Poly(lambda n: self[n+1]*(n+1))

の様になります。簡単ですね!

グラフ

これも元の記事には書いてないのですがせっかく無限次多項式を作ったのだからグラフを見てみたいと思って実装しました。
描画には Tkinter を利用しています

def pGraph(poly,precision=10,x_min=-1.,x_max=1.,y_min=-1.,y_max=1.,scale=200.,step=0.01):
    from Tkinter import Tk,Canvas
    root = Tk()
    c = Canvas(root, width = scale*(x_max-x_min), height = scale*(y_max-y_min))
    c.pack()

    x = x_min
    y = poly(x,precision=precision)
    while x < x_max:
        x0 = x
        y0 = y
        x += step
        y = poly(x)

        xb = scale*(x0-x_min)
        yb = scale*(y_max-y_min)-scale*(y0-y_min)
        xa = scale*(x-x_min)
        ya = scale*(y_max-y_min)-scale*(y-y_min)
        c.create_line(xb,yb,xa,ya)

    root.mainloop()


モジュールとして使う

以上Pythonで作ってきましたが完成したコードはinfpoly.pyという名前でGistで公開しています
自由に使ってください
infpoly.py

ローカルに保存して

from infpoly import Poly

などのようにすれば使えると思います。

まとめ

Pythonで無限次多項式クラスを実装してみました。
実装したあとで友人に見せると
「それSICPに書いてあるよ」
と言われたのでいずれSICPも読みたいと思っています…

せっかく多項式のクラスを作ったので根探索や数値積分などを実装してもいいかもしれません
応用して有理関数クラスとかも考えられるかなぁ

← 追跡と迎撃のアルゴリズム 『21世紀の新しい数学』読みました →
 
comments powered by Disqus