Z[φ]

ZPhi.hs

module ZPhi where

-- | When @a@ is the integers, this represents the ring Z[φ] of numbers of the form a + bφ, where φ = (1 + √5)/2 is the golden ratio. @'ZPhi' a b@ corresponds to the number a + bφ.
data ZPhi a = ZPhi !a !a
         deriving (Eq,Show)

-- | @'fib' n@ computes the n'th Fibonacci number in Θ(log |n|) time. Negative arguments are permitted.
fib :: (Integral i, Num a) => i -> a
fib n = let ZPhi _ r = phiPow n in r

instance Functor ZPhi where
    fmap f (ZPhi a b) = ZPhi (f a) (f b)

instance Num a => Num (ZPhi a) where
    fromInteger n = ZPhi (fromInteger n) 0
    ZPhi a b + ZPhi c d = ZPhi (a + c) (b + d)
    ZPhi a b * ZPhi c d = ZPhi (a*c + b*d) (a*d + b*c + b*d)
    -- the rest is not really needed
    ZPhi a b - ZPhi c d = ZPhi (a - c) (b - d)
    abs = error "no abs for ZPhi"
    signum = error "no signum for ZPhi"

instance Fractional a => Fractional (ZPhi a) where
    recip a = fmap (/ norm a) (conj a)
    fromRational n = ZPhi (fromRational n) 0

instance Integral a => Ord (ZPhi a) where
    ZPhi a b <= ZPhi c d = case compare b d of
        EQ -> a <= c
        GT -> gtPhiDiv (c-a) (b-d)
        LT -> not $ gtPhiDiv (c-a) (b-d)

instance Integral a => Real (ZPhi a) where
    toRational (ZPhi a b) = toRational a + toRational b * phi'
        where phi' = toRational $ (1 + sqrt 5) / 2

-- | φ = (1 + √5)/2
phi :: Num a => ZPhi a
phi = ZPhi 0 1

-- | φ' = 1 - φ is the conjugate of φ.
phi' :: Num a => ZPhi a
phi' = ZPhi 1 (-1)

-- | 1/φ = φ - 1 is the reciprocal of φ.
rPhi :: Num a => ZPhi a
rPhi = ZPhi (-1) 1

-- | The n'th power of φ for an integer n, possibly negative.
phiPow :: (Integral i, Num a) => i -> ZPhi a
phiPow n | n < 0     = rPhi ^ (-n)
         | otherwise = phi ^ n

-- | The conjugate mapping, sending a + bφ to a + bφ'.
conj :: Num a => ZPhi a -> ZPhi a
conj (ZPhi a b) = ZPhi (a+b) (-b)

-- | The norm mapping, sending a + bφ to the integer (a + bφ)(a + bφ').
norm :: Num a => ZPhi a -> a
norm (ZPhi a b) = a^2 + a*b - b^2

toFloating :: (Real a, Floating b) => ZPhi a -> b
toFloating (ZPhi a b) = realToFrac a + (1+sqrt 5)/2 * realToFrac b

-- | 'getPhiDiv' a b is True iff a/b >= 'φ'
gtPhiDiv :: (Num a, Ord a) => a -> a -> Bool
gtPhiDiv a b = case compare b 0 of
    EQ -> error "gtPhiDiv: zero denominator"
    LT -> gtPhiDiv (-a) (-b)
    GT | a <= b    -> False
       | otherwise -> not (gtPhiDiv b (a-b))

-- | @'binarySearch' p a b@ returns the largest number in [a,b] where p is False,
-- assuming @p@ is monotonic, @p a = False@, and @p b = True@
binarySearch :: Integral a => (a -> Bool) -> a -> a -> a
binarySearch p a b | b <= a+1  = a
                   | p c       = binarySearch p a c
                   | otherwise = binarySearch p c b
    where c = a + (b - a) `div` 2

-- | Given n, computes the floor of nφ.
floorPhiTimes :: Integral a => a -> a
floorPhiTimes n = case compare n 0 of
    LT -> -floorPhiTimes (-n)
    EQ -> 0
    GT -> binarySearch ((>= ZPhi 0 n) . (`ZPhi` 0)) n (2*n)
-- I'm too lazy to figure out something faster.

-- | Floor and ceiling.
floorZPhi, ceilingZPhi :: Integral a => ZPhi a -> a
floorZPhi (ZPhi a b) = a + floorPhiTimes b
ceilingZPhi n = -floorZPhi (-n)