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を指数回かけ算してるんじゃないのかね。うーむ...
