2の1000乗

Scheme

28日にあった第二回YMPAミーティングで、2の1000乗の計算について話しが出た。長整数でもオーバーフローしてしまうので求められなかったという声が聞こえたので、Gaucheでやってみようと思った。

(define pow
  (lambda (x n)
    (if (= n 0)
        1
        (* x (pow x (- n 1))))))

gosh> (pow 2 1000)
gosh> 10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376

ふむ、あっさりできた。

しかし、Schemer見習いとしては、末尾再帰にしたいところ。末尾再帰にすることでどんだけ早くなるのかも気になる。っちゅーことで、やってみた。手続き名の末尾の "/tr" は "with tail recursion" のつもり。

(define pow/tr
  (lambda (x n)
    (let loop ((m n) (p 1))
      (if (= m 0)
          p
          (loop (- m 1) (* p x))))))

gosh> (time (pow/tr 2 1000))
;(time (pow/tr 2 1000))
; real   0.000
; user   0.000
; sys    0.000
10715086071862673209484250490600018105614048117055336074437503883703510511249361224931983788156958581275946729175531468251871452856923140435984577574698574803934567774824230985421074605062371141877954182153046474983581941267398767559165543946077062914571196477686542167660429831652624386837205668069376

...一瞬すぎてわからん！1000乗では足りないようなので、10000乗してみる。あと、いちいち結果を出力されるのはウザイので、マクロtimeに引数追加してやった。

gosh> (time (pow/tr 2 10000) #t)
;(time (pow/tr 2 10000) #t)
; real   0.034
; user   0.030
; sys    0.000
#t

もうちょっと長い方がいいな。ってことで、100000乗してみる。

gosh> (time (pow/tr 2 100000) #t)
;(time (pow/tr 2 100000) #t)
; real   3.209
; user   3.110
; sys    0.050
#t

これくらいでいいか。じゃぁ、powでもやってみる。

gosh> (time (pow 2 100000) #t)
;(time (pow 2 100000) #t)
; real   2.910
; user   2.870
; sys    0.020
#t

なんですと！？非末尾再帰版の方が速いやん...。どういうこと？なんか間違ってるのかな。複数回やって平均とか取った方がいいのかな...。

話しは変わるけど

28日のYMPAで「2の1000乗の値が計算できましたよ」って結果を見せると、schemeのソースを欲しいと言われた。だけど渡す手段がなかったので、とりあえずtweetしといたら（pow/trはこれを後ほどnamed-letで書き直したもの）、SaitoAtsushiさんから返信を頂いた。

@yagiey_tw (define(pow-t x n)(if(zero? n)1(*(if(zero?(modulo n 2))1 x)(pow-t(* x x)(quotient n 2))))) 指数を2のべき乗に分解するとO(log(n))で計算できることが知られてます

http://twitter.com/SaitoAtsushi/status/22356470814

ということなので、これで時間を計測してみる。

gosh> (time (pow/tr2 2 100000) #t)
;(time (pow/tr2 2 100000) #t)
; real   0.027
; user   0.020
; sys    0.000
#t

さすがO(log(n))、桁違いに速えぇ。元のソースはO(n)なんだよな、たぶん。

あれ？

(pow/tr 2 10000)と(pow/tr 2 100000)にかかった時間から判断すると、O(n^2)かな？でも末尾再帰に書き換えてループと等価になるんだったら、単に2を指数回かけ算してるんじゃないのかね。うーむ...