yukiのブログ

作ったものなど

私のemacs環境をコピーする方法

emacsの設定をあれこれいじった結果それなりに使えるものになったので、emacsがインストールされていないPCで同じ設定のemacsを利用する方法をメモしておきます。未来の自分用です。

環境

Ubuntu 14.04 LTS

できるようになること

環境をコピーして拡張機能をすべて使えるようにすれば、C, OCaml, Verilog HDLでシンタックスハイライト、補完、リアルタイムなエラーチェックができるようになります。

コピーする手順

順に説明します。

emacsの最新版をインストー

Ubuntu 14.04 LTSでは普通にapt-getしても最新版のEmacsが入らないので、手動でppaを追加する必要があります。

sudo apt-add-repository -y ppa:adrozdoff/emacs
sudo apt update
sudo apt install emacs25

~/.emacsと~/.emacs.elを削除

そのままです。もしかしたらどちらも作られて無いかも。そうだったら何もしないでいいです。

リポジトリをクローンして、~/.emacs.d/におく

リポジトリここ。もともと存在するファイルは全て削除してください。

Caskをインストールする

Caskはemacs用のパッケージマネージャです。emacsとは別にインストールする必要があります。

インストール方法は公式サイトを参照。現在のバージョンでは以下のコマンドです。

curl -fsSL https://raw.githubusercontent.com/cask/cask/master/go | python

~/.emacs.d で cask install を実行

~/.emacs.d/Cask に使用するパッケージが羅列されており、cask install を実行することでそれらのパッケージがすべてインストールされます。

通常のインストールではcaskの実行ファイルにPATHが通らないっぽいので、その場合は

~/.emacs.d$ ~/.cask/bin/cask install

としてください。

emacsを起動して、エラーを潰していく

cask install が終わったあとにemacsを起動すれば基本的な拡張機能が読み込まれますが、一部は正常に動きません。一部の拡張機能には追加のパッケージが必要なためです(clang-formatterとかironyとかGNU Globalとか)。

具体的にどの拡張機能にどのパッケージが必要なのかは忘れたので、エラーを見てその都度パッケージ名で検索してください。そんなにアクロバティックなことはやってないので、ググって公式サイトのやり方に従えば大丈夫なはずです。

拡張機能の読み込みに関するエラーは*init log*タブか*minibuffer*タブに表示されます。Ctrl-x + Ctrl-b でタブの一覧を表示してから見に行ってください。

emacs拡張機能に必要となるパッケージは基本的にはapt-getで入るはずですが、Ubuntuのバージョンが最新でないために面倒な手順が必要になることがあります。特に GNU Global はソースからビルドする必要があるはずです。

多分この部分が一番めんどくさいでしょう。

(必要なら)フォントをインストールする

フォントにはrictyを使っています。インストールするには公式サイトからダウンロードして自分でビルドする必要があります。

別のフォントを使う場合には、~/.emacs.d/inits/98_themes.elを編集してください。

Sympyで方程式を解いてみる

Sympyで方程式を解いてみる

Sympyがまだ面白そうだったのでもっと色々やりました。

とりあえず前に受けた微分方程式講義ノートを引っ張り出してそこに書いてある練習問題を解いていきました。

微分方程式はすぐ解ける形もあればあれこれ手を焼いてやらないと解けないものもあるようです。すぐに解けないときは小手先の変更(hintを変えるとか)だけではうまくいかないことが多いので、さっさと地道な方法で解いたほうがいいと思います。

この記事はjupyter notebookで書いたのをmarkdownに変換したものです。notebookはここにあります。

# 初期処理
%matplotlib inline
from sympy import *
from sympy.plotting import plot
from IPython.display import display
from  sympy.solvers.ode import *
init_printing()

# デフォルトだと微分方程式を解いたときに誤差関数(error function)を使うのでそれを書き直す関数
def expand_erfs(expr):
    """ 誤差関数を開く """
    return expr.replace(erf, lambda x: 2/sqrt(pi) * Integral(exp(-t**2), (t, 0, x)))\
               .replace(erfi, lambda x: 2/sqrt(pi) * Integral(exp(t**2), (t, 0, x))).simplify()

# 積分定数を潰す関数constantsimp()は積分定数を引数として
# 渡さなければいけないので、C1,C2,...をシンボルとして定義しなければいけない。
# それだとめんどくさいので、積分定数でないシンボルを渡せばすむ関数があると便利。
def constsimp(expr, variable):
    """ variable以外の自由変数をconstantsimp()に渡す """
    return constantsimp(expr, [c for c in ans.free_symbols if c != variable])

x = symbols("x")
a, b, c, t = symbols("a b c t")
y, f, g, P, Q, R = symbols('y f g P Q R', cls=Function)

方程式

expr = Eq(x**2+3*x-4)
display(expr, solveset(expr, x)) # solve()という関数もあるがsolveset()が推奨されている

\(\displaystyle x^{2} + 3 x - 4 = 0\)

\(\displaystyle \left\{-4, 1\right\}\)

expr = x**3+3*x**2+3*x-4
display(Eq(expr), solveset(expr))

\(\displaystyle x^{3} + 3 x^{2} + 3 x - 4 = 0\)

\(\displaystyle \left\{-1 + \sqrt[3]{5}, -1 - \frac{\sqrt[3]{5}}{2} - \frac{\sqrt{3} i}{2} \sqrt[3]{5}, -1 - \frac{\sqrt[3]{5}}{2} + \frac{\sqrt{3} i}{2} \sqrt[3]{5}\right\}\)

a, b, c = symbols("a b c")
solve(a*x**2+b*x+c, x)

\(\displaystyle \left [ \frac{1}{2 a} \left(- b + \sqrt{- 4 a c + b^{2}}\right), \quad - \frac{1}{2 a} \left(b + \sqrt{- 4 a c + b^{2}}\right)\right ]\)

solveset(Gt(x**2, 0), domain=S.Reals)

\(\displaystyle \left(-\infty, 0\right) \cup \left(0, \infty\right)\)

solveset(sin(x) - sqrt(3)/2, domain=S.Reals)

\(\displaystyle \left\{2 n \pi + \frac{2 \pi}{3}\; \middle|\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \frac{\pi}{3}\; \middle|\; n \in \mathbb{Z}\right\}\)

微分方程式

変数分離型

eq = Eq(Derivative(y(x)), g(x)*g(y(x)))
display(eq, dsolve(eq, y(x)))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = g{\left (x \right )} g{\left (y{\left (x \right )} \right )}\)

\(\displaystyle \int^{y{\left (x \right )}} \frac{1}{g{\left (y \right )}}\, dy = C_{1} + \int g{\left (x \right )}\, dx\)

例題: ロジスティック方程式

# シンボルを定義
t, r, K = symbols("t r K")
N = symbols("N", cls=Function)
# ロジスティック方程式を定義
eq = Eq(Derivative(N(t)), r*(1-(N(t)/K))*N(t))
display(eq)
# 一般解を求める
ans = dsolve(eq)
display(ans)

# 特殊化する
args = [(K, 100), (r, 0.75)]
init_val = 1.0
func = ans.rhs.subs(args)
# 初期値を方程式の解として求める
c = solve(Eq(func.subs([(t, 0)]), init_val), domain=S.Reals)
# プロットする
C1 = symbols("C1") # 積分定数用のシンボルを定義しなきゃいけないのがめんどくさい
plot(func.subs([(C1, c[0])]), (t, -10, 15))

\(\displaystyle \frac{d}{d t} N{\left (t \right )} = r \left(1 - \frac{1}{K} N{\left (t \right )}\right) N{\left (t \right )}\)

\(\displaystyle N{\left (t \right )} = \frac{K e^{C_{1} K + r t}}{e^{C_{1} K + r t} - 1}\)

png

<sympy.plotting.plot.Plot at 0x23613929b70>

LiouVille(リウビル)型

LiouVille = Eq(y(x).diff(x,2) + P(x)*y(x).diff(x) + Q(y(x))*y(x).diff(x)**2)
display(LiouVille, dsolve(LiouVille, y(x), hint="Liouville"))

\(\displaystyle P{\left (x \right )} \frac{d}{d x} y{\left (x \right )} + Q{\left (y{\left (x \right )} \right )} \frac{d}{d x} y{\left (x \right )}^{2} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle C_{1} + C_{2} \int e^{- \int P{\left (x \right )}\, dx}\, dx + \int^{y{\left (x \right )}} e^{\int Q{\left (y \right )}\, dy}\, dy = 0\)

例題

eq = Eq(y(x).diff(x,2) + 2*x*y(x).diff(x) + 2*y(x)*y(x).diff(x) ** 2)
ans = dsolve(eq, hint="Liouville")
display(eq, ans, constsimp(expand_erfs(ans), x))

\(\displaystyle 2 x \frac{d}{d x} y{\left (x \right )} + 2 y{\left (x \right )} \frac{d}{d x} y{\left (x \right )}^{2} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle C_{1} + C_{2} \operatorname{erf}{\left (x \right )} + \frac{\sqrt{\pi}}{2} \operatorname{erfi}{\left (y{\left (x \right )} \right )} = 0\)

\(\displaystyle C_{1} + C_{2} \int_{0}^{x} e^{- t^{2}}\, dt + \int_{0}^{y{\left (x \right )}} e^{t^{2}}\, dt = 0\)

eq = Eq(y(x).diff(x), sin(x)*tan(y(x)))
display(eq, dsolve(eq, y(x)))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \sin{\left (x \right )} \tan{\left (y{\left (x \right )} \right )}\)

\(\displaystyle \left [ y{\left (x \right )} = - \operatorname{a{}c{}o{}s}{\left (- \sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )} + 2 \pi, \quad y{\left (x \right )} = - \operatorname{a{}c{}o{}s}{\left (\sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )} + 2 \pi, \quad y{\left (x \right )} = \operatorname{a{}c{}o{}s}{\left (- \sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )}, \quad y{\left (x \right )} = \operatorname{a{}c{}o{}s}{\left (\sqrt{C_{1} e^{- 2 \cos{\left (x \right )}} + 1} \right )}\right ]\)

同次型

eq = Eq(y(x).diff(x), f(y(x)/x))
display(eq, dsolve(eq, y(x)))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = f{\left (\frac{1}{x} y{\left (x \right )} \right )}\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- \int^{\frac{x}{y{\left (x \right )}}} \frac{f{\left (\frac{1}{u_{2}} \right )}}{u_{2} f{\left (\frac{1}{u_{2}} \right )} - 1}\, du_{2}}\)

eq = Eq(y(x).diff(x) , exp(y(x)/x) + y(x)/x)
display(eq, dsolve(eq))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = e^{\frac{1}{x} y{\left (x \right )}} + \frac{1}{x} y{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \log{\left (\left(\frac{1}{C_{1} - \log{\left (x \right )}}\right)^{x} \right )}\)

分数型

eq = Eq(y(x).diff(x), (x-2*y(x)+3)/(2*x+y(x)-4))
display(eq, dsolve(eq))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \frac{x - 2 y{\left (x \right )} + 3}{2 x + y{\left (x \right )} - 4}\)

\(\displaystyle \left [ y{\left (x \right )} = - 2 x - \sqrt{C_{1} + 5 x^{2} - 10 x} + 4, \quad y{\left (x \right )} = - 2 x + \sqrt{C_{1} + 5 x^{2} - 10 x} + 4\right ]\)

eq = Eq(y(x).diff(x), x*y(x)/(x**2+y(x)))
display(eq, dsolve(eq, hint="lie_group"))

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \frac{x y{\left (x \right )}}{x^{2} + y{\left (x \right )}}\)

\(\displaystyle \left [ y{\left (x \right )} = C_{1} \left(- \sqrt{- 2 C_{1} x^{2} + 1} - 1\right), \quad y{\left (x \right )} = C_{1} \left(\sqrt{- 2 C_{1} x^{2} + 1} - 1\right)\right ]\)

# 自動で解けないのもあるので、そのときは答を代入する
eq = Eq(y(x).diff(x), x**2*y(x)/(x**3+y(x)))
# この答えは手で解いたもの
ans_implicit = Eq(y(x)**3/(2*x**3+3*y(x)), 6*C1)
# 実数解だけ取る(虚数解でも成り立つだろうけどとりあえず無視)
ans = solve(ans_implicit, y(x))[0]
# これがTrueなのでans_implicitが正しいと分かる
display(eq, ans_implicit, Eq(y(x), ans), eq.subs([(y(x), ans)]).doit().simplify())

\(\displaystyle \frac{d}{d x} y{\left (x \right )} = \frac{x^{2} y{\left (x \right )}}{x^{3} + y{\left (x \right )}}\)

\(\displaystyle \frac{y^{3}{\left (x \right )}}{2 x^{3} + 3 y{\left (x \right )}} = 6 C_{1}\)

\(\displaystyle y{\left (x \right )} = - \frac{18 C_{1}}{\sqrt[3]{- 162 C_{1} x^{3} + \frac{1}{2} \sqrt{- 629856 C_{1}^{3} + 104976 C_{1}^{2} x^{6}}}} - \frac{1}{3} \sqrt[3]{- 162 C_{1} x^{3} + \frac{1}{2} \sqrt{- 629856 C_{1}^{3} + 104976 C_{1}^{2} x^{6}}}\)

\(\displaystyle \mathrm{True}\)

定数変化法

eq = Eq(y(x).diff(x) + P(x)*y(x), Q(x))
display(eq, dsolve(eq, y(x), hint="almost_linear"))

\(\displaystyle P{\left (x \right )} y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} = Q{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + \int Q{\left (x \right )} e^{\int P{\left (x \right )}\, dx}\, dx\right) e^{- \int P{\left (x \right )}\, dx}\)

eq = Eq(x*y(x).diff(x) - 2*y(x), x**3*cos(x))
display(eq, dsolve(eq, y(x)))

\(\displaystyle x \frac{d}{d x} y{\left (x \right )} - 2 y{\left (x \right )} = x^{3} \cos{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = x^{2} \left(C_{1} + \sin{\left (x \right )}\right)\)

# 一般解から解いてみる
eq = Eq(y(x).diff(x) + P(x)*y(x), Q(x))
display(eq.subs([(P(x), -2/x), (Q(x), x**2*cos(x))]),
        dsolve(eq, hint="almost_linear").subs([(P(x), -2/x), (Q(x), x**2*cos(x))]).doit())

\(\displaystyle \frac{d}{d x} y{\left (x \right )} - \frac{2}{x} y{\left (x \right )} = x^{2} \cos{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = x^{2} \left(C_{1} + \sin{\left (x \right )}\right)\)

Bernoulli(ベルヌーイ)型

n = symbols("n")
eq = Eq(y(x).diff(x) + P(x)*y(x), Q(x)*y(x)**n)
display(eq, dsolve(eq, hint="Bernoulli"))

\(\displaystyle P{\left (x \right )} y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} = Q{\left (x \right )} y^{n}{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \left(\left(C_{1} - \left(n - 1\right) \int Q{\left (x \right )} \left(e^{- n \int P{\left (x \right )}\, dx}\right) e^{\int P{\left (x \right )}\, dx}\, dx\right) e^{- \left(- n + 1\right) \int P{\left (x \right )}\, dx}\right)^{\frac{1}{- n + 1}}\)

eq = Eq(y(x).diff(x)+x*y(x), x*y(x)**3)
display(eq, dsolve(eq, hint="Bernoulli"))

\(\displaystyle x y{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} = x y^{3}{\left (x \right )}\)

\(\displaystyle y{\left (x \right )} = \frac{1}{\sqrt{\left(C_{1} + e^{- x^{2}}\right) e^{x^{2}}}}\)

Riccati(リッカチ)型

# 自動で解くのは無理っぽい
eq = Eq(y(x).diff(x) - y(x)**2 - (2-x)*y(x) + 2*x - 1)
spec_ans = x
display(eq.subs([(y(x), spec_ans)]).doit().simplify())
bern = eq.subs([(y(x), P(x)+spec_ans)]).doit().simplify()
ans = dsolve(bern, hint="Bernoulli")
display(eq, bern, ans, ans.subs([(P(x), y(x)-x)]))

\(\displaystyle \mathrm{True}\)

\(\displaystyle 2 x - \left(- x + 2\right) y{\left (x \right )} - y^{2}{\left (x \right )} + \frac{d}{d x} y{\left (x \right )} - 1 = 0\)

\(\displaystyle - x P{\left (x \right )} - P^{2}{\left (x \right )} - 2 P{\left (x \right )} + \frac{d}{d x} P{\left (x \right )} = 0\)

\(\displaystyle P{\left (x \right )} = \frac{e^{x \left(\frac{x}{2} + 2\right)}}{C_{1} - \frac{\sqrt{2} \sqrt{\pi}}{2 e^{2}} \operatorname{erfi}{\left (\sqrt{2} \left(\frac{x}{2} + 1\right) \right )}}\)

\(\displaystyle - x + y{\left (x \right )} = \frac{e^{x \left(\frac{x}{2} + 2\right)}}{C_{1} - \frac{\sqrt{2} \sqrt{\pi}}{2 e^{2}} \operatorname{erfi}{\left (\sqrt{2} \left(\frac{x}{2} + 1\right) \right )}}\)

Clairaur(クレロー)型

eq = Eq(y(x), x*y(x).diff(x) - y(x).diff(x)**2)
dsolve(eq)

\(\displaystyle y{\left (x \right )} = \frac{x^{2}}{4} - \frac{1}{4} \left(C_{1} + x\right)^{2}\)

減衰振動

m, k, u = symbols("m k u")
eq = Eq(m*y(x).diff(x,2), -k*y(x) - u * y(x).diff(x))
dsolve(eq)

\(\displaystyle y{\left (x \right )} = C_{1} e^{\frac{x}{2 m} \left(- u - \sqrt{- 4 k m + u^{2}}\right)} + C_{2} e^{\frac{x}{2 m} \left(- u + \sqrt{- 4 k m + u^{2}}\right)}\)

定数係数斉次線形

# 一番簡単なやつ
eq = Eq(y(x).diff(x,2) - y(x).diff(x) - 6*y(x))
display(eq, dsolve(eq))

\(\displaystyle - 6 y{\left (x \right )} - \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{3 x}\)

# 重解のパターン
eq = Eq(y(x).diff(x, 2) - 6*y(x).diff(x) + 9*y(x))
display(eq, dsolve(eq))

\(\displaystyle 9 y{\left (x \right )} - 6 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{3 x}\)

eq = Eq(y(x).diff(x,3) - 5*y(x).diff(x,2) + 8*y(x).diff(x) - 4*y(x))
display(eq, dsolve(eq))

\(\displaystyle - 4 y{\left (x \right )} + 8 \frac{d}{d x} y{\left (x \right )} - 5 \frac{d^{2}}{d x^{2}} y{\left (x \right )} + \frac{d^{3}}{d x^{3}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + \left(C_{2} + C_{3} x\right) e^{x}\right) e^{x}\)

eq = Eq(y(x).diff(x,4) - 8*y(x).diff(x,2) + 16*y(x))
display(eq, dsolve(eq))

\(\displaystyle 16 y{\left (x \right )} - 8 \frac{d^{2}}{d x^{2}} y{\left (x \right )} + \frac{d^{4}}{d x^{4}} y{\left (x \right )} = 0\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{- 2 x} + \left(C_{3} + C_{4} x\right) e^{2 x}\)

定数係数非斉次線形

eq = Eq(y(x).diff(x, 2) + 3*y(x).diff(x) + 2*y(x), exp(3*x))
display(eq, dsolve(eq))

\(\displaystyle 2 y{\left (x \right )} + 3 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = e^{3 x}\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{- x} + \frac{e^{3 x}}{20}\)

eq = Eq(y(x).diff(x, 2) + 2*y(x).diff(x) + y(x), exp(-x))
display(eq, dsolve(eq))

\(\displaystyle y{\left (x \right )} + 2 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = e^{- x}\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x + \frac{x^{2}}{2}\right) e^{- x}\)

eq = Eq(y(x).diff(x, 2) + 3*y(x).diff(x) + 2*y(x), x)
display(eq, dsolve(eq))

\(\displaystyle 2 y{\left (x \right )} + 3 \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = x\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{- x} + \frac{x}{2} - \frac{3}{4}\)

eq = Eq(y(x).diff(x, 2) - y(x).diff(x) - 2*y(x), x**2)
display(eq, dsolve(eq))

\(\displaystyle - 2 y{\left (x \right )} - \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = x^{2}\)

\(\displaystyle y{\left (x \right )} = C_{1} e^{- x} + C_{2} e^{2 x} - \frac{x^{2}}{2} + \frac{x}{2} - \frac{3}{4}\)

eq = Eq(y(x).diff(x, 2) + 2*a*y(x).diff(x) + a**2*y(x), x)
display(eq, dsolve(eq))

\(\displaystyle a^{2} y{\left (x \right )} + 2 a \frac{d}{d x} y{\left (x \right )} + \frac{d^{2}}{d x^{2}} y{\left (x \right )} = x\)

\(\displaystyle y{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{- a x} + \frac{x}{a^{2}} - \frac{2}{a^{3}}\)

Sympyを試す

 

Sympyを試す

Sympyが面白かったので色々試しました。

このページはjupyter notebookを作ったものをmarkdownに変換したものです。もととなったnotebookはここにあります。

初期設定

# ここはとりあえずコピペすればよさそう
%matplotlib inline
from sympy import *
from sympy.plotting import plot
from IPython.display import display
init_printing()

シンボルの定義

x = Symbol("x")
y = Symbol("y")            # これでもいいけど...
x, y, z = symbols("x y z") # こっちのほうが楽なのでこっちでいいと思う
f, g = symbols('f g', cls=Function)

# displayはnotebook側の関数で、引数を良い感じに出力してくれる。
display(x, y, z)

\(\displaystyle x\)

\(\displaystyle y\)

\(\displaystyle z\)

式の定義

expr = (x + y) ** 5

# display()を使うとnotebookを通じて出力できる
display(expr)
display(x * expr * y) # 変数の並び順は簡略化されるらしい

\(\displaystyle \left(x + y\right)^{5}\)

\(\displaystyle x y \left(x + y\right)^{5}\)

式の展開/因数分解

expr = (x + y) ** 5

display(expr.expand())                                # expand()で式を展開できる
display((x**3 + 3*x**2*y + 3*x*y**2 + y**3).factor()) # factor()はexpand()の逆

\(\displaystyle x^{5} + 5 x^{4} y + 10 x^{3} y^{2} + 10 x^{2} y^{3} + 5 x y^{4} + y^{5}\)

\(\displaystyle \left(x + y\right)^{3}\)

式の比較

expr = (x + y) ** 2
expanded = expand(expr)
display(expr)
display(expanded)

\(\displaystyle \left(x + y\right)^{2}\)

\(\displaystyle x^{2} + 2 x y + y^{2}\)

display(expr == expanded)              # この比較は代数的な比較ではない
display(Eq(expr, expanded))            # Eq()を使えば代数的な比較を作れる
display(Eq(expr, expanded).simplify()) # さらにSimplify()で恒等的に成り立つかどうかをチェックできる
False

\(\displaystyle \left(x + y\right)^{2} = x^{2} + 2 x y + y^{2}\)

\(\displaystyle \mathrm{True}\)

expr.equals(expanded) # equals()はランダムな点でテストをして等しいかどうかを判定するらしい
True

代入

expr = sin(x) + 1
expr.subs(x, pi / 3)

\(\displaystyle \frac{\sqrt{3}}{2} + 1\)

expr = x**2 + 3*x + 4
expr.subs(x, expr) # 式の代入も可能

\(\displaystyle 3 x^{2} + 9 x + \left(x^{2} + 3 x + 4\right)^{2} + 16\)

expr = x**3 + 5*x*y + x*z + 5*y**2 + y*z**2 + z**2
display(expr)
display(expr.subs([(x, 1), [y, 3], [z, 2]]))
display(expr.subs([(x, x), [y, x], [z, x]]))

\(\displaystyle x^{3} + 5 x y + x z + 5 y^{2} + y z^{2} + z^{2}\)

\(\displaystyle 79\)

\(\displaystyle 2 x^{3} + 12 x^{2}\)

無理数の計算

pi.evalf(30) # 1000ぐらい余裕

\(\displaystyle 3.14159265358979323846264338328\)

expr = sin(1)**2 + cos(1)**2
display((expr - 1).evalf())
display((expr - 1).evalf(chop = True))

\(\displaystyle -4.0 \cdot 10^{-124}\)

\(\displaystyle 0\)

numpy関数の利用

import numpy 
a = numpy.arange(100) 
expr = sin(x)**2 + cos(x)**3
display(expr)
f = lambdify(x, expr, "numpy") # numpyを使わないと遅い
f(a)

\(\displaystyle \sin^{2}{\left (x \right )} + \cos^{3}{\left (x \right )}\)

array([ 1.        ,  0.86580202,  0.75475425, -0.95036208,  0.29348079,
        0.94236043,  0.96327991,  0.86012577,  0.97574947, -0.58654075,
       -0.29478182,  0.9999805 ,  0.88881004,  0.92378616,  0.98385952,
       -0.01556116, -0.79539238,  0.90345118,  0.85189192,  0.98895835,
        0.90142733,  0.53566986, -0.99980414,  0.5648118 ,  0.89639378,
        0.99135691,  0.85223411,  0.88972231, -0.81857033,  0.02180437,
        0.97987667,  0.92866036,  0.88463129,  0.99982139, -0.33110279,
       -0.55466852,  0.98152993,  0.86256579,  0.95901972,  0.94785945,
        0.25853532, -0.93733553,  0.7760188 ,  0.86290782,  0.99984336,
        0.86900634,  0.73250099, -0.96191188,  0.32789382,  0.93680443,
        0.96737779,  0.85798037,  0.96910399, -0.61757906, -0.25811662,
        0.99952124,  0.89314651,  0.91891613,  0.98748892, -0.0530218 ,
       -0.77101524,  0.91618973,  0.85189896,  0.98629159,  0.90661841,
        0.50571095, -0.99823778,  0.59310782,  0.89154121,  0.99347747,
        0.85292682,  0.87499481, -0.84050791,  0.05901876,  0.97557656,
        0.93351789,  0.88062925,  0.99901083, -0.36702009, -0.5220171 ,
        0.98646955,  0.8652878 ,  0.95461493,  0.95326866,  0.22310345,
       -0.92285576,  0.79628569,  0.86033431,  0.99937407,  0.87250858,
        0.70927031, -0.97196403,  0.3617301 ,  0.93122345,  0.97129634,
        0.85614111,  0.96157091, -0.64772986, -0.22116672,  0.99847744])

極限/微分/積分

expr = sin(x)/x
# 極限をその場で計算しないならLimitを、するならlimitを使う(以下同様に2種類の関数がある)
display(Eq(Limit(expr, x, 0), limit(expr, x, 0)))
# 極限の方向は片方しか指定できないみたい
display(Eq(Limit(expr, x, 0, dir="-"), limit(expr, x, 0, dir="-"))) 

\(\displaystyle \lim_{x \to 0^+}\left(\frac{1}{x} \sin{\left (x \right )}\right) = 1\)

\(\displaystyle \lim_{x \to 0^-}\left(\frac{1}{x} \sin{\left (x \right )}\right) = 1\)

expr = exp(sin(x))
Eq(Derivative(expr), diff(expr)) 

\(\displaystyle \frac{d}{d x} e^{\sin{\left (x \right )}} = e^{\sin{\left (x \right )}} \cos{\left (x \right )}\)

expr = x * (exp(x)*sin(x) + exp(x)*cos(x))
Eq(Integral(expr), integrate(expr))

\(\displaystyle \int x \left(e^{x} \sin{\left (x \right )} + e^{x} \cos{\left (x \right )}\right)\, dx = x e^{x} \sin{\left (x \right )} - \frac{e^{x}}{2} \sin{\left (x \right )} + \frac{e^{x}}{2} \cos{\left (x \right )}\)

expr = exp(-x**2)
# 無限大は"oo"で指定
Eq(Integral(expr, (x, -oo, oo)), integrate(expr, (x, -oo, oo)))

\(\displaystyle \int_{-\infty}^{\infty} e^{- x^{2}}\, dx = \sqrt{\pi}\)

級数

expr = exp(x)
Eq(expr, expr.series(x=x, x0=0, n=10))

\(\displaystyle e^{x} = 1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + \frac{x^{6}}{720} + \frac{x^{7}}{5040} + \frac{x^{8}}{40320} + \frac{x^{9}}{362880} + \mathcal{O}\left(x^{10}\right)\)

expr = exp(sin(x)**2)
Eq(expr, expr.series(x=x, x0=0, n=10))

\(\displaystyle e^{\sin^{2}{\left (x \right )}} = 1 + x^{2} + \frac{x^{4}}{6} - \frac{11 x^{6}}{90} - \frac{71 x^{8}}{2520} + \mathcal{O}\left(x^{10}\right)\)

プロット

x = symbols('x')
plot(x*x)
plot(sin(x))

png

png

<sympy.plotting.plot.Plot at 0x271b9e4d828>

正規表現の実装(python)

前回の記事で作ったεNFAを使って正規表現を実装できたのでまとめます。ソースは前回と同じくこのリポジトリに。

実装した正規表現

正規表現のすべての機能を実装するのは大変なので、実装するのは基本的なものだけにしました。

実装したのは

  • 一文字 (“a"など)
  • 文字の連結 (“ab"など)
  • 文字の和集合 (“(a|b)"など)
  • 表現の繰り返し(“a*"など)
  • 任意の一文字 (“."のこと)

のみです。教科書には定義に"Φ"(空文字列)がありましたが使い所がなさそうなので代わりに".“を実装しました。

“a+”(aの一回以上の繰り返し)や"(^a|b)“(aでもbでもない文字全て)などの機能は上の表現を組み合わせれば作れます。

他の特徴として入力される文字集合を指定しなければならないというのがあります。これは正規表現の都合というよりもこれまで作ってきたεNFAが文字集合がないと定義できないという都合のためです。

インターフェース

関数regex_to_eNFA()で正規表現をεNFAに変換し、eNFA.run()で文字列が正規表現にマッチするかどうかを調べます。

enfa = regex_to_eNFA("10*1", {"0", "1"})
enfa.run("101")   # -> True
enfa.run("10101") # -> False
enfa.run("10001") # -> True

実装の手順

正規表現をεNFAに変換するプログラムは大きく分けて3つの部分からなります。

  • 正規表現を一つづつに分ける
  • 一つづつに分けた正規表現をεNFAに変換する
  • 変換されたεNFAをε遷移で繋いで一つのεNFAを作る

この部分をコードにすると

def regex_to_eNFA(string, alphabet):
    """ 正規表現stringを認識するeNFAを返す """
    return eNFA.serial_connect([single_regex_to_eNFA(single_regex, alphabet) for single_regex in concat_split(string)])

となります。一行にまとめているので読みにくいですが、concat_split()で正規表現を一つづに分け、single_regex_to_eNFA()で一つづつに分けた正規表現をεNFAに変換し、最後にeNFA.serial_connect()でεNFAを繋いでいます。

順に説明します。

concat_split()

concat_split()は正規表現を一つづつに分ける関数です。一つづつと言うのは連結されていないという意味です。例えば、"a(b|c)de“は["a”, “(a|b)”, “d”, “e”]と分かれます。

ここでの詰まったのは、処理を始める前に一番外側のはずせるカッコを外しておく必要がある点です。もしカッコを外さないと、例えば"(a(b|c)d)“という正規表現は[”(a(b|c)d)“]としか分解されず、無限ループに陥ってしまうからです。対策としてremove_outer_branket()を作り、これに通してから処理を始めるようにしました。

def remove_outer_branket(string):
    """ 一番外側のカッコが外せるなら外して、外せないならそのまま返す """
    if not (string[0] == "(" and string[-1] == ")"):
        return string

    blanket_count = 0
    for char in string[:-1]:
        if char == "(":
            blanket_count += 1
        elif char == ")":
            blanket_count -= 1
        if blanket_count == 0:
            return string
        if blanket_count == 1 and char == "|":
            return string
    return remove_outer_branket(string[1:-1])

def concat_split(string):
    """ 正規表現の連結stringを正規表現一つづつに切って返す """
    ans = []
    right = left = 0
    string = remove_outer_branket(string)
    while left < len(string):
        if string[left] == "(":
            blanket_count = 1
            while blanket_count > 0:
                left += 1
                if string[left] == "(":
                    blanket_count += 1
                elif string[left] == ")":
                    blanket_count -= 1
        if left + 1 < len(string) and string[left + 1] == "*":
            left += 1
        ans.append(string[right:left + 1])
        right, left = left + 1, left + 1
    return ans

テストしてみます。なおgituhubのソースではdoctestを使ってテストをコメントに書いています。改変したときにあってるかどうかがすぐに調べられて便利です。

>>> remove_outer_branket("(abc)")
'abc'
>>> remove_outer_branket("(((abc))")
'abc'
>>> remove_outer_branket("(a|b)(c|d)")
'(a|b)(c|d)'
>>> concat_split("ab")
['a', 'b']
>>> concat_split("a*")
['a*']
>>> concat_split("(abc)*(a|b|c)(cba)*abc")
['(abc)*', '(a|b|c)', '(cba)*', 'a', 'b', 'c']
>>> concat_split("...")
['.', '.', '.']
>>> concat_split("((a|A)(b|B)(c|C))")
['(a|A)', '(b|B)', '(c|C)']
>>> concat_split("((a|A)(b|B)(c|C))*")
['((a|A)(b|B)(c|C))*']

single_regex_to_eNFA()

single_regex_to_eNFA()は連結のない正規表現を受け取り、それを認識するεNFAを返す関数です。

連結のない正規表現

  • 一文字 (“a"など)
  • 文字の和集合 (“(a|b)"など)
  • 表現の繰り返し(“a*"など)
  • 任意の一文字 (“."のこと)

のどれかなので、場合分けして処理します。

def is_dot(string):
    """ stringが任意の一文字を表す記号"."ならTrue """
    return string == "."

def is_atom(string):
    """ stringが通常の文字一つならTrue、それ以外はFalse """
    return True if len(string) == 1 and string not in "()|*." else False

def is_repeat(string):
    """ stringが"X*"ならTrue、それ以外はFalse """
    return len(string) > 1 and string[-1] == "*"

def is_union(string):
    """ stringが文字の集合を表す正規表現ならTrue """
    return string[0] == "(" and string[-1] == ")"

def single_regex_to_eNFA(string, alphabet):
    """ 一項からなる正規表現stringを認識するeNFAを返す """
    if is_dot(string):
        return eNFA.any_word(alphabet)
    if is_atom(string):
        return eNFA.one_word(string, alphabet)
    elif is_repeat(string):
        return eNFA.repeat(regex_to_eNFA(string[:-1], alphabet))
    elif is_union(string):
        return eNFA.parallel_connect([regex_to_eNFA(sub_regex, alphabet) for sub_regex in union_split(string)])
    else:
        assert False, "invalid regular expression: %s" % string

わかりやすいようにis_…()という関数を作っていますが、判定は文脈に依存しています。例えばis_union()は場合分けの最後に置かれなければ正しく動きません。なので関数にしなくても良かったかもしれません。

εNFAを作る関数

any_word(), one_word(), repeat(), parallel_connect(), serial_connact() は全てεNFAを作る関数です。一例としてserial_connect()を示し、解説は省略します。

def serial_connect(automaton):
    """ 複数のオートマトンを横並びに一つにまとめる """
    alphabet = set()
    states = set()
    transitions = {}
    for i, automata in enumerate(automaton):
        # alphabetを更新
        alphabet = alphabet.union(automata.alphabet)

        # statesを更新
        # 状態とautomataのインデックスを組にすることで唯一性を確保
        states = states.union(set([(i, state) for state in automata.states]))

        # automata内部の遷移をtransitionsに追加
        for final_state in automata.final_states:
            transitions[(i, final_state)] = {}
        for state, trans_dict in automata.transitions.items():
            transitions[(i, state)] = {}
            for char, states_transit_to in trans_dict.items():
                transitions[(i, state)][char] = frozenset([(i, state) for state in states_transit_to])

        # 一つ前のautomataからのイプシロン遷移を追加
        if i == 0:
            continue
        for pre_final_state in automaton[i - 1].final_states:
            transitions[(i - 1, pre_final_state)][-1] = frozenset([(i, automata.init_state)])

    # 始状態と終状態を作る
    init_state = (0, automaton[0].init_state)
    final_states = set([(len(automaton) - 1, f) for f in automaton[-1].final_states])
    return NFAWithEpsilonTransition(states, alphabet, transitions, init_state, final_states)

DFAへの変換と最小化

εNFAのDFAへの変換とDFAの最小化はこれまでに作ったので、正規表現から作ったεNFAを最小DFAに変換することができます。

これによって正規表現の性質が分かります。

>>> from Regex import regex_to_eNFA
>>> any = regex_to_eNFA("(0|1)*", {"0", "1"})
>>> any_mindfa = any.convert_to_DFA().minimized()
Automata
    states      : {{{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}}
    alphabet   : {'0', '1'}
    transitions : {{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}} : {'0': {{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}, '1': {{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}}
    init_state  : {{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}
    final_states: {{{(0, (0, (1, (0, 0)))), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (0, (0, 1))))}, {(0, 1), (0, -1), (0, (0, 1)), (0, (0, (0, (0, 0)))), (0, (0, (1, (0, 0))))}, {(0, -1), (0, (0, 1)), (0, (0, (1, (0, 1)))), (0, (0, (1, (0, 0)))), (0, (0, (0, (0, 0))))}}}

とても見にくいですが、よく見ればこのDFAは状態が一つしか無い、どんな文字列も認識するDFAだとわかります。正規表現"(0|1)“が表すのは任意の文字列なので、このDFAが正しく変換されていることがわかります。

最小DFAは全て同型なので、DFAが同型かどうかを判定する関数を作れば正規表現が同値かどうかを判定できるようになるはずです(やってないけど)。

またDFAの表記がみにくい問題については、DFAのコンストラクタで状態の名前を振り直す処理を加えればマシになるはずです。デバッグがしにくくなるので実装はしませんでしたが。

pythonで無限リスト

 

 

pythonで無限リスト

この前pythonのジェネレータを使って作った無限リストには、一度計算した値を覚えておく機能がないため実行速度がとても遅いという問題がありました。SICPを参考にして試行錯誤した結果ラムダ式を使った遅延評価を使うことでメモ化付きの無限リストを作ることができたので書きます。

参考文献

参考にしたのは前の記事と同じでstructure and interpretation of computer programsの3.5節です。メモ化関数(memo_proc)を使った高速化が参考になりました。また後半の級数についての関数は練習問題3.53-3.61を参考にしました。

遅延評価

遅延評価とは式の値を定義されたときには計算せず、必要になってから初めて計算するという計算方法です。つまりa = f(x)と変数aを定義してもそれだけではf(x)が計算されるとは限らず、別の部分(たとえば"b = a + 1")でaが参照されて初めてf(x)が計算されます。aが参照されなければ、f(x)が計算されないでプログラムが終わることもありえます。

pythonの評価順序

pythonのプログラムを全く変更しないで遅延評価を使うことはできません。なぜなら、pythonが「変数に式の値を代入するときは式を計算する」と言う決まりを持っているためです。

def f(x):
    """ 呼ばれるとf calledとprintする関数 """
    print("f called.")
    return x + 1

a = f(0) # こう書いた時点でf(x)は計算される
print("a is", a)
b = a + 10
print("b is", b)
f called.
a is 1
b is 11

"f called"が"a is 1"よりも先に出力されていることから、f(0)が実行されているのは"b = a + 10"が実行されているときではなくa = f(0)が実行されているときだとわかります。

遅延評価の実装

そんなpythonで遅延評価を行うために、遅延評価させたい式を呼ぶゼロ引数関数を変数に代入し、値が必要になったときはその関数を呼ぶという工夫をします。pythonには無名関数を作るラムダ記法があるので、以下のように書けます。

a = lambda: f(0)     # aを「f(0)を呼び出す引数を取らない関数」とする
print("a is", a)
b = a() + 10         # aを呼び出す、すなわちf(0)の値を取得する
print("b is", b)
print("a() is", a()) # もう一度aを呼び出す
a is <function <lambda> at 0x000001FBDE38FBF8>
f called.
b is 11
f called.
a() is 1

出力の一行目から、aがf(0)の値ではなくて関数であることがわかります。2, 4行目から、aの値が参照されるとそのたびにf(0)が計算されていることがわかります。

遅延評価の組み合わせ

遅延評価される値を組み合わせて別の値を遅延評価されるようにすることもできます。

a = lambda: f(0)
b = lambda: a() + 10
c = lambda: a() + b() + 10

print("definition end.")
print("c() is", c())
definition end.
f called.
f called.
c() is 22

b()を計算するときにはa()が一回呼ばれ、c()を計算するときにはa()とb()が一回ずつ呼ばれるので、c()を計算するときにf(0)が二回呼ばれます。このことは"f called "が二回呼ばれていることから確かめられます。またa,b,cを定義しただけではf(0)は呼ばれないので、出力の最初の行は"definition end."となっています。

無限リスト

この遅延評価の仕組みを使うと、値を無限に持つリストが作れます。

簡単な無限リスト

リストが値を計算するのは先頭の要素だけで、2番目以降の要素は遅延評価する、という考え方が基本となります。つまり、

some_list = (first_value, lambda: list_tail()) # list_tail()はリストを返す

と言う形を用い、list_tail()の中でsome_listを再帰的に用いればsome_listの長さが無限になるという仕組みです。一番簡単な例を示します。

ones = (1, lambda: ones)

for _ in range(5):
    print(ones)
    print(ones[0])
    ones = ones[1]() # リストを一つ進める
(1, <function <lambda> at 0x000001FBDE803E18>)
1
(1, <function <lambda> at 0x000001FBDE803E18>)
1
(1, <function <lambda> at 0x000001FBDE803E18>)
1
(1, <function <lambda> at 0x000001FBDE803E18>)
1
(1, <function <lambda> at 0x000001FBDE803E18>)
1

「onesの第二項以降はonesそのものである」と定義しているので、for文の中の"ones = ones[1]()"でonesを一つ進めているにもかかわらず、どれだけ項を進めてもones[0]で取り出されるのは1です。

ユーティリティ関数

便利のためにユーティリティ関数を定義しておきます。

def car(seq):
    """ seqの初項を返す """
    return seq[0]

def cdr(seq):
    """ seqの第2項以降を返す """
    return seq[1]()

def seq_ref(seq, n=0):
    """ seqの第n項を表示する """
    for _ in range(n):
        seq = cdr(seq)
    print(car(seq))

def seq_print(seq, n=10, end=", "):
    """ seqの最初のn項を表示する """
    for i in range(n):
        print(car(seq), end=end)
        seq = cdr(seq)
    print("")
        
def add(s, t):
    """ 数列としての和s+tを返す """
    return (car(s) + car(t), lambda: add(cdr(s), cdr(t)))

無限リストの例

# 自然数
integer = (1, lambda: add(integer, ones))
seq_print(integer)
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
# フィボナッチ数列
fib = (1, lambda: (1, lambda: add(cdr(fib), fib)))
seq_print(fib)
1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 

メモ化

問題点

ここまで使ってきた遅延評価の仕組みには同じ計算を不必要に何度もしてしてしまうという問題点があります。たとえば、この例ではf(0)が2回計算されます。

a = lambda: f(0)
print(a())
print(a())
f called.
1
f called.
1

ここではf(0)の計算にかかる時間が僅かなのであまり気になりませんが、再帰的な定義などによって何回も呼ばれるときは計算の負荷が無視できなくなります。前に定義したfib()で第30項を計算するのは厳しいでしょう。

メモ化の実装

そこで一回目の呼び出しのときに値を保存しておいて、2回目以降の呼び出しでは保存した値をそのまま返すようにします。

そのために、今まで考えてきた形の遅延評価のための関数を受取り、値を保存するようにした関数を返す関数memoize()をつくります。pythonではメソッドが変数を持つことができるので、memoize()の中で新しく関数を定義してその関数自身に関数の値をresultとしてもたせます。

def memoize(func):
    """ funcに値を保存する機能を追加したものを返す """
    def temp():
        if temp.result is None:  # もしresultが計算されていなかったら...
            temp.result = func() # funcの値を計算し、保存する
        return temp.result       # 保存された値を返す
    temp.result = None
    return temp # 新しく作った関数を返す

memoize()を使うと、遅延評価を何度行っても関連づいた関数が呼ばれるのは一度だけです。

a = memoize(lambda: f(0))
print(a())
print(a())
print(a())
print(a())
f called.
1
1
1
1

メモ化付きの無限リスト

無限リストにもメモ化を組み込むために、関数を新たに一つ定義しadd()を書き直します。

def cons(a, b):
    """ 初項a, 第二項以降がbであるリストを返す"""
    return (a, memoize(b))


def add(s, t):
    """ 数列としての和s+tを返す """
    return cons(car(s) + car(t), lambda: add(cdr(s), cdr(t)))

このconsとaddを使えばメモ化付きの無限リストが作れます。

メモ化によってfibも常識的な速度で計算できます。

fib = cons(1, lambda: cons(1, lambda: add(cdr(fib), fib)))
seq_ref(fib, 1000)
70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501

無限リストの例

もういくつか無限リストの例を示します。

# 2の階乗
two_powered = cons(1, lambda: add(two_powered, two_powered))
seq_print(two_powered)
1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 
def mul(s, t):
    """ 数列としての積s*tを返す """
    return cons(car(s) * car(t), lambda: mul(cdr(s), cdr(t)))

# 階乗
fact = cons(1, lambda: mul(fact, integer))
seq_print(fact)
1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 
def sum(s):
    """ sの累積和 """
    ans = cons(0, lambda: add(s, ans))
    return ans

# 自然数の累積和
seq_print(sum(integer))
0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 
def expand(num, den, radix):
    """ x_n = num/denのradix進表示の小数点n桁目 """
    ans = cons((num * radix) // den, lambda: expand((num*radix)%den, den, radix))
    return ans

print("1/7の10進表示")
seq_print(expand(1, 7, 10))
print("1/7の小数点第1000桁")
seq_ref(expand(1, 7, 10), 1000)
print("1/3の2進表示")
seq_print(expand(1, 3, 2))
1/7の10進表示
1, 4, 2, 8, 5, 7, 1, 4, 2, 8, 
1/7の小数点第1000桁
5
1/3の2進表示
0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 

幾何級数を計算する

無限といえば級数です。

\(e\) の計算

手始めに\(e\)を計算します。

\[e = \frac{1}{1!} + \frac{1}{2!} + \frac{1}{3!} + \frac{1}{4!} + ...\]

なので、\(e\)に収束する数列approximate_e()を定義できます。

def map(s, f):
    """ sの各項にfを適用する """
    return cons(f(car(s)), lambda: map(cdr(s), f))

approximate_e = sum(map(fact, lambda x: 1/x))
print("eの近似値")
seq_print(approximate_e, end="\n")
print("第300項")
seq_ref(approximate_e, 300)
eの近似値
0
1.0
2.0
2.5
2.6666666666666665
2.708333333333333
2.7166666666666663
2.7180555555555554
2.7182539682539684
2.71827876984127

第300項
2.7182818284590455

\(e^x\) の計算

さらに一般化して、\(x\)を受け取り\(e^x\)に収束する数列を返す関数eを定義できます。

def scale(s, k):
    """ sをk倍したものを返す """
    return map(s, lambda x: k * x)

def x_powers(x):
    """ 初項init, 公比xの等比級数を返す """
    return cons(1, lambda: scale(x_powers(x), x))

def series(s):
    """ sを係数に持つ幾何級数の各項を計算する関数を返す """
    return lambda x: mul(s, x_powers(x))

def series_sum(s):
    """ sを係数に持つ幾何級数の累積和を計算する関数を返す """
    return lambda x: sum(series(s)(x))

e_coeff = map(fact, lambda x: 1/x) # e^xの係数
e = series_sum(e_coeff)
seq_ref(e(1), 30)
seq_ref(e(2), 30)
seq_ref(e(3), 30)
seq_ref(e(4), 30)
2.7182818284590455
7.389056098930649
20.08553692318766
54.59815003314426

直接的な \(e^x\) の計算

上のようにすれば\(e^x\)の計算ができるのですが、この定義の仕方だとはっきり言って普通に定義したほうが簡単にできます。つまり"e_coeff = map(fact, lambda x: 1/x)"などと定義しなくても

def fact_func(n):
        return n * fact_func(n-1) if n > 1 else 1

def e(x, n):
    ans = 0
    for i in range(n):
        ans += x ** i / fact_func(i)
    return ans

print(e(1, 30))
print(e(2, 30))
print(e(3, 30))
print(e(4, 30))
2.7182818284590455
7.389056098930649
20.08553692318766
54.59815003314426

とできるのです。こうした方が読みやすい上に遅延評価がどうだとか考える必要もありません。このままだと遅延評価の意味が無いので、\(e^x\)をもう少し面白い方法で定義し直します。

\(e^x\) の陰な定義

注目するのは、「\(e^x\)積分すると定数部分を除いて\(e^x\)に等しい」という数学的な事実です。幾何級数に対して積分を行うと項が一つ後ろにずれるので、もと級数積分した級数の間に関係があれば、級数を陰に定義することができます。

まず、幾何級数の係数を受け取ってそれを積分した級数の係数を返す関数を書きます。

def integrate_series(s):
    """ 級数sを積分した級数を返す、ただし定数項は(不定なので)除く """
    return mul(s, map(integer, lambda x: 1/ x))

seq_print(fact)
e_coeff = map(fact, lambda x: 1/x) # e^xの係数
print("e^xの係数")
seq_print(e_coeff, 7)
print("e^xを積分した係数(初項を除く)")
seq_print(integrate_series(e_coeff), 7)
1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 
e^xの係数
1.0, 1.0, 0.5, 0.16666666666666666, 0.041666666666666664, 0.008333333333333333, 0.001388888888888889, 
e^xを積分した係数(初項を除く)
1.0, 0.5, 0.16666666666666666, 0.041666666666666664, 0.008333333333333333, 0.0013888888888888887, 0.0001984126984126984, 

\(e^x\)積分したときの係数が\(e^x\)の係数と初項を除いて同じことがわかります。

このことを用いて、e_coeffを再帰的に(陰に)定義します。

# e_coeffは初項が1で、それ以降は「e_coeffを積分した級数から初項を除いたもの」である
e_coeff = cons(1, lambda: integrate_series(e_coeff))
e = series_sum(e_coeff)
seq_ref(e(1), 30)
seq_ref(e(2), 30)
seq_ref(e(3), 30)
seq_ref(e(4), 30)
2.7182818284590455
7.389056098930649
20.08553692318766
54.59815003314426

陰に定義できる他の関数

級数の陰な定義の便利な点は\(e^x\)以外の関数を定義するときにも使える点です。たとえばsinとcosは次のように定義できます。

from math import pi

cos_coeff = cons(1, lambda: scale(integrate_series(sin_coeff), -1))
sin_coeff = cons(0, lambda: integrate_series(cos_coeff))
cos = series_sum(cos_coeff)
sin = series_sum(sin_coeff)

print("cos")
seq_ref(cos(pi / 6), 30)
seq_ref(cos(pi / 4), 30)
seq_ref(cos(pi / 3), 30)
print("sin")
seq_ref(sin(pi / 6), 30)
seq_ref(sin(pi / 4), 30)
seq_ref(sin(pi / 3), 30)
cos
0.8660254037844386
0.7071067811865475
0.5000000000000001
sin
0.49999999999999994
0.7071067811865475
0.8660254037844386

級数の積

ところで、\({sin}^2(x) + {cos}^2(x) = 1\)\(x\)がいくらであってもいつでも成り立つはずです。このことを確認してみます。

そのためには\({sin}^2, {cos}^2\)級数展開を求める必要があるので、級数の掛け算をする関数mul_coeff()を定義します。mul_coeff()の計算には

\[(a_0+a_1x+a_2x^2+\cdots)(b_0+b_1x+b_2x^2+\cdots) = a_0b_0 + a_0(b_1x+b_2x^2+\cdots) + (a_1x+a_2x^2+\cdots)(b_0+b_1x+b_2x^2+\cdots)\]

という関係を利用します。

def mul_coeff(s, t):
    """ sとtの多項式としての積 """
    return cons(car(s) * car(t), lambda: add(scale(cdr(t), car(s)), mul_coeff(cdr(s), t)))

maybe_one_coeff = add(mul_coeff(sin_coeff, sin_coeff), mul_coeff(cos_coeff, cos_coeff)) 
maybe_one = series_sum(maybe_one_coeff)

seq_ref(maybe_one(0.1), 10)
seq_ref(maybe_one(0.5), 10)
seq_ref(maybe_one(1.0), 10)
seq_ref(maybe_one(2.0), 10)
seq_ref(maybe_one(3.0), 10)
seq_ref(maybe_one(4.0), 10)
seq_ref(maybe_one(10.0), 10)
1.0
1.0
1.0
1.0000000000000004
0.9999999999999993
0.9999999999999858
0.9999999999936162

\({sin}^2(x) + {cos}^2(x)\)がほぼ1になっていることが確認できます。

級数の割り算

掛け算ができたので割り算もやります。割り算は掛け算より複雑です。

級数の逆数

まず多項式の逆数を考えます。\(S\)級数として

\[ 1 \div S = X \]

を求めることを考えます。

\(S=s_c + S_r\)(\(s_c\)は定数、\(S_r\)級数)として

\[ S \cdot X = 1\] \[ (s_c + S_r) \cdot X = 1\] \[ s_c \cdot X + S_r \cdot X = 1\] \[ X = \frac{1}{s_c}(1 - S_r \cdot X) \]

となります。\(S_r\)には定数項が含まれないので、最後の式の右辺にある\(S_r \cdot X\)にも定数項が含まれないことがポイントです。これによって最初にXの定数項を計算し、それ以外の項は\(-S_r \cdot X\)の部分から逐次に計算することができます。

コードにします。

def reciprocal(s):
    """ 級数としての逆数1/sを返す """
    return cons(1/car(s), lambda: scale(mul_coeff(cdr(s), reciprocal(s)), -1))

あっているか確認するために\(sec(x)\)を計算してみます。

sec_coeff = reciprocal(cos_coeff)
sec = series_sum(sec_coeff)

maybe_one_coeff = mul_coeff(cos_coeff, sec_coeff) 
maybe_one = series_sum(maybe_one_coeff)

print("cosの係数")
seq_print(cos_coeff, 10)
print("secの係数")
seq_print(sec_coeff, 10)
print("maybe_oneの係数")
seq_print(maybe_one_coeff, 10)
cosの係数
1, -0.0, -0.5, 0.0, 0.041666666666666664, -0.0, -0.0013888888888888887, 0.0, 2.4801587301587298e-05, -0.0, 
secの係数
1.0, 0.0, 0.5, -0.0, 0.20833333333333334, -0.0, 0.08472222222222223, -0.0, 0.034350198412698416, -0.0, 
maybe_oneの係数
1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 

maybe_oneの計算はせずに係数を表示させました。cosとsecは係数に値を持っているにも関わらずそれらの積は係数が0なことから、たしかにsecがcosの逆数なことがわかります。

一般の級数の割り算

逆数の計算ができれば一般の割り算も簡単です。

def div_coeff(s, t):
    """ 級数としての逆数s/tを返す """
    return mul_coeff(s, reciprocal(t))

\(tan(x) = sin(x) \div cos(x)\)で確認します。

tan_coeff = div_coeff(sin_coeff, cos_coeff)
tan = series_sum(tan_coeff)

seq_ref(tan(0.0), 30)
seq_ref(tan(pi/6), 30)
seq_ref(tan(pi/4), 30)
seq_ref(tan(pi/3), 30)
0.0
0.5773502691896235
0.9999999992094688
1.7320428394905787

感想

SICPではこのあとにエイトケンのΔ2乗加速法を使って円周率を求める話があるのですがこの記事には書きませんでした。

疲れた。

ジェネレータで遊ぶ(python)

 

pythonにはジェネレータという組み込み型が存在し、値を逐次に計算する配列として使うことができます。ジェネレータを使って遊んだので書きます。

参考文献

structure and interpretation of computer programsの3.5節を参考にしました。

無限ジェネレーターに関するユーティリティ関数

この文章では無限に値を返すジェネレーターを扱うので、最初に便利な関数を定義しておきます。

# gの最初のn項を表示する
def g_print(g, n=10):
    for i in range(n):
        print(next(g))

# gをn項飛ばす
def g_delayed(g, n):
    for _ in range(n):
        next(g)
    return g

# gの第n項を表示する
def g_ref(g, n):
    return next(g_delayed(g, n))

無限ジェネレータ

pythonではジェネレータを使って無限に値を返し続けるものを書けます。

# 自然数を0から返す
def integers():
    i = 0
    while True:
        yield i
        i += 1
>>> g_print(integers())
0
1
2
3
4
5
6
7
8
9

同様に少し凝ったものも書けます。

# 5の倍数でない自然数を小さい順に返す
def no_five():
    i = 0
    while True:
        if i % 5 != 0:
            yield i
        i += 1
>>> g_print(no_five())
1
2
3
4
6
7
8
9
11
12

このno_fiveではローカル変数iを1づつ増やすことで全ての自然数をチェックしています。

no_fiveで行っている「ローカル変数に無限に自然数を代入して何かする」ということはまさにintegersがfor文で使われたときに行うことなので、次のように書くことができます。

def no_five():
    for i in integers() # while文の代わりにfor文とintegersを使った
        if i % 5 != 0:
            yield i

integersからno_fiveが作られる関係は、高階関数を使ってまとめることができます。

# ジェネレータg(無限でも良い)からcondを満たすものだけを抜き取る
def g_filter(g, cond):
    for i in g:
        if cond(i): # "i % 5 != 0" をcondに一般化した
            yield i

no_five = g_filter(integers(), lambda x: x % 5 != 0)

なおpythonの組み込み関数にfilterがありますが、無限ジェネレータに対しては上手く動きません。要素をすべて読みに行ってしまうようです。

g_filterと同様に、ジェネレータに対する高階関数を定義できます。

# gの全要素にfを適用する
def g_map(g, f):
    for i in g:
        yield f(i)

# gの要素とhの要素を足す
def g_add(g, h):
    for a, b in zip(g, h):
        yield a + b

# gの要素とhの要素をかける
def g_mul(g, h):
    for a, b in zip(g, h):
        yield a * b

# gの要素の累積和
def g_sum(g):
    sum = 0
    for i in g:
        sum += i
        yield sum

これらを使って様々なジェネレータを定義できます。

# 偶数全部
evens = g_map(integers(), lambda x: 2*x)

# 7で割ると3余る偶数全部
Nseven_plus_three = g_filter(evens, lambda x: x % 7 == 3)

# 7で割ると3余る偶数に自然数を足したもの
Nseven_plus_three_plus_N = g_add(Nseven_plus_three, integers())

再帰ジェネレータ

ジェネレータは関数の一種なので、再帰的に定義することができます。

def ones():
    yield 1
    yield from ones() # 次の要素は自分自身から取る

定義が言っているのは「onesは最初1を返し、それ以降はonesが返すものを最初から返す」ということです。つまり、1を無限に返します。

>>> g_print(ones())
1
1
1
1
1
1
1
1
1
1

ただしこの例では、onesはyieldするたびに新しくonesを作成するので、onesはメモリーを無限に食います。そのためこのonesは実用的ではありません。

もう少し面白い例がこれです。

def fib():
    yield 1
    yield 1
    yield from g_add(g_delayed(fib(), 1), fib())

fibの定義を解釈すると、

  1. fibの最初は1
  2. その次は1
  3. その次以降は「fibと、fibを一つずらしたジェネレータの和」というジェネレータの返す値を返す

3.が言っているのは\(x_n = x_{n-1} + x_{n-2}\)ということなので、fibはフィボナッチ数列を表します。

>>> g_print(fib())
1
1
2
3
5
8
13
21
34
55

メモ化していないので低速ですが、ifを使わないフィボナッチ数列の定義というのは珍しいと思います。

再帰的に定義されたジェネレータの例をいくつかあげます。他では見ないような形で再帰が使われているので面白いです。

# 1, 2, 3, 4, 5, ...
def integers():
    yield 1
    yield from g_add(ones(), integers())

# 1, 2, 4, 8, 16, ...
def two_powered():
    yield 1
    yield from g_add(s(), s())

# 1, 2, 6, 24, 120, 720, ...
def fact():
    yield 1
    yield from g_mul(g_delayed(integers(), 2), fact())

感想

本当はこのあと再帰ジェネレータをメモ化して速くして終わる予定だったんだけどメモ化ができなかったのでここで終わる。こういう細かいところはpythonじゃ無理だ。

オートマトンをpythonで書く(続き)

この前作ったDFAをいくらか拡張したのでまとめます。ソースはgithubに。

NFA

講義ではDFAの次に非決定性有限オートマトン(Nondeterministic Finite Automata, NFA)が出てきたので、最初の拡張としてNFAを作りました。

NFAは文字列を読み込んでいる途中に状態を複数取ることができます。状態を複数取っているときに文字を読んだ場合、次の状態は現在の各状態の遷移先の和集合です。また文字列を認識するかどうかは最後まで文字を読み込んだあとの状態に終状態が一つでもあるかどうかで判断します。

こう書くと難しく感じますが、DFAのプログラムと骨格はほとんど同じなので特に迷わずに書くことができました。

class NondeterministicFiniteAutomata(Automata):
    """ 非決定性有限オートマトン """

    def run(self, string):
        """ stringを認識するかチェックする """
        states = {self.init_state}
        for char in string:
            if char not in self.alphabets:
                print("ERROR : invalid character \"%s\"" % char)
                return False
            else:
                states = self.next_states(states, char)
        return len(states.intersection(self.final_states)) != 0

    def next_states(self, states, char):
        """ statesからcharを読んだ時の次の状態集合を返す """
        next_states = set()
        for state in states:
            next_states = next_states.union(self.transitions[state].get(char, frozenset()))
        return frozenset(next_states)

参考までに、前の記事で書いたDFAのrun関数も載せておきます。NFAのrun関数と同じような形をしているのが分かるはずです。

class DeterministicFiniteAutomata(Automata):
    """ 決定性有限オートマトン """

    def run(self, string):
        """ stringを認識するかチェックする """
        state = self.init_state
        for char in string:
            if char not in self.alphabets:
                print("ERROR : invalid character \"%s\"" % char)
                return False
            else:
                state = self.transitions[state][char]
        return state in self.final_states

具体的なDFAは次のように作ります。

from Automata import NondeterministicFiniteAutomata as NFA
# 010を部分列に持つ文字列を認識するNFA
a, b, c, d = range(4)
states = {a, b, c, d}
alphabets = {"0", "1"}
transitions = {
    a: {"0": {a, b}, "1": {a}},
    b: {"1": {c}},
    c: {"0": {d}},
    d: {"0": {d}, "1": {d}},
}
init_state = a
final_states = {d}
has_010 = NFA(states, alphabets, transitions, init_state, final_states)

作った上でrun()すれば動作を確認できます。

# pythonコンソール
>>> has_010.run("0001011")
True
>>> has_010.run("0011")
False
>>> has_010.run("001100")
False
>>> has_010.run("100")
False
>>> has_010.run("0100")
True

NFAからDFAへの変換

任意のNFAを等価なDFAに変換するアルゴリズムがあるのでそれも作りました。状態の冪集合を考えることで、NFAが文字を読む間に取りうる状態の集合をDFAの一つの状態としてしまうことがこのアルゴリズムのポイントです。

下のコードでは目標となるDFAのstates, alphabets, transitions, init_state, final_statesを作ることでDFAを作っています。これらの中でalphabetsとinit_stateは自明であり、final_statesはstatesがわかってしまえばすぐに分かるので、作るときに苦労するのはstatesとtransitionsです。

    def convert_to_DFA(self):
        """ 等価なDFAに変換する """
        # 状態と遷移関数を作る
        states = {frozenset({self.init_state})}
        transitions = {}
        to_search = {frozenset({self.init_state})}
        searched = set()
        while len(to_search) != 0:
            searching = to_search.pop()
            searched.add(searching)
            for char in self.alphabets:
                # 実際に遷移させて考える
                reachables = self.next_states(searching, char)
                # 遷移先の状態集合そのものが一つの状態となる
                states.add(reachables)
                if reachables not in searched:
                    # 遷移させた先がまだ調べていない状態だったらさらに調べる必要がある
                    to_search.add(reachables)
                if not transitions.get(searching):
                    # 二重辞書の要素を一気に作ることはできないので注意
                    transitions[searching] = {}
                # 遷移を追加する
                transitions[searching][char] = frozenset(reachables)

        # DFAの終状態は終状態を一つでも含む状態全体からなる集合
        final_states = set([x for x in states if len(self.final_states.intersection(x)) != 0])
        # DFAの始状態はNFAの始状態だけからなる集合
        init_state = frozenset({self.init_state})
        # DFAのアルファベットはNFAと同じ
        alphabets = self.alphabets
        return DeterministicFiniteAutomata(states, alphabets, transitions, init_state, final_states)

次のように使います。同じ機能を持つDFAができていますが、変換によって状態数が増え、複雑になっていることがわかります。

# has_010を作った状態のpythonコンソール
>>> print(has_010)

Automata
    states      : {0, 1, 2, 3}
    alphabets   : {'1', '0'}
    transitions : 0 : {'1': {0}, '0': {0, 1}}
                  1 : {'1': {2}}
                  2 : {'0': {3}}
                  3 : {'1': {3}, '0': {3}}
    init_state  : 0
    final_states: {3}
>>> has_010_DFA = has_010.convert_to_DFA()
>>> print(has_010_DFA)
Automata
    states      : {{0, 2, 3}, {0, 2}, {0, 1, 3}, {0}, {0, 3}, {0, 1}}
    alphabets   : {'1', '0'}
    transitions : {0, 2, 3} : {'1': {0, 3}, '0': {0, 1, 3}}
                  {0, 2} : {'1': {0}, '0': {0, 1, 3}}
                  {0} : {'1': {0}, '0': {0, 1}}
                  {0, 1, 3} : {'1': {0, 2, 3}, '0': {0, 1, 3}}
                  {0, 3} : {'1': {0, 3}, '0': {0, 1, 3}}
                  {0, 1} : {'1': {0, 2}, '0': {0, 1}}
    init_state  : {0}
    final_states: {{0, 2, 3}, {0, 3}, {0, 1, 3}}

>>> has_010_DFA.run("010")
True
>>> has_010_DFA.run("0110")
False
>>> has_010_DFA.run("100110")
False
>>> has_010_DFA.run("1001010")
True

全てのNFAはDFAに変換できることと全てのDFAはNFAなことから、DFAとNFAが同じ計算能力を持つことがわかります。

εNFA

NFAをさらに拡張してイプシロン動作付き非決定性有限オートマトン(εNFA)も作りました。

εNFAとNFAの違いは、イプシロン遷移があるかどうかです。イプシロン遷移が状態\(p\)から状態\(q\)に張られていた場合、状態\(p\)に到達したら無条件で状態\(q\)にも到達したとみなします。さらに状態\(q\)から別の状態にイプシロン遷移が張られていれば、その状態についても同様です。

違いがこれだけなのでNFAをεNFAにするのは簡単かと思いましたがそれなりに手こずりました。イプシロン遷移を一回たどるだけでは不十分で、飽和するまでたどりきらないといけないためです。

class NFAWithEpsilonTransition(Automata):
    """ ε動作付き非決定性有限オートマトン """

    def reachables_with_a_epsilon_from(self, state):
        """ stateから一回のε遷移だけで到達可能な状態の集合を返す """
        return frozenset(self.transitions[state].get(-1, {}))

    def reachables_with_epsilons_from(self, state):
        """ stateからε遷移だけで到達可能な状態の集合を返す """
        ans = frozenset([state])
        reachables = self.reachables_with_a_epsilon_from(state)
        while not reachables.issubset(ans): # 飽和するまで繰り返す
            ans = ans.union(reachables)
            new_reachables = frozenset()
            for reachable in reachables:
                new_reachables = new_reachables.union(self.reachables_with_a_epsilon_from(reachable))
            reachables = new_reachables
        return ans

    def next_states(self, states, char):
        """ statesからcharを読んだ時の次の状態集合を返す """
        next_states = frozenset()
        for state in states:
            for reachable in self.transitions[state].get(char, frozenset()):
                next_states = next_states.union(self.reachables_with_epsilons_from(reachable))
        return frozenset(next_states)

    def run(self, string):
        """ stringを認識するかチェックする """
        states = self.reachables_with_epsilons_from(self.init_state)
        for char in string:
            if char not in self.alphabets:
                print("ERROR : invalid character \"%s\"" % char)
                return False
            else:
                states = self.next_states(states, char)
        return len(states.intersection(self.final_states)) != 0

以下のようにして動かします。

>>> from Automata import NFAWithEpsilonTransition as eNFA
>>> a, b, c, d, e = range(5)
>>> states = {a, b, c, d, e}
>>> alphabets = {"0", "1", "2", "3", "4"}
>>> transitions = {
...     a: {"0": {a}, -1: {b}},
...     b: {"1": {b}, -1: {c}},
...     c: {"2": {c}, -1: {d}},
...     d: {"3": {d}, -1: {e}},
...     e: {"4": {e}}
>>> }
>>> init_state = a
>>> final_states = {e}
>>> is_increasing = eNFA(states, alphabets, transitions, init_state, final_states)

>>> is_increasing.run("01144")
True
>>> is_increasing.run("01231")
False
>>> is_increasing.run("00004")
True
>>> is_increasing.run("33444")
True

NFAと同様の方法でεNFAもDFAに変換することができます(コードはほとんど同じなので省略します)。

DFAの最小化

任意のDFAは最小化することができます。つまり、任意のDFAに対して、そのDFAが同じ言語を認識するDFAのうち状態数が最小のDFAを作ることができます。

重要となるのは以下の区別できないという関係です。

DFA \(M(Q,\varSigma,\delta,q_0,F)\) について、状態 \(p,q \in Q\)区別できない

\(\Longleftrightarrow\) 全ての文字列\(s \in \varSigma^\ast\)について\(\delta^\ast (p,s) \in F \Leftrightarrow \delta^\ast (p,s) \in F\)

DFAの2つの状態が区別できないとことを言い換えると、2つの状態がDFAの中で同じ役割をしていると言う事ができます。DFAがどちらの状態にあってもそれから読む文字列に対して全く同じ反応を示すからです。「区別できない」という名前もここからきたのでしょう。

あるいは、区別できない2つの状態はDFAの中で無駄な状態だともいえます。なぜなら、区別できない状態の遷移を上手くいじれば2つの状態を1つにまとめてしまえることを示せる定理があるからです。いわゆるMyhill-Nerodeの定理(の一部)です。

DFA \(M(Q,\varSigma,\delta,q_0,F)\)が最小DFAである。

\(\Longleftrightarrow\) 任意の異なる状態の組\((p,q) \in Q \times Q (p \neq q)\)が区別できる。

すごくざっくり言うと、「最小DFAには"無駄な"状態が存在しない / ”無駄な”状態が無いDFAは最小」ということになります。

この定理の証明は区別できないという関係の同値類を使って遷移の具体的な張り方を示すことで行います。ここでは省略します。

というわけでDFAを最小化するには、状態のうち区別できるペアを見つけ、そのペアを一つの状態として良い感じに遷移を張ればよいということになります。以下の性質を使います。どちらも定義から明らかです。

  • \(p \in F\)\(q \notin F\)について、\(p\)\(q\)は区別できる。
  • \((p, q) \in Q \times Q\)について、ある\(c \in \varSigma\)について\(\delta(p, c)\)\(\delta(q, c)\)が区別できるならば、\(p\)\(q\)も区別できる。

これをプログラムにすると以下のようになります。

class DeterministicFiniteAutomata(Automata):
    """ 決定性有限オートマトン """
    ...
    (略)
    ...
    def minimized(self):
        """ 最小化されたDFAを返す """

        # markedとunmarkedを初期化
        # markedは区別できるペアの集合
        # unmarkedは区別できないペアの集合
        marked = set()
        unmarked = set()
        checked = set()
        for p in self.states:
            for q in self.states:
                if frozenset({p, q}) in checked or p == q:
                    continue
                if (p in self.final_states and q not in self.final_states) or \
                        (q in self.final_states and p not in self.final_states):
                    marked.add(frozenset({p, q}))
                else:
                    unmarked.add(frozenset({p, q}))
                checked.add(frozenset({p, q}))
                checked.add(frozenset({q, p}))

        # markedを限界まで増やす
        flag = True
        while flag:
            flag = False
            for p, q in unmarked:
                for s in self.alphabets:
                    # (p,q)をsで遷移させてmarkedにはいるなら(p,q)もmarkedに入る
                    if frozenset({self.transitions[p][s], self.transitions[q][s]}) in marked:
                        flag = True
                        marked.add(frozenset({p, q}))
                        unmarked.remove(frozenset({p, q}))
                        break
                if flag:
                    break

        # まず最小DFAの状態がどうなるか計算する
        states_dict = {}
        for p in self.states:
            states_dict[p] = {p}
        for p in self.states:
            for q in self.states:
                if frozenset({p, q}) in unmarked:
                    states_dict[p].add(q)
                    states_dict[q].add(p)

        # 次にstates_dictからtransitions, init_state, final_statesを作る
        states = set()
        init_state = None
        final_states = set()
        transitions = {}
        for p in self.states:
            state = frozenset(states_dict[p])
            if state in states:
                continue
            states.add(state)
            transitions[state] = {}
            for s in self.alphabets:
                # next(iter(X))は「Xから適当に一つ取る」という意味
                transitions[state][s] = states_dict[self.transitions[next(iter(state))][s]]
            if self.init_state in state:
                init_state = state
            if len(self.final_states.intersection(state)) != 0:
                final_states.add(state)
        alphabets = self.alphabets
        return DeterministicFiniteAutomata(states, alphabets, transitions, init_state, final_states)

以下のように使います。たしかに状態数が減っているのがわかります。

>>> from Automata import DeterministicFiniteAutomata as DFA
>>> a, b, c, d, e, f, g, h, = range(8)
>>> states = {a, b, c, d, e, f, g, h}
>>> alphabets = {"0", "1"}
>>> transitions = {
...     a: {"0": b, "1": f},
...     b: {"0": g, "1": c},
...     c: {"0": a, "1": c},
...     d: {"0": c, "1": g},
...     e: {"0": h, "1": f},
...     f: {"0": c, "1": g},
...     g: {"0": g, "1": e},
...     h: {"0": g, "1": c}
... }
>>> init_state = a
>>> final_states = {c}
>>> automata = DFA(states, alphabets, transitions, init_state, final_states)
>>> print(automata)

Automata
    states      : {0, 1, 2, 3, 4, 5, 6, 7}
    alphabets   : {'1', '0'}
    transitions : 0 : {'1': 5, '0': 1}
                  1 : {'1': 2, '0': 6}
                  2 : {'1': 2, '0': 0}
                  3 : {'1': 6, '0': 2}
                  4 : {'1': 5, '0': 7}
                  5 : {'1': 6, '0': 2}
                  6 : {'1': 4, '0': 6}
                  7 : {'1': 2, '0': 6}
    init_state  : 0
    final_states: {2}

>>> minimized = automata.minimized()
>>> print(minimized)

Automata
    states      : {{2}, {6}, {0, 4}, {1, 7}, {3, 5}}
    alphabets   : {'1', '0'}
    transitions : {6} : {'1': {0, 4}, '0': {6}}
                  {2} : {'1': {2}, '0': {0, 4}}
                  {0, 4} : {'1': {3, 5}, '0': {1, 7}}
                  {3, 5} : {'1': {6}, '0': {2}}
                  {1, 7} : {'1': {2}, '0': {6}}
    init_state  : {0, 4}
    final_states: {{2}}

感想

正規表現やりたいんだけどできるかな...