Fib.hs

{-# LANGUAGE NumericUnderscores #-}

import Data.Int

modulus = 1_000_000_007

-- M = Z / 1_000_000_007 Z
newtype M = M Int64 deriving (Eq, Ord)

-- Ext = M[x]/(x^2 - 5)
-- Note: x^2 - 5 is irreducible over M, so this is a field
data Ext = Ext M M deriving (Eq)

-- Functions used for printing M and Ext
instance Show M where
  show (M n) = show n

instance Show Ext where
  show (Ext a 0) = show a
  show (Ext a b) = (shows a . showString " + " . shows b . showString " * root5") ""

-- Implementations of arithmetical operations
instance Num M where
  M a + M b = M $ (a + b) `mod` modulus
  M a * M b = M $ (a * b) `mod` modulus

  negate (M a) = M $ negate a `mod` modulus
  fromInteger i = M $ fromInteger i `mod` modulus

  abs = id
  signum = error "unimplemented"

instance Num Ext where
  Ext a b + Ext c d = Ext (a + c) (b + d)
  Ext a b * Ext c d = Ext (a * c + b * d * 5) (a * d + c * b)

  negate (Ext a b) = Ext (negate a) (negate b)
  fromInteger i = Ext (fromInteger i) 0

  abs = id
  signum = error "undefined"

-- Modular inverse
-- Note that a ^ b does exponentiation by squaring
inverse :: M -> M
inverse x = x ^ (modulus - 2)

-- Conjugation
conjugate :: Ext -> Ext
conjugate (Ext a b) = Ext a (negate b)

-- Some constants
root5, root5inv, inv2, goldenRatio, goldenRatioConj :: Ext
root5 = Ext 0 1
root5inv = Ext 0 (inverse 5)
inv2 = Ext (inverse 2) 0
goldenRatio = (1 + root5) * inv2
goldenRatioConj = conjugate goldenRatio

-- Fibonacci!!!
fib :: Int -> M
fib n = case (goldenRatio ^ n - goldenRatioConj ^ n) * root5inv of
  Ext a 0 -> a
  _ -> error "this won't happen"