Julia で遊ぼうフィボナッチ数列

Julia で遊ぼう フィボナッチ数列 ===   2019/10/05 JuliaTokai \#04 ＠道の駅とよはし antimon2（後藤俊介） Note: Julia でフィボナッチ数と戯れようっ！ ----  ## お品書き + お前誰よ？ + 簡単なJuliaの紹介 + Julia でフィボナッチ ---  # お前誰よ？ ----  ## 自己紹介 + 名前：後藤俊介 + 所属：**[有限会社来栖川電算](https://www.kurusugawa.jp)** + コミュニティ：**[JuliaTokai](https://juliatokai.connpass.com/)**, **[機械学習名古屋](https://machine-learning.connpass.com/)**, [Python東海](https://connpass.com/series/292/), Ruby東海, … + 言語：**[Julia](https://julialang.org)**, Python, Scala（勉強中）, Ruby, … + ![Twitter](https://i.imgur.com/HqouMIg.png) [@antimon2](https://twitter.com/antimon2) / ![Facebook](https://i.imgur.com/01nPd37.png) [antimon2](https://www.facebook.com/antimon2) + ![Github](https://i.imgur.com/yBKtii5.png) [antimon2](https://github.com/antimon2/) / ![Qiita](https://i.imgur.com/FxHMi64.png) [@antimon2](http://qiita.com/antimon2) / [ @antimon2](https://hackmd.io/@antimon2) Note: 今日は(も？)大っぴらに Julia の話ができるっ ----  [![有限会社来栖川電算](https://i.imgur.com/8Kuhfel.png) https://www.kurusugawa.jp](https://www.kurusugawa.jp) Note: 普段は Python で仕事してますっ ----  [![Python東海第40回勉強会](https://i.imgur.com/aRbEU0t.png) https://connpass.com/event/144750/](https://connpass.com/event/144750/) Note: ちょうど1週間後に開催！ロングLTやりますっ ----  [![機械学習名古屋第22回勉強会 AnnoFab ハンズオン](https://i.imgur.com/YiPWCbQ.png) https://machine-learning.connpass.com/event/148266/](https://machine-learning.connpass.com/event/148266/) Note: 来月開催！ ---  # 簡単な Julia の紹介 Note: Julia の紹介っ ----  [![Julia](https://raw.githubusercontent.com/JuliaLang/julia-logo-graphics/master/images/julia-logo-color.svg?sanitize=true)](https://julialang.org) Note: 実は何気にロゴマイナーチェンジしてますっ ----  ## Julia とは？(1) + [The Julia Language](https://julialang.org) + 最新 v1.2.0（2019/08/20） + LTS v1.0.5（2019/09/09） + pre v1.3.0-rc3（2019/10/04） + 科学技術計算に強い！ + 動作が速い！（LLVM JIT コンパイル） Note: ググるときはなるべく [julialang](https://www.google.co.jp/search?q=julialang) で！ ----  ## Julia とは？(2) > + Rのように中身がぐちゃぐちゃでなく、 > + Rubyのように遅くなく、 > + Lispのように原始的またはエレファントでなく、 > + Prologのように変態的なところはなく、 > + Javaのように硬すぎることはなく、 > + Haskellのように抽象的すぎない > > ほどよい言語である  引用元：http://www.slideshare.net/Nikoriks/julia-28059489/8  ----  ## Julia とは？(3) > + C のように高速だけど、 Ruby のような動的型付言語である > + Lisp のようにプログラムと同等に扱えるマクロがあって、しかも Matlab のような直感的な数式表現もできる > + Python のように総合的なプログラミングができて、 R のように統計処理も得意で、 Perl のように文字列処理もできて、 Matlab のように線形代数もできて、 shell のように複数のプログラムを組み合わせることもできる > + 超初心者にも習得は簡単で、超上級者の満足にも応えられる > + インタラクティブにも動作して、コンパイルもできる  （[Why We Created Julia](http://julialang.org/blog/2012/02/why-we-created-julia) から抜粋・私訳）  Note: いろんな言語の「いいとこどり」言語！ってことでっ ----  ## 要するに  + 動的言語なのに速い！ + 文法も覚えやすい！ + 数値計算に強い！  Note: 機械学習とかにも持って来いっ！ ----  ## 主な機能  + [多重ディスパッチ](https://docs.julialang.org/en/v1/manual/methods/) + [動的型システム](https://docs.julialang.org/en/v1/manual/types/) + [並行・並列処理](https://docs.julialang.org/en/v1/manual/parallel-computing/)、コルーチン + [組込パッケージマネージャ](https://docs.julialang.org/en/v1/stdlib/Pkg/)  Note: っ ----  ## 文法・関数 Note: 以降、ほぼ過去スライドからのコピペ。すっ飛ばして先へ進んで戴いてもOKっ ----  ### 基本的な演算 ```julia julia> 1 + 2 - 3 * 4 # 四則演算（除算以外） -9 julia> 7 / 5 # `整数 / 整数` の結果は浮動小数 1.4 julia> 7 ÷ 5 # `整数 ÷ 整数` の結果は整数 1 julia> 2 ^ 10 # 冪乗は `^` 1024 julia> 123 & 234 | 345 # 論理積 / 論理和 376 julia> 123 ⊻ 234 # 排他的論理和(==`xor(123, 234)`) 145 ```  Note: 整数同士の除算は実数になりますっ整数除算演算子 `÷` が別に存在します（Python の `//` 相当）っまた冪乗も（`**` ではなく）`^` ですっ `⊻` は `\xor`+<kbd>Tab</kbd> または `\veebar`+<kbd>Tab</kbd> で変換できますっちなみに先ほどの `÷` も `\div`+<kbd>Tab</kbd>で（基本的に ${\rm \TeX}$ の書式）っ ----  ### 配列 ```julia julia> a = [1, 2, 3, 4, 5] 5-element Array{Int64,1}: 1 2 3 4 5 julia> a[1] # Julia は 1-origin 1 julia> println(a[2:3]) # 範囲指定は両端含む [2, 3] ```  Note: 1-origin であることに注意すればあとは普通の配列っあと `a:b` は範囲（`Range`）の記法。両端を含む（Ruby の `a..b` と同じ）っ ----  ### 配列の内包表記 (1) ```julia julia> a = [n^2 for n=1:5] 5-element Array{Int64,1}: 1 4 9 16 25 julia> A = [x+10y for y=1:3, x=1:3] 3×3 Array{Int64,2}: 11 12 13 21 22 23 31 32 33 ```  Note: 内包表記の記法は Python に類似っかつ、`for` にカンマ区切りで複数のイテレータを渡すことで2次元以上の配列も作成可能っ ----  ### 配列の内包表記 (2) ```julia julia> [(a,b,c) for c=1:15,b=1:15,a=1:15 if a^2+a*b+b^2==c^2] 6-element Array{Tuple{Int64,Int64,Int64},1}: (3, 5, 7) (5, 3, 7) (6, 10, 14) (7, 8, 13) (8, 7, 13) (10, 6, 14) ``` Note: Python と同様に `if` で条件を指定することも可能っあと Python と同様、`[○ for ○=○]` を `(○ for ○=○)` と書くと配列ではなくて `Generator` が返りますっ ----  ### ベクトル ```julia julia> x = [1., 2., 3.]; y = [3., 1., 2.]; julia> x + y # `x .+ y` と書いても同じ（elementwise operation） [4., 3., 5.] julia> x .* y # これは `x * y` と書くとNG [3., 2., 6.] julia> using LinearAlgebra julia> x ⋅ y # 内積（dot積、`dot(x, y)` と書いても同じ） 11.0 julia> x × y # 外積（cross積、`cross(x, y)` と書いても同じ） [1., 7., -5.] ```  Note: Julia では実は1次元配列がベクトルの扱いっ `⋅` は `\cdot`+<kbd>Tab</kbd>、`×` は `\times`+<kbd>Tab</kbd>（これらを利用するには `using LinearAlgebra` 必要）っあとこれらや先ほどの `÷` や `⊻` などのように、ASCIIの範囲を超えたUnicode文字の演算子（そのほとんどが $\TeX$ 由来）が Julia にはたくさんあります（他には例えば比較演算子の `≤` `≥` や、集合の要素 `∈` や包含関係 `⊆` などなど） ----  ### 行列 ```julia julia> A = [1 2; 3 4] # この記法は MATLAB/Octave 由来 2×2 Array{Int64,2}: 1 2 3 4 julia> A' # `○'` は転置行列の記法（これも MATLAB/Octave 由来） 2×2 LinearAlgebra.Adjoint{Int64,Array{Int64,2}}: 1 3 2 4 julia> transpose(A) # 正確には転置行列はこっち 2×2 LinearAlgebra.Transpose{Int64,Array{Int64,2}}: 1 3 2 4 ```  Note: Julia では2次元配列が行列の扱いっあと `○.'` という書式は廃止されました（`transpose(A)` 使ってね）っ ----  ### 行列の演算 ```julia julia> A = [1 2; 3 4]; B = [3 0; 0 6]; julia> A + B # A .+ B でも同様 2×2 Array{Int64,2}: 4 2 3 10 julia> A * B # matrix multiply 2×2 Array{Int64,2}: 3 12 9 24 julia> A .* B # elementwise multiply 2×2 Array{Int64,2}: 3 0 0 24 ```  Note: 行列は `*` で通常の行列積になりますっこれ便利っ ----  ### ブロードキャスト ```julia julia> sin(0.1) 0.09983341664682815 julia> sin.([0.1, 0.2, 0.3, 0.4]) 4-element Array{Float64,1}: 0.0998334 0.198669 0.29552 0.389418 julia> [0.1, 0.2, 0.3, 0.4] .^ 2 # => [0.01, 0.04, 0.09, 0.16] ```  Note: 関数名と `(` の間に `.` を置くと、普通の関数を配列に拡張してくれる（ブロードキャスト）っ `.^` のように演算子の前に `.` を書いても同様（先ほど出た `.+` `.*` もブロードキャスト）っ ----  ### 関数定義 ```julia julia> f(x) = x^2 + 2x - 1 f (generic function with 1 method) julia> f(1) 2 julia> f.(1:5) # => [2, 7, 14, 23, 34] ``` Note: 数学のように直感的な記述で関数を定義可能っ `2x` は `2*x` の省略形、曖昧さがなければリテラルと他の識別子が続く場合などに勝手に乗算と解釈してくれるっまたユーザ定義関数も `.` をつけて自動的にブロードキャスト対応っ ----  ### 有理数・複素数 ```julia julia> 1//2 == 0.5 true julia> 1//2 - 1//3 1//6 julia> 1im ^ 2 == -1 true julia> (1.0 + 0.5im) * (2.0 - 3.0im) 3.5 - 2.0im ``` Note: 有理数・複素数を標準サポート。 `//` は有理数除算（結果は有理数） `im` は虚数単位。どちらも四則演算も普通に書けますっ ---  # Julia で フィボナッチ Note: みんな大好きフィボナッチ数列っ ----  ## フィボナッチ数参考：[Wikipedia: フィボナッチ数](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A3%E3%83%9C%E3%83%8A%E3%83%83%E3%83%81%E6%95%B0) Note: 知らない方のために説明っ ---  ## 1. アルゴリズム Note: まずは一般項を求めてみよう ----  $$ F_n = \begin{cases} 0 & (n = 0)\\ 1 & (n = 1)\\ F_{n-2} + F_{n-1} & (n \ge 2) \end{cases} $$ Note: とりあえずよくある定義っ ----  ### 0. ナイーブ実装 Note: よくある定義の通りに実装っ ----  ```julia function fib0(n::Integer) if n < 2 return big(n) end return fib0(n - 2) + fib0(n - 1) end ```  Note: 簡単のために `n ≥ 0` と仮定したコードになっていますっあと `n::Integer` で引数を受けているので、整数以外はエラーになりますっあとあと、戻り値の型はどこにも明示してないですが、型推論により `BigInt` が返るようになっていますっ ----  ```julia julia> [fib0(n) for n=0:10] # => [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55] julia> @time fib0(20) # => 0.001952 seconds (54.73 k allocations: 941.000 KiB) # => 6765 julia> @time fib0(100) # => ... ``` Note: 100番目のフィボナッチ数求めようとすると全然返ってきません :persevere: ----  `fib0` $\in O(2^n)$ Note: つまり指数オーダーっ ----  ### 1. メモ化 Note: 厳密なメモ化ではなくてキャッシュ化とでも言うべき？ ----  ```julia function fib1(n::Integer) d = Dict(zero(n)=>big"0", one(n)=>big"1") fib1(n, d) end function fib1(n, d) if haskey(d, n) return d[n] end result = fib1(n - 2, d) + fib1(n - 1, d) d[n] = result return result end ```  Note: っ ----  ```julia julia> [fib1(n) for n=0:10] # => [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55] julia> @time fib1(20) # => 0.000026 seconds (75 allocations: 3.391 KiB) # => 6765 julia> @time fib1(100) # => 0.000024 seconds (311 allocations: 11.422 KiB) # => 354224848179261915075 ``` Note: 100番目のフィボナッチ数も割と一瞬で返ってくるっ ----  ```julia julia> @time fib1(1000) # => 0.000248 seconds (3.02 k allocations: 177.641 KiB) # => 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875 julia> @time fib1(5000) # => 0.005538 seconds (15.03 k allocations: 1.600 MiB, 61.86% gc time) # => 38789684543883256337019163083259053120821277146462451061605972148955501390440370970108229164622106694792934528588829738134831020089549829403614301569114789383642165639441069102145056341337065586562382546567… ``` Note: $F_{5000}$ の値はばかでかいので省略っあとなんとなく比例してそうなのは見て取れるかな？ ----  `fib1` $\in O(n)$ Note: 理論上はっもっと厳密には、辞書への登録と読み出しも関わってくるのでもう少し変わってくるっ ----  ### 2. 末尾再帰 Note: メモ化を使わずに $O(n)$ を実現するテクニック ----  $$ \begin{eqnarray*} F_0 + F_1 &=& F_2 \\ F_1 + F_2 &=& F_3 \\ &\vdots& \\ F_{n-2} + F_{n-1} &=& F_n \end{eqnarray*} $$ Note: 考え方：これを順々に実行っ ----  ```julia function fib2(n::Integer) if n < 2 return big(n) end fib2(n, big"0", big"1") end function fib2(n, a, b) if n == 1 return b end fib2(n - 1, b, a + b) end ```  Note: っ ----  ```julia julia> [fib2(n) for n=0:10] # => [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55] julia> @time fib2(20) # => 0.000017 seconds (64 allocations: 1.266 KiB) # => 6765 julia> @time fib2(100) # => 0.000024 seconds (301 allocations: 4.844 KiB) # => 354224848179261915075 ``` Note: 100番目のフィボナッチ数も割と一瞬で返ってくるっ ----  ```julia julia> @time fib2(1000) # => 0.000249 seconds (3.00 k allocations: 85.328 KiB) # => 43466557(ry julia> @time fib2(5000) # => 0.001513 seconds (15.00 k allocations: 1.244 MiB) # => 38789684(ry ``` Note: やっぱりだいたい比例してそうっあとメモリ消費量が少ない分実は高速だったりもするっ ----  `fib2` $\in O(n)$ Note: やはり理論上っ ----  ### 3. 行列計算利用 Note: いよいよここから $O(\log n)$ の世界っ ----  $$ \begin{pmatrix} F_{n+1} & F_{n} \\ F_{n} & F_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} ^ n $$ 参考：[Wikipedia: フィボナッチ数](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A3%E3%83%9C%E3%83%8A%E3%83%83%E3%83%81%E6%95%B0)  Note: こんな公式があるのさっ [Wikipedia: フィボナッチ数](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A3%E3%83%9C%E3%83%8A%E3%83%83%E3%83%81%E6%95%B0) にも載ってます。 ----  ```julia function fib3(n::Integer) if n < 2 return big(n) end F₁ = BigInt[1 1 1 0] Fₙ = F₁ ^ n return Fₙ[2, 1] # `Fₙ[1, 2]` でもOK end ```  Note: `F₁` が行列になっている。それを `^ n` で簡単に n乗できるっ（裏で BLAS が動きますたぶん） ----  ```julia julia> [fib3(n) for n=0:10] # => [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55] julia> @time fib3(20) # => 0.000028 seconds (190 allocations: 3.820 KiB) # => 6765 julia> @time fib3(100) # => 0.000025 seconds (298 allocations: 5.648 KiB) # => 354224848179261915075 ``` Note: ここまでは `fib2` までとあんま変わらないっ ----  ```julia julia> @time fib3(1000) # => 0.000053 seconds (520 allocations: 12.023 KiB) # => 43466557(ry julia> @time fib3(5000) # => 0.000060 seconds (594 allocations: 25.398 KiB) # => 38789684(ry ``` Note: 速度もメモリ消費量もそんなに増えてないことに注目っ ----  `fib3` $\in O(\log n)$ Note: これは行列の累乗が内部で $O(\log n)$ となるように実装されているからっ ----  ### 4. バイナリ公式 Note: ＋メモ化っ ----  $$ \begin{eqnarray*} F_{2n} &=& (2F_{n-1} + F_n) F_n \\ F_{2n+1} &=& F_{n}^2 + F_{n+1}^2 \end{eqnarray*} $$ Note: こんな公式もあるのさってかさっきの行列の累乗と $F_n$ の関係が証明できればそこから簡単に導き出せます。 ----  ```julia function fib4(n::Integer) d = Dict(zero(n)=>big"0", one(n)=>big"1") fib4(n, d) end function fib4(n, d) if haskey(d, n) return d[n] end m = n ÷ 2 result = if iseven(n) (2 * fib4(m - 1, d) + fib4(m, d)) * fib4(m, d) else fib4(m, d) ^ 2 + fib4(m + 1, d) ^ 2 end d[n] = result return result end ```  Note: ちょっと複雑になってきましたっ ----  ```julia julia> [fib4(n) for n=0:10] # => [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55] julia> @time fib4(20) # => 0.000021 seconds (64 allocations: 1.805 KiB) # => 6765 julia> @time fib4(100) # => 0.000017 seconds (127 allocations: 3.945 KiB) # => 354224848179261915075 ``` Note: ここまでは `fib3` と比較してメモリ消費量半分くらいっ ----  ```julia julia> @time fib4(1000) # => 0.000027 seconds (202 allocations: 5.555 KiB) # => 43466557(ry julia> @time fib4(5000) # => 0.000038 seconds (262 allocations: 8.867 KiB) # => 38789684(ry ``` Note: 速度もメモリ消費量もそんなに増えてないことに注目っ ----  `fib4` $\in O(\log n)$ Note: これでも十分良い感じっ ----  ### 5. 和の公式＋末尾再帰 Note: 最強っ ----  $$ \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} \\ (F_{2^k+m} &=& F_{2^k}F_m + F_{2^k−1}F_m + F_{2^k}F_{m−1}) \end{eqnarray*} $$ Note: こんな公式もあるのさってかやはりさっきの行列の累乗と $F_n$ の関係から簡単に導き出せます。 ----  ```julia function fib5(n::Integer) if n < 2 return big(n) end fib5(n, big"0", big"1", big"0", big"1") end function fib5(n, a, b, p, q) if n == 1 return a * p + b * q end if isodd(n) a, b = a * p + b * q, a * q + b * q + b * p end fib5(n ÷ 2, a, b, p ^ 2 + q ^ 2, (2p + q) * q) end ```  Note: かなり複雑になってきましたっどうしてあの公式からこんなコードが生まれるのか、は読者の宿題にいたしますっw ----  ```julia julia> [fib5(n) for n=0:10] # => [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55] julia> @time fib5(20) # => 0.000021 seconds (96 allocations: 1.797 KiB) # => 6765 julia> @time fib5(100) # => 0.000015 seconds (150 allocations: 2.477 KiB) # => 354224848179261915075 ``` Note: ここまでは `fib4` とほぼ同速度でメモリ消費量ちょっと多めっ ----  ```julia julia> @time fib5(1000) # => 0.000028 seconds (267 allocations: 5.219 KiB) # => 43466557(ry julia> @time fib5(5000) # => 0.000036 seconds (300 allocations: 10.164 KiB) # => 38789684(ry ``` Note: っ ----  `fib5` $\in O(\log n)$ Note: こちらもまぁまぁ良い感じっだけど結果 `fib4` の方が良さげ ---  ## 2. 数列 Note: 型を定義して数列を表現しよう！ ----  ### やりたいこと Note: こんな風に書けるようにしてみよう ----  ```julia julia> Fib[10] # => 55 julia> collect(Fib[-5:5]) # => BigInt[5, -3, 2, -1, 1, 0, 1, 1, 2, 3, 5] julia> sum(Fib[1:2:19]) == Fib[20] # => true ``` Note: 最後のは、 $F_1 + F_3 + \dots + F_{2n-1} = F_{2n}$ という公式の確認（n=10）。 ----  {%gist antimon2/578f2417f8a3310d94eedb459721d422%} [FibType.jl.ipynb (nbviewer)](https://nbviewer.jupyter.org/gist/antimon2/00f9c2bd8e1bc30b3217bd4234973ac5/FibType.jl.ipynb) Note: 時間なかったので gist そのまま貼っときますｗ ---  # まとめ ----  + 関数と型を使いこなせ！ + 多重ディスパッチ使いこなせ！ + Julia 楽しい！ Note: っ ---  # リンク ----  ## 実験Note + [Fibonacci.jl.ipynb](https://gist.github.com/antimon2/00f9c2bd8e1bc30b3217bd4234973ac5#file-fibonacci-jl-ipynb) ([nbviewer](https://nbviewer.jupyter.org/gist/antimon2/00f9c2bd8e1bc30b3217bd4234973ac5/Fibonacci.jl.ipynb)) + [FibType.jl.ipynb](https://gist.github.com/antimon2/00f9c2bd8e1bc30b3217bd4234973ac5#file-fibtype-jl-ipynb) ([nbviewer](https://nbviewer.jupyter.org/gist/antimon2/00f9c2bd8e1bc30b3217bd4234973ac5/FibType.jl.ipynb)) ----  ## 参考リンク + [Wikipedia: フィボナッチ数](https://ja.wikipedia.org/wiki/%E3%83%95%E3%82%A3%E3%83%9C%E3%83%8A%E3%83%83%E3%83%81%E6%95%B0) + [Collections and Data Structures](https://docs.julialang.org/en/v1/base/collections/) - [Julia Documentation](https://docs.julialang.org/en/v1) ---  # おまけ ----  ## 機械学習名古屋 ----  [![機械学習名古屋第22回勉強会 AnnoFab ハンズオン](https://i.imgur.com/YiPWCbQ.png) https://machine-learning.connpass.com/event/148266/](https://machine-learning.connpass.com/event/148266/) Note: 来月開催！（再） ---  ご清聴ありがとうございます。 Note: ご清聴ありがとうございますっ！