前回、無理数の取り扱いを複素数的な考えで取り扱ってフィボナッチ数列の一般項を求めた。意味的にはさほど遅くないはずの実装ではあるが、それでも 100,000,000 番目を求めようとすると結構時間がかかってしまう。実際に計測するとこんな感じだ。
90.67s user 0.21s system 99% cpu 1:30.94 totalチューニングの余地があるのは恐らく素朴に書いた累乗計算、fibExp だ。素朴と言っても一応アルゴリズム的な考慮はしたつもりだが、とはいえ言語について詳しく知っているわけでもなし、色々素人が小手先のチューニングを施すよりも、標準ライブラリに任せたほうがこういうものは余程パフォーマンスがいいものと相場が決まっている。
Haskell での累乗計算は (^) を使う。GHCi に聞いてみると、
> :i (^)
(^) :: (Num a, Integral b) => a -> b -> a -- Defined in `GHC.Real'
infixr 8 ^ということなので、a が Num クラスに属してさえいれば、a はなんだっていい。ということは、前回の FibNum を Num クラスのインスタンスにしてしまえばいいわけだ。それさえできれば標準の (^) が独自の型にも適用できるようになる。素晴らしい。
標準ライブラリを活用するために Num クラスのインスタンスにする
さて、前回は手抜き(?)でタプルを利用したが、タプルのまま型クラスのインスタンスにしようとすると下記のように叱られてしまう。
Illegal instance declaration for `Num FibNum'
(All instance types must be of the form (T t1 ... tn)
where T is not a synonym.というわけできちんと型を定義しよう。今考えている a + b√5 は黄金比にちなんだ数なので、名前も付け直して GoldenExp とでもしようか。…と思ったがモノに似合わず壮大な名前だなぁ。
「ふん GoldenExp というのかい。ぜいたくな名だね。今からお前の名前は GExp だ。いいかい GExp だよ。わかったら返事をするんだ、GExp！」
import Data.Ratio
data GExp = Rational :+: Rationalというわけでこんな感じで。:+: という表記は見慣れないかもしれないが、値コンストラクタの中置記法だ。こんな感じ：
> :t (:+:)
(:+:) :: Rational -> Rational -> GExp要するに (:) とかと同じものを自前で定義できるというわけ。ただし : で始まる記号列に限るという制限があるが。もちろん data GExp = GExp Rational Rational みたいにして書いてもいい。
> :t GExp
GExp :: Rational -> Rational -> GExp同じことだ。単に今回は a + b√5 を表すデータ型なので、a :+: b と書けたほうが少しは見やすかろうという程度である。
さて Num 型クラスはこうなっている。
> :i Num
class Num a where
(+) :: a -> a -> a
(*) :: a -> a -> a
(-) :: a -> a -> a
negate :: a -> a
abs :: a -> a
signum :: a -> a
fromInteger :: Integer -> a加算、乗算、減算、符号反転、絶対値、符号、整数からの変換、の 7 種類の定義を提供せよ、と言っている。「とにかくこういう操作が出来るヤツは Num である」というわけだ。もちろん満たすべき性質はあるのだが、ここでは数学の常識をそのまま当てはめれば問題ない。
乗算、減算は既に前回定義している。加算も符号反転も整数からの変換も難しくない。さっさと書いてしまおう。
instance Num GExp where
(a :+: b) + (c :+: d) = (a + c) :+: (b + d)
(a :+: b) - (c :+: d) = (a - c) :+: (b - d)
(a :+: b) * (c :+: d) = (a*c + 5*b*d) :+: (a*d + b*c)
negate (a :+: b) = negate a :+: negate b
fromInteger a = (a % 1) :+: 0絶対値と符号は？ 今回は使わないので別に未定義でもいいが：
signum = undefined
abs = undefinedせっかくなので書こう。まず符号調査。GExp は a + b√5 を表現しているのだから、a b 共に符号が同じであればその符号で良い。問題は違っている場合で、-2 + √5 は正なの？ 負なの？ と言われると迷う。が、よくよく考えると、全体の符号は a と b√5 で絶対値が大きい方の符号に従うことが分かる。√5 は無理数だが二乗すれば整数だ。絶対値の大小比較は二乗の大小比較に等しい。ということは：
signum (a :+: b)
| a*a < 5*b*b = signum b :+: 0
| otherwise = signum a :+: 0こうか。
これができれば絶対値は難しくない。
abs a = signum a * aこうだね。ちなみに Haskell では関数適用はどんな演算子よりも優先順位が高い。だから (signum a) * a と読む。
計算してみる
これで GExp は Num クラスのインスタンスとなった。嬉しいことに (^) が適用できる。累乗計算を書かなくて済む上に、表記もかなり自然になる。
fibMain :: Integer -> GExp
fibMain n = (((1 % 2) :+: (1 % 2))^n - ((1 % 2) :+: (-1 % 2))^n) * (0 :+: (1 % 5))
fib :: Integer -> Integer
fib n = numerator a
where (a :+: _) = fibMain n素晴らしい。
実行時間を見てみよう。
50.13s user 0.19s system 99% cpu 50.362 totalお、さすが標準ライブラリ。大分速くなった。きっと素人が書くよりも効率の良い実装が行われているのだろう。
最後の除算をやめる
ライブラリ活用とは逸れるが、ついでに。
最後の除算――今回は乗算の形にしてある部分、つまり分母の √5 のことだが、この直前の形は (a + b√5)/√5 となっている。そして、ここで a は絶対に 0 のはずである。なぜなら a が 0 以外なら、√5 で割ったときに √5 の項が残ってしまうからだ。
ということは結局 b√5/√5 ということで、この除算は単に b を左の項に移動するだけの効力しかないことになる。だがそんな移動だけならわざわざ計算ををする必要はない。パターンマッチで取り出せるのだから取り出してしまえば良い。
fib :: Integer -> Integer
fib n = numerator a
where (_ :+: a) = ((1 % 2) :+: (1 % 2))^n - ((1 % 2) :+: (-1 % 2))^nたかが乗算一回減ったくらいじゃ大差もあるまいが、
46.58s user 0.32s system 99% cpu 46.925 totalまあ桁数が多いと結構バカにならぬらしい。3 秒くらいは短縮された。
それにしても元々 90 秒ほどかかっていたことを考えると約倍速であるし、実に標準ライブラリ様様である。
完成品
signum やら abs やら一部不要なものも定義したが、自分で累乗を書いていた前回のコードよりもすっきりした感じがする。ま、main は相変わらず不親切設計だが。
import System.Environment (getArgs)
import Data.Ratio
data GExp = Rational :+: Rational
instance Num GExp where
(a :+: b) + (c :+: d) = (a + c) :+: (b + d)
(a :+: b) - (c :+: d) = (a - c) :+: (b - d)
(a :+: b) * (c :+: d) = (a*c + 5*b*d) :+: (a*d + b*c)
negate (a :+: b) = negate a :+: negate b
fromInteger a = (a % 1) :+: 0
signum (a :+: b)
| a*a < 5*b*b = signum b :+: 0
| otherwise = signum a :+: 0
abs a = signum a * a
fib :: Integer -> Integer
fib n = numerator a
where (_ :+: a) = ((1 % 2) :+: (1 % 2))^n - ((1 % 2) :+: (-1 % 2))^n
main :: IO ()
main = mapM_ (print . fib . read) =<< getArgs
＃人気の記事
＃タグ
