初めに
この記事は、2019/10/05 に行われた JuliaTokai #04 での発表 Julia で遊ぼう フィボナッチ数列 を補完するものです1。
主に、実験結果だけ示してほとんと説明できなかった後半の「2. 数列」を中心に、実験もし直したので新しい結果とともに紹介・解説をします。
Topics
-
○に対して○[n]のように書けるようにする
⇒Base.getindex(○, n::Integer)の（多重）定義 -
○に対して○[n:m]や○[n:s:m]のように書けるようにする
⇒Base.getindex(○, r::AbstractRange)の（多重）定義
＋ Iteration プロトコルの実装 -
○に対して○[n:end]や○[n:s:end]のように書けるようにする
⇒Base.lastindex(○)の（多重）定義
＋Base.:(:)(n, □)やBase.:(:)(n, m, □)の（多重）定義
＋Base.getindex(○, ▼)の（多重）定義
＋ Iteration プロトコルの実装（無限列挙）
0. フィボナッチ数列とは？
勉強会での発表時に参加者に聞いたら、フィボナッチ数列を知らな人が1人だけいました(!)。
ので念のため、この記事でも簡単に定義を書いておきます。
「大丈夫わかってる早く本題みたい」という方は 飛ばしてもらってOK です。
フィボナッチ数列の定義をざっくり書くと、以下のような三項間漸化式となります。
\begin{eqnarray*}
F_0 &=& 0 \\
F_1 &=& 1 \\
F_n &=& F_{n-1} + F_{n-2} \, (n \ge 2)
\end{eqnarray*}
つまり「直前の2つの数を足すと次の数になる数列」です。
最初の何件かを書き下すとこんな感じ：
0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, \dots
あと、フィボナッチ数列の一般項 $F_n$ のことを「フィボナッチ数」と呼びます。
詳細は、Wikipedia: フィボナッチ数 などを見てみてください。
1. 一般項を求めるアルゴリズム
なし崩し的ですがこちらは途中経過をすっ飛ばして↓を採用します。
function fib_n(n::Integer)
d = Dict(zero(n)=>big"0", one(n)=>big"1")
fib_n(n, d)
end
function fib_n(n, d)
if haskey(d, n)
return d[n]
end
if n < 0
result = iseven(n) ? -fib_n(-n, d) : fib_n(-n, d)
d[n] = result
return result
end
m = n ÷ 2
result = if iseven(n)
(2 * fib_n(m - 1, d) + fib_n(m, d)) * fib_n(m, d)
else
fib_n(m, d) ^ 2 + fib_n(m + 1, d) ^ 2
end
d[n] = result
return result
end
これはフィボナッチ数の「バイナリ公式」↓と「メモ化再帰」を利用した、一般の n 番目のフィボナッチ数を求める関数になっています。計算量オーダーは $O(\log n)$ です。
\begin{eqnarray*}
F_{2n} &=& (2F_{n-1} + F_n) F_n \\
F_{2n+1} &=& F_{n}^2 + F_{n+1}^2
\end{eqnarray*}
あと何気にマイナス方向にも拡張してたりします。数式で書くと↓のようになります。
F_{-n} = (-1)^{n+1} F_n
詳細は、英語版 Wikipedia の Fibonacci number-Sequence properties などを参照してください。
2. 数列
やりたいこと
こんな風に書けてこんな風に動くものを作ることを目指します。
-
Fib[n]でn番目のフィボナッチ数を取得できるようにする。- 例：
Fib[10] == 55
- 例：
-
Fib[n:m]でn〜m番目のフィボナッチ数列を取得できるようにする。- 例：
collect(Fib[-5:5]) == BigInt[5, -3, 2, -1, 1, 0, 1, 1, 2, 3, 5]
- 例：
-
Fib[n:end]でn番目以降のフィボナッチ数列を無限に取得できるようにする。- 例：
for v in Fib[5:end]; println(v); end # => 5, 8, 13, 21, 34, 55, 89, 144, …
- 例：
- ↑を
Fib[n:s:m]とかFib[n:s:end]とかでステップ指定して同様に処理できるように実装- 例：
sum(Fib[1:2:19]) == Fib[20] - 例：
collect(Iterators.take(Fib[5:-1:end], 10)) == BigInt[5, 3, 2, 1, 1, 0, 1, -1, 2, -3]
- 例：
Fib の定義
Fib は以下のようにしました。
struct FibType end
const Fib = FibType()
immutable な型 FibType を宣言し、そのインスタンスを const Fib として定義。なおメンバ何もないので FibType() インスタンスはシングルトンになります。関数で Fib を引数に受け取るときは単に ::FibType と書くだけでOKです（以降のコードを参照）。
Fib[n] の実装
○[〜] という記法（インデクシング）を実装するには、Base.getindex(○, 〜) を多重定義します2。
Fib[n] は「n 番目のフィボナッチ数を取得」という仕様なので、そのまま fib_n(n) を呼べばOK。簡単ですね！
Base.getindex(::FibType, index::Integer) = fib_n(index)
julia> Fib[10]
# => 55
julia> [Fib[n] for n=1:10]
# => BigInt[1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
julia> Fib[5000]
# => 3878968454388325633701916308325905312082127714646245(ry
$O(\log n)$ なので n=5000 とかでも一瞬で答えが返ってきます。1045桁もあるので途中で省略してあります。
Fib[n:m] の実装
Juliaでは、普通の配列でも a[1:5] と書くことで「1番目から5番目を取得」などということができます。
これを踏襲して、 Fib[n:m] と書くことで「n 〜 m 番目のフィボナッチ数列を取得」という仕様にしました。
fib_n(n) 〜 fib_n(m) を順に呼んで並べても良い（例：[fib_n(i) for i=n:m]）のですが、後々のことを考えて、ここでは「フィボナッチ数列の範囲を表すオブジェクト」を返すようにしてみます。
FibRange
struct FibRange{R<:AbstractRange{<:Integer}}
range::R
end
Base.getindex(::FibType, r::AbstractRange{<:Integer}) = FibRange(r)
getindex() の第2引数、r::AbstractRange{<:Integer} としていますが、範囲オブジェクトの要素の数値型は実質 Int しか来ない（コード中に普通に 1:10 とか書いたら isa(1:10, AbstractRange{Int}) == true になる）のですが、まぁ念のため。
とにかくこれで、「Fib[n:m] で 『n 〜 m 番目のフィボナッチ数列の範囲を表すオブジェクト』を取得」できます。
julia> Fib[1:10]
# => FibRange{UnitRange{Int64}}(1:10)
まだこのままだと数列の各項が取得できないので、Iteration プロトコル を実装します。
# IteratorEltype
Base.IteratorEltype(::Type{<:FibRange}) = Base.HasEltype()
# Base.eltype(fr::FibRange) = BigInt
Base.eltype(::Type{<:FibRange}) = BigInt
# IteratorSize
Base.IteratorSize(::Type{<:FibRange}) = Base.HasShape{1}()
Base.length(fr::FibRange) = length(fr.range)
Base.size(fr::FibRange) = size(fr.range)
# iterate
function Base.iterate(fr::FibRange{<:UnitRange})
index = fr.range.start
start = Fib[index]
prev = Fib[index - 1]
len = length(fr)
iterate(fr, (start, prev, len))
end
function Base.iterate(fr::FibRange{<:UnitRange}, (current, prev, len)::Tuple{BigInt, BigInt, Int})
len <= 0 && return nothing
(current, (current + prev, current, len - 1))
end
ちょっとだけ解説します。
-
IteratorEltype- イテレータの要素の型についての情報。対象の型を引数に受け取って、以下のいずれかを返す。
-
EltypeUnknown()：列挙する各要素の型が不定の場合。 -
HasEltype()：型が決まっている場合。この時Base.eltype()も実装する必要があります。↓参照
-
- デフォルトは
HasEltype()を返す仕様になっています。
- イテレータの要素の型についての情報。対象の型を引数に受け取って、以下のいずれかを返す。
-
::Type{<:FibRange}- 「引数として
FibRangeの subtype（派生型）を受け取る」という指定。FibRange{UnitRange{Int64}}もFibRange{StepRange{Int64,Int64}}も受け取れます。
- 「引数として
-
eltype- 実際に返ってくる各要素の型を返す。今回の場合、フィボナッチ数はすべて
BigIntで表現しているのでBigIntを返す。 -
eltype(x)は、xが型でない場合にeltype(typeof(x))に委譲されます（例：eltype([1,2,3]) == eltype(Array{Int,1}) == Int）。ので::Type{<:FibRange}を引数に受け取るものだけ実装すればOK。
逆にBase.eltype(fr::FibRange) = BigIntだけを実装すると後々問題が起きるので注意（詳細後述）
- 実際に返ってくる各要素の型を返す。今回の場合、フィボナッチ数はすべて
-
IteratorSize- イテレータのサイズについての情報。ここで言うサイズとは、列挙する長さや次元が決まっているか、ということ。対象の型を引数に受け取って、以下のいずれかを返す。
-
SizeUnknown()：列挙するサイズが不定の場合、「実行する度に長さが変わる可能性がある場合」など。 -
HasLength()：長さが有限で固定の場合。この時Base.length()も実装する必要があります。 -
HasShape{N}()：長さが分かっていて、かつ配列と互換性のある『次元数』の概念がある場合。この時Base.size()も実装する必要があります。 -
IsInfinite()：無限列挙する場合。
-
- デフォルトは
HasLength()を返す仕様になっています。
- イテレータのサイズについての情報。ここで言うサイズとは、列挙する長さや次元が決まっているか、ということ。対象の型を引数に受け取って、以下のいずれかを返す。
-
size- 各次元のサイズからなる tuple を返す。今回の場合
(length(fr),)という1要素からなる tuple を返せば良いけれど、中で持っているrangeのsizeと一致するのでそのような実装（Base.size(fr::FibRange) = size(fr.range)）にしています。
- 各次元のサイズからなる tuple を返す。今回の場合
-
iterate
結果：
julia> fib_1_10 = Fib[1:10]
# => FibRange{UnitRange{Int64}}(1:10)
julia> for v in fib_1_10
println(v)
end
# => 1
# 1
# 2
# 3
# 5
# 8
# 13
# 21
# 34
# 55
julia> collect(fib_1_10)
# => BigInt[1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
julia> sum(fib_1_10)
# => 143
julia> collect(Fib[-5:5])
# => BigInt[5, -3, 2, -1, 1, 0, 1, 1, 2, 3, 5]
パフォーマンス確認：
julia> collect(Fib[1:100]) == [Fib[n] for n=1:100]
# => true
julia> @time [Fib[n] for n=1:100];
# => 0.043813 seconds (57.95 k allocations: 2.765 MiB)
julia> @time collect(Fib[1:100]);
# => 0.000079 seconds (315 allocations: 7.023 KiB)
collect(Fib[1:100]) の方が速度もメモリ使用量も圧倒的に少ない！
Fib[n:s:m] の実装
Fib[n:s:m] は「n 〜 m 番目のフィボナッチ数を s個ごと（＝s-1個おき）に取得」という仕様。
実は先ほどの FibRange の定義で、そのようなオブジェクトを生成することはもうできてます。Base.eltype() も Base.size() もそのまま使用できるので、後は iterate() を実装するだけ。
ここで先ほどの FibRange{<:UnitRange} の時は、フィボナッチ数の定義 $F_{n}=F_{n-1}+F_{n-2}$ を利用していました。今回はそれを s 回繰り返すことでも次の要素を算出することは可能ですが、もう一工夫して以下の「和の公式」を利用します。
\begin{eqnarray*}
F_{n+m−1} &=& F_nF_m + F_{n−1}F_{m−1} \\
F_{n+m} &=& F_n F_m + F_{n−1} F_m + F_n F_{m−1}
\end{eqnarray*}
具体的には以下のようになります：
function Base.iterate(fr::FibRange{<:StepRange})
n = fr.range.start
m = oftype(n, fr.range.step)
d = Dict(zero(n)=>big"0", one(n)=>big"1")
fₙ = fib_n(n, d)
fₙ₋₁ = fib_n(n - one(n), d)
fₘ = fib_n(m, d)
fₘ₋₁ = fib_n(m - one(n), d)
len = length(fr)
iterate(fr, (fₙ, fₙ₋₁, fₘ, fₘ₋₁, len))
end
function Base.iterate(fr::FibRange{<:StepRange}, (fₙ, fₙ₋₁, fₘ, fₘ₋₁, len)::Tuple{BigInt, BigInt, BigInt, BigInt, Int})
len <= 0 && return nothing
fₙ₊ₘ = fₙ * fₘ + fₙ₋₁ * fₘ + fₙ * fₘ₋₁
fₙ₊ₘ₋₁ = fₙ * fₘ + fₙ₋₁ * fₘ₋₁
(fₙ, (fₙ₊ₘ, fₙ₊ₘ₋₁, fₘ, fₘ₋₁, len - 1))
end
公式と、先ほどと共通の部分以外で少しだけ追加解説：
-
fr::FibRange{<:StepRange}- 先程は
fr::FibRange{<:UnitRange}でした。isa(Fib[1:10], FibRange{<:UnitRange}) == trueです。今回はisa(Fib[0:2:20], FibRange{<:StepRange}) == trueとなり、よってfr::FibRange{<:StepRange}で引数を受け取れます。
これ何気に便利。つまり「型パラメータの差異」で実装を変えることができるのです。これはJuliaの多重ディスパッチの大きな特徴の一つ。
- 先程は
-
oftype(n, m)-
T = typeof(n)として、convert(T, m)することと同じ。
StepRangeは仕様上はstartとstepの型は同じとは限らないので、ここで合わせているわけです。実質0:2:20と書けばisa(0:2:20, StepRange{Int,Int}) == trueなのでたいていの場合startとstepの型は一致するのですが、念のため。
-
-
fₙとかfₘ₋₁とか- Julia は変数名にUnicode文字列を自由に使えます。定義されていればこのような所謂「下付き文字」も使えます。数式をそのままコードで表現できるので、かなり便利。
なおJuliaのREPLや対応するエディタなら、例えば\_n+Tab で下付き文字のₙに変換・入力ができます。
- Julia は変数名にUnicode文字列を自由に使えます。定義されていればこのような所謂「下付き文字」も使えます。数式をそのままコードで表現できるので、かなり便利。
確認
julia> fr2 = Fib[0:2:20]
# => FibRange{StepRange{Int64,Int64}}(0:2:20)
julia> for v in fr2
println(v)
end
# => 0
# 1
# 3
# 8
# 21
# 55
# 144
# 377
# 987
# 2584
# 6765
julia> collect(Fib[10:-1:1])
# => BigInt[55, 34, 21, 13, 8, 5, 3, 2, 1, 1]
julia> sum(Fib[1:2:19]) == Fib[20] # F_1 + F_3 + … + F_19 == F_20
# => true
julia> collect(Fib[0:10:100])
# => BigInt[ 0,
# 55,
# 6765,
# 832040,
# 102334155,
# 12586269025,
# 1548008755920,
# 190392490709135,
# 23416728348467685,
# 2880067194370816120,
# 354224848179261915075]
パフォーマンス確認：
julia> collect(Fib[0:10:1000]) == [Fib[n] for n=0:10:1000]
# => true
julia> @time [Fib[n] for n=0:10:1000];
# => 0.038576 seconds (65.97 k allocations: 2.921 MiB)
julia> @time collect(Fib[0:10:1000]);
# => 0.000831 seconds (2.38 k allocations: 70.695 KiB)
良い感じですね♪
Fib[n:end]、Fib[n:s:end] の実装準備
Julia では、普通の配列で a[n:end] と書くことで「n番目以降最後まで取得」などということができます5。
フィボナッチ数列は無限数列なので、Fib[n:end] または Fib[n:s:end] とすると、終わりなく無限に列挙できると良いですよね！
無限列挙できるようにしましょう！
ポイントは以下：
-
○[△:end]はgetindex(○, △:lastindex(○))に変換されて評価される
なので、これを実現するのは以下の4ステップとなります：
-
Base.lastindex(○)を実装（多重定義）して適切なオブジェクト（＝□）を返すようにする -
△:□を実装（多重定義）して適切なオブジェクト（＝▼）を返すようにする -
Base.getindex(○, ▼)を実装（多重定義）して適切なオブジェクト（無限列）（＝◎）を返すようにする -
Base.iterate(◎, 〜)を適切に実装（多重定義）して無限列挙できるようにする
1. Base.lastindex(::FibType) の実装
通常の配列や文字列なら、lastindex は文字通り「最後のインデックス番号」返せばOKです。
でもフィボナッチ数列は無限数列なので、「最後」はありません。
そこで、特殊な型を定義してその型のインスタンスを返すようにしましょう。
具体的には、以下のようにしました：
struct EndOfFib end
Base.lastindex(::FibType) = EndOfFib()
つまり仮想的な「Fibの終端」オブジェクトを生成してそれを返すようにするだけ。
簡単ですね！
2. n:(::EndOfFib) の実装
n:m とか n:s:m は通常、Rangeオブジェクト（UnitRange/StepRange）を返します。
これは : 演算子がそのように実装されているから、なのです。
今回は通常の Rangeオブジェクト を返してもらっては困るので、やはり新しい型を定義してそれを返してもらうようにします。
具体的にはこんな感じ：
abstract type AbstractSequence{I<:Integer} end
struct UnitSequence{I<:Integer} <: AbstractSequence{I}
start::I
end
struct StepSequence{I<:Integer} <: AbstractSequence{I}
start::I
step::I
end
Base.:(:)(start::Integer, ::EndOfFib) = UnitSequence(start)
Base.:(:)(start::I, step::I, ::EndOfFib) where {I<:Integer} = StepSequence(start, step)
Base.:(:)(start::Integer, step::Integer, ::EndOfFib) = StepSequence(promote(start, step)...)
ちょっとだけややこしいので、少し解説。
-
AbstractSequence/UnitSequence/StepSequence-
AbstractRange/UnitRange/StepRangeに倣って似たような型を定義してみました。それらとの違いは、「終端」の情報がない、いうこと。
-
-
Base.:(:)(〜, ::EndOfFib)-
:演算子の多重定義。特に最後がEndOfFibとなるような場合についての定義。
引数2つの場合はUnitSequenceを、3つの場合はStepSequenceを返しています。 - ちなみに
:演算子の多重定義の書き方は他にもあります。
事前にimport Base: (:)と1行書いておけば、(:)(start, stop) = 〜のように書くこともできます。
またtypeof(:) == Colonなのですがそれを利用して(::Colon)(start, stop) = 〜とも書けます。6
-
-
promote(start, step)...-
startとstepが異なる型（だけどどちらもIntegerの subtype）の場合に、より上位の型（正確にはよりビット幅の大きい方の型）に変換（昇格）して型を揃えるのがpromote()関数の役目。その戻り値はtupleなので、...で引数展開しています。
-
3. Base.getindex(::FibType, ::AbstractSequence) の実装
やっとここまで来ました。Fib[n:end] や Fib[n:s:end] でフィボナッチ数列の無限列を取得できるようにしましょう。やはり新しい型を定義して対応します：
struct FibSequence{S<:AbstractSequence}
sequence::S
end
Base.getindex(::FibType, s::AbstractSequence) = FibSequence(s)
Base.IteratorEltype(::Type{<:FibSequence}) = Base.HasEltype()
# Base.eltype(fr::FibSequence) = BigInt
Base.eltype(::Type{<:FibSequence}) = BigInt
Base.IteratorSize(::Type{<:FibSequence}) = Base.IsInfinite()
まずはここまで。さきほどの FibRange とだいたい同じ。違うのは最後の Base.IteratorSize(::Type{<:FibSequence}) = Base.IsInfinite() の部分。先程は Base.HasShape{1}() を返していましたが、今回は無限列挙なので Base.IsInfinite() を返します。あと無限列なので Base.length() も Base.size() も定義（実装）していません（むしろ定義してはいけない）。
4. Base.iterate(::FibSequence) の実装
続いて iterate()：
function Base.iterate(fr::FibSequence{<:UnitSequence})
index = fr.sequence.start
start = Fib[index]
prev = Fib[index - 1]
iterate(fr, (start, prev))
end
function Base.iterate(fr::FibSequence{<:UnitSequence}, (current, prev)::Tuple{BigInt, BigInt})
(current, (current + prev, current))
end
function Base.iterate(fr::FibSequence{<:StepSequence})
n = fr.sequence.start
m = oftype(n, fr.sequence.step)
d = Dict(zero(n)=>big"0", one(n)=>big"1")
fₙ = fib_n(n, d)
fₙ₋₁ = fib_n(n - one(n), d)
fₘ = fib_n(m, d)
fₘ₋₁ = fib_n(m - one(n), d)
iterate(fr, (fₙ, fₙ₋₁, fₘ, fₘ₋₁))
end
function Base.iterate(fr::FibSequence{<:StepSequence}, (fₙ, fₙ₋₁, fₘ, fₘ₋₁)::Tuple{BigInt, BigInt, BigInt, BigInt})
fₙ₊ₘ = fₙ * fₘ + fₙ₋₁ * fₘ + fₙ * fₘ₋₁
fₙ₊ₘ₋₁ = fₙ * fₘ + fₙ₋₁ * fₘ₋₁
(fₙ, (fₙ₊ₘ, fₙ₊ₘ₋₁, fₘ, fₘ₋₁))
end
こちらも実は FibRange とほとんど同じ。違うのは、長さのチェックをしていないこと。
len <= 0 && return nothing がない（いらない）のです。これをなくすだけで、終了条件チェックがなくなる＝ 無限列挙 になります。
簡単ですね！
確認
動作確認してみましょう。
julia> fibseq1 = Fib[0:end]
# => FibSequence{UnitSequence{Int64}}(UnitSequence{Int64}(0))
julia> for v in fibseq1
println(v)
end
# => 0
# 1
# 1
# 2
# 3
# 5
# 8
# 13
# 21
# 34
# 55
# 89
# 144
# 233
# 377
# 610
# : ※適当なところで止めてください！
julia> collect(Iterators.take(fibseq1, 21))
# => BigInt[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]
julia> collect(Iterators.take(Fib[0:-1:end], 21))
# => BigInt[0, 1, -1, 2, -3, 5, -8, 13, -21, 34, -55, 89, -144, 233, -377, 610, -987, 1597, -2584, 4181, -6765]
まとめ？
実験はここまでですが、せっかく無限列挙できるようにしたので、それを利用してなにか追加実験もしてみたいですね。
あと逆に、無限じゃないほう（FibRange）は AbstractArray の subtype にできる（UnitRange や StepRange のように）ので、Abstract Arrays に示されている他の関数を実装したり、Broadcasting に対応させたり、まだいろいろ遊べますね。
また余裕ができたら追加実験の結果を記事にしたいと思います。
参考
-
いろいろあって2週間空いてしまいました…。 ↩
-
さらに
○[〜] = ×という代入を実装したい場合はBase.setindex!(○, ×, 〜)を実装する必要があります。今回は代入はないのでこちらには触れません。 ↩ -
ちなみに「
for文（のin節）に渡す」ことだけが目的なら、このiterate()だけを実装（多重定義）すればOKです。IteratorEltype()とかIteratorSize()とかは、collect()とかIterators.take()などに渡す場合に必要になってくるものです。 ↩ -
Julia のイテレーションについては、過去にも記事にしています（→ Julia 0.7-DEV の新しい Iteration に触れてみた。）。こちらも参照してみてください。 ↩
-
Julia の
endは予約語です。ブロックの終端と、インデックスの終端の、2つの使い途があります。 ↩ -
あと実は
n:(::EndOfFib) = 〜とかn:s:(::EndOfFib) = 〜みたいに書いても定義できちゃいます。できるんですけれどなんか逆に難易度高いので（解説めんどいので）興味のある人はご自分で試してみてください。 ↩