本日は
先日 速度でJuliaに勝つのが難しかった（というか無理な）話．
という記事を書きましたが，多くの方に読んていただいたみたいです．ありがとうございます．
Twitter で Gen Kuroki さんがJuliaネタ，数学ネタを投下していましたので，このネタでも単純ですが速度で挑戦してみましょう．
元ネタはこちらから
#数楽 続き。添付画像のように、
— 黒木玄 Gen Kuroki (@genkuroki) 2018年1月23日
(-N/2～N/2-1の整数の組(a,b)で最大公約数が1のものの個数)/N^2
は6/π^2=1/ζ(2)に近くなります。
二つの整数(a,b)の最大公約数が1になる「確率」は6/π^2なのだ。 pic.twitter.com/6GevEmVccU
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の実装がこちらにあります：
#Julia言語
— 黒木玄 Gen Kuroki (@genkuroki) 2018年1月23日
-N/2からN/2-1までの整数の対(a,b)で最大公約数が1であるものの個数cを単純にカウントする。そして、sqrt(6N^2/c) を計算する。円周率に近い値が得られる。
私のパソコンではN=20000で2.8秒かかりました。4億回のループ。https://t.co/2NgE06gJu2 pic.twitter.com/f9bnHL8s5V
ノートブックはこちらから 互いに素な整数の対の数え上げによる円周率の計算
ノートブックには並列で実行することで高速化するコードもあります．いずれにせよ本質的なコードはおおまかに次の通りです（ループは１から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
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が速いのか・・・．
#Julia言語 の gcd の定義はhttps://t.co/4pCXCZIUGu
— 黒木玄 Gen Kuroki (@genkuroki) 2018年1月24日
にあります。Juliaで書かれています。これをNumbaに移植すれば速くなる可能性があります。trailing_zeros(n)=(nを割り切る最大の2のべきの指数)=ord_2(n)です。
私も知らないアルゴリズムでした。コードをぱっと見ても理解不能。解説が必要。
# 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がんば．．．
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
);
}
[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に勝てました．