Python
数学
Julia
Nim
Numba

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

本日は

先日 速度でJuliaに勝つのが難しかった(というか無理な)話.
という記事を書きましたが,多くの方に読んていただいたみたいです.ありがとうございます.

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

元ネタはこちらから

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

6/π^2 はどこから来るんだ!

これは 1以上の整数の逆数の二乗の無限和

\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}$ になるというところからきています.確率はその逆数というわけです.

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

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

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

から由来します.証明は リーマンゼータ関数 が参考になります.

ほかにも別証明が
Probability of Two Integers Being Coprime
に書いてあります.数学の英語はプログラミングの英語に比べると読みやすいと思うのでチャレンジしてみてください.

試してみよう.(Julia版)

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

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

ノートブックはこちらから 互いに素な整数の対の数え上げによる円周率の計算

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

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

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を使えるので並列実行もできます.
ちなみに実行結果は下記のとおりです:

> python gcdprob.py
3.1415283804654663
Elapsed time= 41.65773940086365
3.1415283804654663
Elapsed time= 1.2866055965423584

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

いや,むしろJuliaの組み込みのgcdが速いのか・・・.

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

# 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

へーこういうのがあるのね・・・.なにやってるのかry・・・.

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

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

Nim

暇だったので書きました.Nimの math.gcd はナイーブな実装を採用しているようですね.上で書いたPython+Numbaの実装よりスピードに負けてしまいました.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の実装に追いつきました.ふぅ・・・.

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

とりあえず,実装.

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
    );
}
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'
> cargo run --release
Pi = 3.1415283804654663
Elapsed time = 17.015

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