Internal.hs

{-# LANGUAGE NoImplicitPrelude, FlexibleInstances, UndecidableInstances, DefaultSignatures #-}
module Numeric.Domain.Internal where

import Data.Maybe(fromJust)
import Numeric.Additive.Group
import Numeric.Algebra.Class
import Numeric.Algebra.Commutative
import Numeric.Algebra.Division
import Numeric.Natural (Natural)
import Numeric.Semiring.ZeroProduct
import Numeric.Algebra.Unital.UnitNormalForm
import Numeric.Ring.Class
import Numeric.Decidable.Zero
import Numeric.Decidable.Units

import Prelude (Integer, Maybe (..), Bool(..),
                otherwise, fst, snd, ($), (.))
import qualified Prelude                 as P

infixl 7 `quot`, `rem`
infix  7 `divide`, `divides`, `maybeQuot`

-- | (Integral) domain is the integral semiring.
class (ZeroProductSemiring d, Ring d) => Domain d
instance (ZeroProductSemiring d, Ring d) => Domain d

-- | An integral domain is a commutative domain in which 1≠0.
class (Domain d, Commutative d) => IntegralDomain d where
    divides :: d -> d -> Bool
    default divides :: (Euclidean d) => d -> d -> Bool
    m `divides` n 
        | isZero m = False
        | otherwise = isZero (n `rem` m)
    maybeQuot :: d -> d -> Maybe d
    default maybeQuot :: (Euclidean d) => d -> d -> Maybe d
    m `maybeQuot` n
        | isZero n = Nothing
        | otherwise = let (q,r) = m `divide` n in
                      if isZero r then Just q else Nothing

instance IntegralDomain Integer

class (IntegralDomain d, UnitNormalForm d, DecidableZero d) => GCDDomain d where
    gcd :: d -> d -> d
    default gcd :: (PID d) => d -> d -> d
    gcd a b = let (r,_,_) = egcd a b in r
    {-# INLINE gcd #-}

    reduceFraction :: d -> d -> (d,d)
    reduceFraction a b =
        let c = gcd a b in
        (fromJust (a `maybeQuot` c), fromJust (b `maybeQuot` c))

    lcm :: d -> d -> d
    lcm p q = fromJust $ (p * q) `maybeQuot` (gcd p q)

instance GCDDomain Integer

class (GCDDomain d) => UFD d

instance UFD Integer

class (UFD d) => PID d where
    egcd :: d -> d -> (d,d,d)
    default egcd :: (Euclidean d) => d -> d -> (d,d,d)
    egcd a b = P.head (euclid a b)
    {-# INLINE egcd #-}

instance PID Integer

class (PID d) => Euclidean d where
  -- | Euclidean (degree) function on @r@.
  degree :: d -> Maybe Natural
  default degree :: (Division d) => d -> Maybe Natural
  degree a | isZero a = Nothing
           | otherwise = Just zero
  -- | Division algorithm. @a `divide` b@ calculates
  --   quotient and remainder of @a@ divided by @b@.
  --
  -- prop> let (q, r) = divide a p in p*q + r == a && degree r < degree q
  divide :: d                   -- ^ elements divided by
         -> d                   -- ^ divisor
         -> (d,d)               -- ^ quotient and remainder
  default divide :: (Division d) => d -> d -> (d,d)
  -- Be strict in order to make sure division by zero gets caught
  divide a b = let q = a/b in (q,P.seq q zero)

  quot :: d -> d -> d
  quot a b = fst $ a `divide` b
  {-# INLINE quot #-}

  rem :: d -> d -> d
  rem a b = snd $ a `divide` b
  {-# INLINE rem #-}

instance Euclidean Integer where
  degree = Just . P.fromInteger . P.abs
  {-# INLINE degree #-}

  divide = P.divMod
  {-# INLINE divide #-}


-- | Extended euclidean algorithm.
--
-- prop> euclid f g == xs ==> all (\(r, s, t) -> r == f * s + g * t) xs
euclid :: (Euclidean d) =>  d -> d -> [(d,d,d)]
euclid f g =
  let (ug, g') = splitUnit g
      Just t'  = recipUnit ug
      (uf, f') = splitUnit f
      Just s   = recipUnit uf
  in step [(g', zero, t'), (f', s, zero)]
  where
    step acc@((r',s',t'):(r,s,t):_)
      | isZero r' = P.tail acc
      | otherwise =
        let q         = r `quot` r'
            (ur, r'') = splitUnit $ r - q * r'
            Just u    = recipUnit ur
            s''       = (s - q * s') * u
            t''       = (t - q * t') * u
        in step ((r'', s'', t'') : acc)
    step _ = P.error "cannot happen!"