multiplicable

{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}


import Text.Printf (printf)
import System.CPUTime (getCPUTime)

class Multiplicable a where
    one  :: a
    (×)  :: a -> a -> a

{-# INLINABLE pow #-}
{-# SPECIALIZE pow :: Integer -> Integer -> Integer #-}
pow :: forall a. (Multiplicable a) => a -> Integer ->  a
pow _ 0 = one
pow m 1 = m
pow m n = let
    k = m × m
    f :: Integer -> [a] -> (a, Integer) -> a
    f _ _ (v,0) = v
    f n (e:es) (v,r)
        | r - u >= 0 = f u (w:e:es) (v × w, r - u) -- double
        | r - n >= 0 = f n (e:es) (v × e, r - n) -- stay
        | otherwise = f (n `div` 2) es (v, r) -- half
        where u = n * 2
              w = e × e
    in f 2 [k, m] (k, n - 2)

-- Tests
time :: a -> IO ()
time a = do
    start <- getCPUTime
    let !r = a
    end <- getCPUTime
    let diff = fromIntegral (end - start) / (10^12)
    printf "Computation time: %0.3f sec\n" (diff :: Double)
    return ()

instance {-# OVERLAPPABLE #-}(Num a) => Multiplicable a where
    one  = 1
    (×)  = (*)

-- test with matrixes
data Matrix a where
    Matrix :: (Num a) => (a, a, a, a) -> Matrix a

instance {-# INCOHERENT #-}(Num a) => Multiplicable (Matrix a) where
    one  = Matrix (1, 0, 0, 1)
    (×) (Matrix (a1, a2, a3, a4)) (Matrix (b1, b2, b3, b4)) = let
        r1 = a1 * b1 + a2 * b3
        r2 = a1 * b2 + a2 * b4
        r3 = a3 * b1 + a4 * b3
        r4 = a3 * b2 + a4 * b4
        in Matrix (r1, r2, r3, r4)

-- matrix powered fibonacci solver
fib :: Integer -> Integer
fib 0 = 0
fib 1 = 1
fib n = r1 + r2 where
    (Matrix (r1, r2, _, _)) = pow (Matrix (1,1,1,0)) (n - 2)

-- clasical tail recursive fibonacci solver
fib' :: Integer -> Integer
fib' = f 0 1  where
    f a _ 0 = a
    f a b n = f b (a + b) (n - 1)

main :: IO()
main = do
    putStr "fib 10^6\t"
    time (fib $ 10^6)
    putStr "fib' 10^6\t"
    time (fib' $ 10^6)
    putStr "pow 2 . (pow 10) $ 9\t"
    time (pow 2  . pow 10  $ 9)
    putStr "2^10^9\t"
    time (2^10^9)
    putStrLn ""