Julia が速すぎるお話．みんなで挑戦しよう！

# 本日は

先日 [速度でJuliaに勝つのが難しかった（というか無理な）話．](https://qiita.com/SatoshiTerasaki/items/5bd3ee78a60f9a02e0fc)
という記事を書きましたが，多くの方に読んていただいたみたいです．ありがとうございます．

Twitter で Gen Kuroki さんがJuliaネタ，数学ネタを投下していましたので，このネタでも単純ですが速度で挑戦してみましょう．

# 元ネタはこちらから

<blockquote class="twitter-tweet" data-lang="ja"><p lang="ja" dir="ltr"><a href="https://twitter.com/hashtag/%E6%95%B0%E6%A5%BD?src=hash&amp;ref_src=twsrc%5Etfw">#数楽</a> 続き。添付画像のように、<br><br>(-N/2～N/2-1の整数の組(a,b)で最大公約数が1のものの個数)/N^2<br><br>は6/π^2=1/ζ(2)に近くなります。<br><br>二つの整数(a,b)の最大公約数が1になる「確率」は6/π^2なのだ。 <a href="https://t.co/6GevEmVccU">pic.twitter.com/6GevEmVccU</a></p>&mdash; 黒木玄 Gen Kuroki (@genkuroki) <a href="https://twitter.com/genkuroki/status/955695607628251137?ref_src=twsrc%5Etfw">2018年1月23日</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>


6/π^2 はどれぐらいの値かといいますと大体0.6079です．つまり，大体60％になります．円周率 $\pi$ が出てきて不思議ですよね．

# 6/π^2 はどこから来るんだ！
これは 1以上の整数の逆数の二乗の無限和

```math 
\sum_{k=1}^\infty \frac{1}{k^2} = 1+ \frac{1}{2^2} + \frac{1}{3^2} + \frac{1}{4^2} + \frac{1}{5^2} + \dots
```

の値が $\frac{\pi^2}{6}$ になるというところからきています．確率はその逆数というわけです．

簡単な説明として書こうとしたら
[高校数学の美しい物語：互いに素の意味と関連する三つの定理](https://mathtrain.jp/tagainiso)
とまったく同じ事を書こうとしていたみたいなので引用するにとどめます．

式の導入において，無限の掛け算から $\sum_{k=1}^\infty \frac{1}{k^2}$ ががでてきます．これは
リーマンゼータ関数 $\zeta(s)$ の $s=2$ でのEuler積表示

```math
\sum_{k=1}^\infty \frac{1}{k^2} = \prod_{p} \frac{1}{1-p^{-2}}
```

から由来します．証明は [リーマンゼータ関数](http://integers.hatenablog.com/entry/2015/12/29/011132) が参考になります．

ほかにも別証明が
[Probability of Two Integers Being Coprime]
(https://www.cut-the-knot.org/m/Probability/TwoCoprime.shtml)
に書いてあります．数学の英語はプログラミングの英語に比べると読みやすいと思うのでチャレンジしてみてください．

# 試してみよう．（Julia版）

厳密な証明が理解できた，できないはともかく，エンジニアの皆さんはとりあえず本当にそうなのかを試したくなる，実感したくなると思います．

Gen KurokiさんのJuliaの実装がこちらにあります：

<blockquote class="twitter-tweet" data-lang="ja"><p lang="ja" dir="ltr"><a href="https://twitter.com/hashtag/Julia%E8%A8%80%E8%AA%9E?src=hash&amp;ref_src=twsrc%5Etfw">#Julia言語</a><br><br>-N/2からN/2-1までの整数の対(a,b)で最大公約数が1であるものの個数cを単純にカウントする。そして、sqrt(6N^2/c) を計算する。円周率に近い値が得られる。<br><br>私のパソコンではN=20000で2.8秒かかりました。4億回のループ。<a href="https://t.co/2NgE06gJu2">https://t.co/2NgE06gJu2</a> <a href="https://t.co/f9bnHL8s5V">pic.twitter.com/f9bnHL8s5V</a></p>&mdash; 黒木玄 Gen Kuroki (@genkuroki) <a href="https://twitter.com/genkuroki/status/955704369151819778?ref_src=twsrc%5Etfw">2018年1月23日</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

ノートブックはこちらから [互いに素な整数の対の数え上げによる円周率の計算](http://nbviewer.jupyter.org/gist/genkuroki/0c39688271a23600a3bf3ab56ae8d27b) 

ノートブックには並列で実行することで高速化するコードもあります．いずれにせよ本質的なコードはおおまかに次の通りです（ループは１からNに変更してあります）．

```julia
function calc_pi_by_gcd(N)
    s = 0
    for a in 1:N
        for b in 1:N
            s += ifelse(gcd(a,b)==1, 1, 0)
        end
        s
    end
    println(sqrt(6N^2/s))
end
```

計算結果は $\pi$ になるように調節しています．実行結果が円周率の値に近づければ正しそうなことは実感できますね．
ちなみにN=20000の場合は3.1415283804654663が出力結果で実行速度は20秒ほどでした．
 
# 試してみよう．（Python版）
また懲りずにPythonの実装をトライしています．PureなPythonは遅いのでNumbaでJITコンパイルできるようにしてあります．最大公約数をもとめるAPIは `math.gcd` で提供されていますが，Numbaによる速度最適化を優先させるために自前の実装をもちいています．

Referenceは下記の通り．
-    https://qiita.com/wrist/items/889696a507fbc8b392e4
-    http://swdrsker.hatenablog.com/entry/2017/12/20/235201
-    https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.remainder.html

```python3:gcdprob.py

import math
from numba import jit, prange
import time
import numpy as np


@jit('i8(i8,i8)', nopython=True)
def gcd(a, b):
    """
    Reference:
    https://qiita.com/wrist/items/889696a507fbc8b392e4
    http://swdrsker.hatenablog.com/entry/2017/12/20/235201
    https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.remainder.html
    """
    while b != 0:
        #a, b = b, a % b
        a, b = b, np.remainder(a, b)
    return a


@jit('f8(i8)', nopython=True)
def calc_pi_by_gcd(N):
    s = 0
    for a in range(1, N+1):
        for b in range(1, N+1):
            if gcd(a, b) == 1:
                s += 1
    pi = math.sqrt(6*N**2/s)
    return pi


@jit('f8(i8)', nopython=True, parallel=True)
def calc_pi_by_gcd_parallel(N):
    s = 0
    for a in prange(1, N+1):
        for b in range(1, N+1):
            if gcd(a, b) == 1:
                s += 1
    pi = math.sqrt(6*N**2/s)
    return pi


def main():
    start = time.time()
    pi = calc_pi_by_gcd(20000)
    end = time.time()
    print(pi)
    print("Elapsed time=", end-start)

    start = time.time()
    pi = calc_pi_by_gcd_parallel(20000)
    end = time.time()
    print(pi)
    print("Elapsed time=", end-start)

if __name__ == '__main__':
    main()
```

今回の場合prangeを使えるので並列実行もできます．
ちなみに実行結果は下記のとおりです：

```shell-session
> python gcdprob.py
3.1415283804654663
Elapsed time= 41.65773940086365
3.1415283804654663
Elapsed time= 1.2866055965423584
```

Julia に比べると2倍遅いですね・・・．

# いや，むしろJuliaの組み込みのgcdが速いのか・・・．

<blockquote class="twitter-tweet" data-lang="ja"><p lang="ja" dir="ltr"><a href="https://twitter.com/hashtag/Julia%E8%A8%80%E8%AA%9E?src=hash&amp;ref_src=twsrc%5Etfw">#Julia言語</a> の gcd の定義は<a href="https://t.co/4pCXCZIUGu">https://t.co/4pCXCZIUGu</a><br><br>にあります。Juliaで書かれています。これをNumbaに移植すれば速くなる可能性があります。trailing_zeros(n)=(nを割り切る最大の2のべきの指数)=ord_2(n)です。<br><br>私も知らないアルゴリズムでした。コードをぱっと見ても理解不能。解説が必要。</p>&mdash; 黒木玄 Gen Kuroki (@genkuroki) <a href="https://twitter.com/genkuroki/status/956170553563348994?ref_src=twsrc%5Etfw">2018年1月24日</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>


https://github.com/JuliaLang/julia/blob/d386e40c17d43b79fc89d3e579fc04547241787c/base/intfuncs.jl#L28-L49 を抜粋：

```julia
# binary GCD (aka Stein's) algorithm
# about 1.7x (2.1x) faster for random Int64s (Int128s)
function gcd(a::T, b::T) where T<:Union{Int64,UInt64,Int128,UInt128}
    a == 0 && return abs(b)
    b == 0 && return abs(a)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = unsigned(abs(a >> za))
    v = unsigned(abs(b >> zb))
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(T) && throw(OverflowError())
    r % T
end
```

へーこういうのがあるのね・・・．なにやってるのかｒｙ・・・．

# 試してみよう(その他の言語）

以下は読者の演習問題とします．
数学的にも面白いですし，アルゴリズムのお勉強としてもよい題材になると思います．
ぜひトライしてみてください．そして速度でJuliaに勝ってくだしゃい．今の私じゃ真正面で戦っても無理です．

## Nim 

暇だったので書きました．Nimの　`math.gcd` はナイーブな実装を採用しているようですね．上で書いたPython＋Numbaの実装よりスピードに負けてしまいました．Nimがんば．．．

```nim:probcoprime.nim
from math import gcd, sqrt
import times

proc calc_pi_by_gcd(N:int):float64=
  var s=0
  for a in 1..N:
    for b in 1..N:
      if gcd(a,b)==1:
        inc(s)
  var pi = sqrt(6*N*N/s)
  return pi

if isMainModule:
  var startTime=cpuTime()
  echo calc_pi_by_gcd(1000)
  var endTime=cpuTime()
  echo endTime-startTime
```

## Nim Stein版
この実装でPython+Numbaの実装に追いつきました．ふぅ・・・．

```nim
from math import sqrt 
import times 

proc trailing_zeros(n:int): int =
  var order = 0
  var n = n
  while (n and 1) == 0:
    inc(order)
    n = n shr 1
  return order

proc stein_gcd(a, b:int):int=
  if a == 0:
    return abs(b)
  if b == 0:
    return abs(a)
  let za = trailing_zeros(a)
  let zb = trailing_zeros(b)
  let k = min(za, zb)

  var u = abs(a shr za)
  var v = abs(b shr zb)

  while u != v:
    if u > v:
        swap u,v
    v -= u
    v = v shr trailing_zeros(v)

  let r = u shl k
  return r

proc calc_pi_by_gcd(N:int): float64 = 
  var s = 0
  for a in 1..N:
    for b in 1..N:
      if stein_gcd(a,b)==1:
        inc(s)
  var pi=sqrt(6*N*N/s)
  return pi 

if isMainModule:
  let N = 20000
  var s=cpuTime()
  var pi=calc_pi_by_gcd(N)
  var e=cpuTime()
  echo "pi=",pi
  echo e-s,"[sec]"
```

## Rust 

とりあえず，実装．

```rust
extern crate num;
use num::Integer;
use std::f64;
use std::time::Instant;

fn calc_pi_by_gcd(n:i64)->f64 {
    let mut s=0;
    for a in 1..n+1 {
        for b in 1..n+1{
            if a.gcd(&b)==1{
                s+=1;
            }
        }
    }

    let pi = ((6*n*n) as f64/s as f64).sqrt();
    return pi;

}


fn main(){
    let n = 20000;
    let start = Instant::now();
    let pi = calc_pi_by_gcd(n);
    let end = start.elapsed();
    println!("Pi = {}",pi);
    println!(
        "Elapsed time = {}.{:03}",
        end.as_secs(),
        end.subsec_nanos() / 1_000_000
    );
}
```

```toml:Cargo.toml
[package]
name = "rustgcd"
version = "0.1.0"

[dependencies]
num="0.1.41"

[profile.release]
opt-level = 3
debug = false
rpath = false
lto = true
debug-assertions = false
codegen-units = 1
panic = 'unwind'
```

```shell-session
> cargo run --release
Pi = 3.1415283804654663
Elapsed time = 17.015
```

おおお・・・．20秒のJuliaに勝てました．