読者です 読者をやめる 読者になる 読者になる

yukiのブログ

このブログの記事はhttps://yuki67.github.io/からコピーしたものです。

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乗加速法を使って円周率を求める話があるのですがこの記事には書きませんでした。

疲れた。