-- | Representation of Euclidean domains. That is integral domains with an 
-- Euclidean functions and decidable division.
--
module Algebra.Structures.EuclideanDomain 
  ( EuclideanDomain(..)
  , propEuclideanDomain
  , modulo, quotient, divides 
  , euclidAlg, genEuclidAlg
  , lcmE, genLcmE
  , extendedEuclidAlg, genExtendedEuclidAlg
  ) where

import Test.QuickCheck

import Algebra.Structures.IntegralDomain
-- import Algebra.Structures.Coherent
import Algebra.Ideal


-------------------------------------------------------------------------------
-- | Euclidean domains
--
-- Given a and b compute (q,r) such that a = bq + r and r = 0 || d r < d b. 
-- Where d is the Euclidean function.

class IntegralDomain a => EuclideanDomain a where
  d :: a -> Integer
  quotientRemainder :: a -> a -> (a,a)

-- Check both that |a| <= |ab| and |a| >= 0 for all a,b
propD :: (EuclideanDomain a, Eq a) => a -> a -> Bool
propD a b = 
  a == zero || b == zero || (d a <= d (a <*> b) && d a >= 0 && d b >= 0)

propQuotRem :: (EuclideanDomain a, Eq a) => a -> a -> Bool
propQuotRem a b = b == zero || (a == b <*> q <+> r && (r == zero || d r < d b))
  where (q,r) = quotientRemainder a b 

propEuclideanDomain :: (EuclideanDomain a, Eq a) => a -> a -> a -> Property
propEuclideanDomain a b c =
  if propD a b 
     then if propQuotRem a b
             then propIntegralDomain a b c
             else whenFail (print "propQuotRem") False
     else whenFail (print "propD") False


-------------------------------------------------------------------------------
-- Operations

modulo :: EuclideanDomain a => a -> a -> a
modulo a b = snd (quotientRemainder a b)

quotient :: EuclideanDomain a => a -> a -> a
quotient a b = fst (quotientRemainder a b)

divides :: (EuclideanDomain a, Eq a) => a -> a -> Bool
divides a b = modulo b a == zero

-- | The Euclidean algorithm for calculating the GCD of a and b.
euclidAlg :: (EuclideanDomain a, Eq a) => a -> a -> a
euclidAlg a b | a == zero && b == zero = error "GCD of 0 and 0 is undefined"
              | b == zero = a
              | otherwise = euclidAlg b (a `modulo` b)

-- | Generalized Euclidean algorithm to compute GCD of a list of elements.
genEuclidAlg :: (EuclideanDomain a, Eq a) => [a] -> a
genEuclidAlg = foldl euclidAlg zero

-- | Lowest common multiple, (a*b)/gcd(a,b).
lcmE :: (EuclideanDomain a, Eq a) => a -> a -> a
lcmE a b = quotient (a <*> b) (euclidAlg a b)

-- | Generalized lowest common multiple to compute lcm of a list of elements.
genLcmE :: (EuclideanDomain a, Eq a) => [a] -> a
genLcmE xs = quotient (foldr1 (<*>) xs) (genEuclidAlg xs)

-- | The extended Euclidean algorithm. 
-- 
-- Computes x and y in ax + by = gcd(a,b).
-- 
extendedEuclidAlg :: (EuclideanDomain a, Eq a) => a -> a -> (a,a)
extendedEuclidAlg a b | modulo a b == zero = (zero,one)
                      | otherwise          = (y, x <-> y <*> (a `quotient` b))
  where (x,y) = extendedEuclidAlg b (a `modulo` b)

-- Specification of extended Euclidean algorithm.
propExtendedEuclidAlg :: (EuclideanDomain a, Eq a) => a -> a -> Property
propExtendedEuclidAlg a b = a /= zero && b /= zero ==> 
  let (x,y) = extendedEuclidAlg a b in a <*> x <+> b <*> y == euclidAlg a b

-- | Generalized extended Euclidean algorithm.
--
-- Solves a_1 x_1 + ... + a_n x_n = gcd (a_1,...,a_n)
--
genExtendedEuclidAlg :: (EuclideanDomain a, Eq a) => [a] -> [a]
genExtendedEuclidAlg [x,y] = let (a,b) = extendedEuclidAlg x y in [a,b]
genExtendedEuclidAlg xs    =
  let (x,y) = extendedEuclidAlg (genEuclidAlg (init xs)) (last xs)
  in map (x<*>) (genExtendedEuclidAlg (init xs)) ++ [y]

-- Specification of generalized extended Euclidean algorithm. 
propGenExtEuclidAlg :: (EuclideanDomain a, Eq a) => [a] -> Property
propGenExtEuclidAlg xs = all (/= zero) xs && length xs >= 2 ==> 
  foldr (<+>) zero (zipWith (<*>) (genExtendedEuclidAlg xs) xs) == genEuclidAlg xs