Hmatrix and `Integers`

#1

Inspired by this post, I wrote a fibonacci implementation using hmatrix

{-# language OverloadedLists #-}
import Numeric.LinearAlgebra
import Data.Semigroup (stimes)

-- can't use stimesMonoid to eliminate 0 case, because matrix mempty has the wrong dimension (1><1)
fib :: Int -> Int
fib 0 = 0
fib n = fromIntegral x
  where
    [x, _] = stimes n fibMatrix #> [0, 1]

    fibMatrix :: Matrix Z
    fibMatrix = (2><2) [1, 1, 1, 0]

main = do
  mapM_ print $ takeWhile (>=0) $ map fib [0..]

Unfortunately, it doesn’t get very far into the sequence before overflowing the size of an Int64 (that’s what Z is an alias for). Unfortunately I can’t use Integer - it seems that hmatrix only supports Storable elements. I guess this is because hmatrix is a wrapper over fast c libraries.

So, my questions:

  • Long shot, but is there a way to get around this issue without switching libraries? I like hmatrix's API a lot.
  • If not, what library would you recommend that can handle arbitrary precision integers?
#2

Want to give laop a try? Doubt it will be as fast but maybe you might be able to write this function using a different perspective on linear algebra and statically typed matrices :grinning:

1 Like
#3

The linear package also allows you to use arbitrary Num types such as Integer.

2 Likes
#4

Also note you do not need a 2x2 matrix. You can do it with just 2 Integers:

import Data.Semigroup

data Fib = Fib !Integer !Integer

instance Semigroup Fib where
  Fib x1 y1 <> Fib x2 y2 = Fib (x2 * x1 + y2 * y1) (y2 * x1 + (x2 + y2) * y1)

fib :: Integer -> Integer
fib 0 = 0
fib n = (\(Fib _ x) -> x) (stimes n (Fib 0 1))
4 Likes
#5

I’ll take a look!

Unfortunately you don’t get semigroup/monoid instances out of the box so you’d have to write them yourself, for a newtype wrapper over M22. Or just write your own exponentiation by squaring function instead of using stimes I guess.

That’s interesting!

1 Like