天泣記

ふと、フィボナッチ数を求めるプログラムを4つ書いてみた。素朴な再帰版、ループ版、行列版、浮動小数点版の4つで、おおざっぱには最初の方が遅いはずである。行列版、浮動小数点版はフィボナッチ数の一般項を利用したものである。

require 'matrix'

def fib_rec(n)
  return n if n <= 1
  fib_rec(n-2) + fib_rec(n-1)
end

def fib_iter(n)
  return n if n <= 1
  f0, f1 = 0, 1
  (n-1).times { f0, f1 = f1, f0 + f1 }
  f1
end

def fib_mat(n)
  return n if n <= 1
  (Matrix[[1,1],[1,0]] ** (n-1))[0,0]
end

Root5 = Math.sqrt(5)
C1 = (1 + Root5) / 2
C2 = (1 - Root5) / 2
def fib_flt(n)
  v = (C1**n - C2**n) / Root5
  v.finite? ? v.round : nil
end

とりあえず小さい値について結果を表示してみると、期待通りどれも一致する。

0.upto(10) {|n|
  p %i[fib_rec fib_iter fib_mat fib_flt].map {|m|
    self.send(m, n)
  }
}
#=>
[0, 0, 0, 0]
[1, 1, 1, 1]
[1, 1, 1, 1]
[2, 2, 2, 2]
[3, 3, 3, 3]
[5, 5, 5, 5]
[8, 8, 8, 8]
[13, 13, 13, 13]
[21, 21, 21, 21]
[34, 34, 34, 34]
[55, 55, 55, 55]

fib_flt は Float を使っていて誤差があるので、そのうち失敗する(正しくない値を返す)わけだが、どこで失敗するか調べてみると、fib(71) で失敗するようだ。そのときの結果の値は49ビットであり、IEEE754倍精度なら(仮数部は53ビットあるので)まだ表現できるはずなので、うまく計算すればもうちょっといけるかもしれない。

0.upto(100) {|n|
  f1 = fib_iter(n)
  f2 = fib_flt(n)
  if f1 != f2
    p [n, f1, f2, f1.bit_length]
    exit
  end
}
#=>
[71, 308061521170129, 308061521170130, 49]

さらに大きくなると fib_flt は内部的に Infinity になってしまって、nil を返すはずである。調べてみると、fib(1475) でそうなるようだ。

0.upto(2000) {|n|
  f = fib_flt(n)
  if f.nil?
    p n
    exit
  end
}
#=>
1475

実行時間を測ってみよう。

require 'csv'
puts 'alg,n,t[s]'
%i[fib_rec fib_iter fib_mat fib_flt].each {|m|
  0.upto(10000) {|n|
    t1 = Process.clock_gettime(Process::CLOCK_THREAD_CPUTIME_ID)
    v = self.send(m, n)
    t2 = Process.clock_gettime(Process::CLOCK_THREAD_CPUTIME_ID)
    break if !v
    t = t2-t1
    puts [m,n,t].to_csv
    break if 1 < t
  }
}
#=>
alg,n,t[s]
fib_rec,0,4.8079999999950385e-06
fib_rec,1,1.582999999999446e-06
...

グラフにしてみると以下のようになる。だいたいは予想通りで、fib_rec はすぐに激しく遅くなり、fib_iter はそこそこ、fib_mat は n が大きければ速い。fib_flt はとても速いが、71以上で誤差があるし、1475以上は概数さえ求められない。

fib_iter より fib_mat が速くなるのはかなり n が大きくなってからで、行列を扱うのはけっこうオーバーヘッドがあるようだ。

あと、なんか fib_iter と fib_mat が二股になっているが、よくわからない。真面目に調べていないが、とりあえず GC のせいかもしれないといっておこう。(可能性はあると思うが、確信はない)

fib.R:

library(ggplot2)
d <- read.csv("2016-03/fib.csv")
p <- ggplot(d, aes(n, t.s., shape=alg, color=alg)) +
  geom_point() +
  scale_x_log10() +
  scale_y_log10("t[s]")
print(p)

天泣記

2016-03-26 (Sat)

2016-03-29 (Tue)