Closed form Fibonacci using Integer

Question

I'm implementing the closed form of Fibonacci numbers. I figure using a custom data type I can do so without leaving Integer. It isn't more efficient than the "standard" approach as I couldn't figure an efficient way to calculate the power of my custom number type.

But sticking to this algorithm I would be interested in best practices to write this. In particular I'm uncertain about the use of Num for my FibNum and the use of fromJust.

module Fibo
    ( 
     fibo
    ) where

import Data.Maybe (fromJust)

-- use a data type similar to complex numbers. 
-- The only irrational number are multiples of sqrt 5 
-- which will cancel away in the closed form.
data FibNum = FibNum { int::Integer, sq5::Integer }
  deriving (Show)

instance Num FibNum where
  (+) (FibNum a1 b1) (FibNum a2 b2) = FibNum (a1+a2) (b1 + b2)
  (*) (FibNum a1 b1) (FibNum a2 b2) = FibNum intSide rootSide
    where intSide = a1*a2 + 5 * b1 * b2
          rootSide = a1*b2 + b1*a2
  negate (FibNum a b) = FibNum (negate a) (negate b)
  fromInteger n = FibNum (fromInteger n) 0
  -- no need for abs and sign for now
  abs _ = undefined
  signum _ = undefined

funit = FibNum 1 0
-- actual phi and psi are (1 + sqrt 5)/2 and (1 - sqrt 5)/2
-- I would prefer to avoid fractions so I divide at the end
fphi = FibNum 1 1
fpsi = FibNum 1 (negate 1)

-- was originally planning to calculate the binomial coefficients
-- but seeing that I just need n relatively simple multiplications 
-- that seem like unnecessary complication
pow :: FibNum -> Integer -> FibNum
pow _ 0 = funit
pow x 1 = x
pow x n = x * pow x (n-1)

-- in the closed form the integer part is always 0
divR5 :: FibNum -> Maybe Integer
divR5 (FibNum 0 n) = Just n
divR5 _ = Nothing

fibo :: Integer -> Integer
fibo n = flip div p2 . fromJust $ divR5 nom
  where nom = pow fphi n - pow fpsi n
        p2 = 2^n

addendum Thanks a lot for the useful comments. If it's interesting I put my code in following repo https://github.com/bdcaf/haskell-fibonacci

Answer 1

Opening comments

I'm not a Haskell expert, but I hope I can offer some helpful comments and critiques.

Overall, I think this is mostly good and idiomatic code. You make good use of pattern matching and have pretty clean definitions of your code.

I'm going to first address style and your questions. I'll then address the performance issue, since it actually isn't that hard to fix.

Style

I think most of your style is fine. I just personally prefer to nest where clauses like so

foo x y = x' + y' 
  where
    x' = 2*x
    y' = 2*y

How you do it is up to you.

Using Num

I don't think there's much wrong with using Num for FibNum to make things easier. I would just go ahead and implement abs and signum if you do, since they aren't too difficult, although you don't have to since you don't export FibNum (which I think is good). However, Num is just an interface and doesn't officially make any promises about how its operations behave (i.e. it doesn't have any laws).

The data type you've defined is the Ring* Z[√5]. I found a library that has a Ring type class which you could implement for your datatype. You could then use this type class's implementation of (^) instead of pow for a speedup (see my comments regarding performance).

Whether you want to be more precise like this is up to you. Num unofficially is a ring, too, and is the more common type class.

As an addendum, since you are using division, you may wish to consider changing your Integers to Rationals, which would make your type the Quadratic Field Q[√5]. There is also a Field typeclass.

* Don't worry too much about the abstract algebra here if you aren't familiar. What's more important is understanding the interfaces and the promises they make/laws they obey. For example, the order of operation of addition can be rearranged (known as commutativity) in a Ring and Field. You don't need to know anything more than how to implement the functions required, make sure they satisfy the laws, and how to use the interface (I don't really remember my rings and fields anyway, I just googled around a bit).

Using fromJust

I think that the usage of fromJust is unnecessary and serves to undermine the fact that divR5 returns a Maybe. Below, I offer two alternative options.

I would also consider either renaming divR5 to indicate that it doesn't actually divide an arbitrary FibNum by √5, leaving a comment to that effect, or defining it in the where clause of fibo (which is what I would do). If you define it in the where clause, you can't accidentally misuse it elsewhere.

Return a Maybe

fiboMaybe :: Integer -> Maybe Integer
fiboMaybe n = (`div` p2) <$> divR5 nom
  where
    nom = pow fphi n - pow fpsi n
    p2  = 2^n

I like to use the infix form of fmap (and also I prefer to use div in infix), but you can also explicitly case on the result or use do notation.

Give a better error message

If you aren't going to return a Maybe, then you might as well provide a more informative error message than fromJust's.

fiboError :: Integer -> Integer
fiboError n = divR5 nom `div` p2
  where
    nom = pow fphi n - pow fpsi n
    p2  = 2^n
    divR5 (FibNum 0 n) = n
    divR5 _            = error "fiboError: got nonzero integer part in divR5 (this shouldn't happen)"

Improving efficiency

It turns out that it's not too hard to improve your efficiency. If you implement the Ring or Field type class and use the functions provided there, you should see a major speedup. But if you don't want to, we can fix the speed in only a couple lines of code.

The inefficiency in your code comes from your implementation of pow, which, while correct, takes linear time. You can reduce this significantly. Here's how we can fix it. We're going to use the Product Monoid and the function mtimesDefault from Data.Semigroup.

If you aren't familiar with the abstract algebra terminology here, ignore that junk for a second. Here's the lowdown: we're going to take advantage of the fact that for your datatype, multiplication is associative. What does that mean? It means that

(x * y) * z == x * (y * z)

i.e. we can move around parentheses in a product without changing its value. If you can do that, you have something known as a Semigroup. That's all a Semigroup is! A Monoid is a Semigroup where you know there's some element that does nothing when you multiply it. In this case, that element is 1:

1 * x == x * 1 == x

If you have a Monoid (which FibNum is with respect to the multiplication operation) and want to multiply a number by its self n times, mtimesDefault :: (Integral b, Monoid a) => b -> a -> a does this more efficiently than the naive solution.

The Product wrapper type takes a Num a and uses its multiplication operation as the Monoid operation. So to get a faster pow, here's all the code we have to write:

import Data.Monoid    (Product(..))
import Data.Semigroup (mtimesDefault)

powFast :: Integral a => FibNum -> a -> FibNum
powFast n exp = getProduct $ mtimesDefault exp (Product n)

If we replace pow with powFast, your implementation becomes much faster than the "standard" approach too!

Of course, it isn't too hard to write a faster pow function by hand. It's just neat that Haskell has built-in machinery that lets you avoid doing so. If you wanted to figure out how to do it faster by hand, I would hint you to think about exponentiation of regular numbers.

Say I asked you to compute 2^50 by hand. I claim you don't need to take 50 multiplications to give me an answer. Try and think about how you would do it efficiently, taking advantage of the fact that you're only ever multiplying by 2 and that you multiplication is associative.

Here's a hint:

Think about multiplying exponents with the same base. 2^x * 2^y = 2^(x+y). What about 2^x * 2^x?

and another:

For example, suppose you get to 2^4 = 16. From here, if you multiply 16 by itself, you'll get 16 * 16 = (2^4) * (2^4) = 2^(2*4) = 2^8. That saved 3 multiplications compared to the naive method of multiplying by 2!

You can find the answer in the implementation of stimesDefault.