{-# LANGUAGE NoImplicitPrelude #-}
module Algebra.PrincipalIdealDomain (
    {- * Class -}
    C,
    extendedGCD,
    gcd,
    lcm,
    coprime,

    {- * Standard implementations for instances -}
    euclid,
    extendedEuclid,

    {- * Algorithms -}
    extendedGCDMulti,
    diophantine,
    diophantineMin,
    diophantineMulti,
    chineseRemainder,
    chineseRemainderMulti,

    {- * Properties -}
    propMaximalDivisor,
    propGCDDiophantine,
    propExtendedGCDMulti,
    propDiophantine,
    propDiophantineMin,
    propDiophantineMulti,
    propDiophantineMultiMin,
    propChineseRemainder,
    propDivisibleGCD,
    propDivisibleLCM,
    propGCDIdentity,
    propGCDCommutative,
    propGCDAssociative,
    propGCDHomogeneous,
    propGCD_LCM,
  ) where

import qualified Algebra.Units          as Units
import qualified Algebra.IntegralDomain as Integral
-- import qualified Algebra.Ring           as Ring
-- import qualified Algebra.Additive       as Additive
import qualified Algebra.ZeroTestable   as ZeroTestable

import qualified Algebra.Laws as Laws

import Algebra.Units          (stdAssociate, stdUnitInv)
import Algebra.IntegralDomain (mod, divChecked, divMod, divides, divModZero)
import Algebra.Ring           (one, (*), scalarProduct)
import Algebra.Additive       (zero, (+), (-))
import Algebra.ZeroTestable   (isZero)

import Data.Maybe.HT (toMaybe, )

import Control.Monad (foldM, liftM)
import Data.List (mapAccumL, mapAccumR, unfoldr)

import Data.Int  (Int,  Int8,  Int16,  Int32,  Int64,  )

import NumericPrelude.Base
import Prelude (Integer, )
import Test.QuickCheck ((==>), Property)



{- |
A principal ideal domain is a ring in which every ideal
(the set of multiples of some generating set of elements)
is principal:
That is,
every element can be written as the multiple of some generating element.
@gcd a b@ gives a generator for the ideal generated by @a@ and @b@.
The algorithm above works whenever @mod x y@ is smaller
(in a suitable sense) than both @x@ and @y@;
otherwise the algorithm may run forever.

Laws:

>   divides x (lcm x y)
>   x `gcd` (y `gcd` z) == (x `gcd` y) `gcd` z
>   gcd x y * z == gcd (x*z) (y*z)
>   gcd x y * lcm x y == x * y

(etc: canonical)

Minimal definition:
 * nothing, if the standard Euclidean algorithm work
 * if 'extendedGCD' is implemented customly, 'gcd' and 'lcm' make use of it
-}
class (Units.C a, ZeroTestable.C a) => C a where
    {- |
    Compute the greatest common divisor and
    solve a respective Diophantine equation.

    >   (g,(a,b)) = extendedGCD x y ==>
    >        g==a*x+b*y   &&  g == gcd x y

    TODO: This method is not appropriate for the PID class,
          because there are rings like the one of the multivariate polynomials,
          where for all x and y greatest common divisors of x and y exist,
          but they cannot be represented as a linear combination of x and y.
    TODO: The definition of extendedGCD does not return the canonical associate.
    -}
    extendedGCD :: a -> a -> (a,(a,a))
    extendedGCD = extendedEuclid divMod

    {- |
    The Greatest Common Divisor is defined by:

    >   gcd x y == gcd y x
    >   divides z x && divides z y ==> divides z (gcd x y)   (specification)
    >   divides (gcd x y) x
    -}
    gcd         :: a -> a -> a
    gcd x y     = fst $ extendedGCD x y

    {- |
    Least common multiple
    -}
    lcm         :: a -> a -> a
    lcm x y     =
       if isZero x
         then x -- avoid computing undefined (gcd 0 0)
         else divChecked x (gcd x y) * y  -- avoid big temporary results
    -- lcm x y     = divChecked (x * y) (gcd x y)


{-
These do only work if zero and one are really identity elements.

gcdMulti :: (C a) => [a] -> a
gcdMulti = foldl gcd zero

lcmMulti :: (C a) => [a] -> a
lcmMulti = foldl lcm one
-}

coprime :: (C a) => a -> a -> Bool
coprime x y =
   Units.isUnit (gcd x y)



{-
We could implement a helper function,
which exposes the temporary results.
This could be re-used for extendedEuclid.
-}
euclid :: (Units.C a, ZeroTestable.C a) =>
   (a -> a -> a) -> a -> a -> a
euclid genMod =
   let aux x y =
          if isZero y
            then stdAssociate x
            else aux y (x `genMod` y)
   in  aux

-- could be implemented in a tail-recursive manner
{-
Unfortunately, with the normalization to the stdAssociate,
@gcd 0@ is no longer the identity function,
since @gcd 0 (-2) = 2@.
-}
extendedEuclid :: (Units.C a, ZeroTestable.C a) =>
   (a -> a -> (a,a)) -> a -> a -> (a,(a,a))
extendedEuclid genDivMod =
   let aux x y =
          if isZero y
            then (stdAssociate x, (stdUnitInv x, zero))
            else
              let (d,m) = x `genDivMod` y   -- x == d*y + m
                  (g,(a,b)) = aux y m       -- g == a*y + b*m
              in  (g,(b,a-b*d))             -- g == a*y + b*(x-d*y)
   in aux


{- |
Compute the greatest common divisor for multiple numbers
by repeated application of the two-operand-gcd.
-}
extendedGCDMulti :: C a => [a] -> (a,[a])
extendedGCDMulti xs =
   let (g,cs) = mapAccumL extendedGCD zero xs
   in  (g, snd $ mapAccumR (\acc (c0,c1) -> (acc*c0,acc*c1)) one cs)

{- |
A variant with small coefficients.
-}


{- |
@Just (a,b) = diophantine z x y@
means
@a*x+b*y = z@.
It is required that @gcd(y,z) `divides` x@.
-}
diophantine :: C a => a -> a -> a -> Maybe (a,a)
diophantine z x y =
   fmap snd $ diophantineAux z x y

{- |
Like 'diophantine', but @a@ is minimal
with respect to the measure function of the Euclidean algorithm.
-}
diophantineMin :: C a => a -> a -> a -> Maybe (a,a)
diophantineMin z x y =
   fmap (uncurry (minimizeFirstOperand (x,y))) $
   diophantineAux z x y

minimizeFirstOperand :: C a => (a,a) -> a -> (a,a) -> (a,a)
minimizeFirstOperand (x,y) g (a,b) =
   if isZero g
     then (zero,zero)
     else
       let xl = divChecked x g
           yl = divChecked y g
           (d,aRed) = divModZero a yl
       in  (aRed, b + d*xl)

diophantineAux :: C a => a -> a -> a -> Maybe (a, (a,a))
diophantineAux z x y =
   let (g,(a,b)) = extendedGCD x y
       (q,r) = divModZero z g
   in  toMaybe (isZero r) (g, (q*a, q*b))


{- |
-}
diophantineMulti :: C a => a -> [a] -> Maybe [a]
diophantineMulti z xs =
   let (g,as) = extendedGCDMulti xs
       (q,r)  = divModZero z g
   in  toMaybe (isZero r) (map (q*) as)

{- |
Not efficient because it requires duplicate computations of GCDs.
However GCDs of neighbouring list elements were not computed before.
It is also quite arbitrary,
because only neighbouring elements are used for balancing.
There are certainly more sophisticated solutions.
-}
diophantineMultiMin :: C a => a -> [a] -> Maybe [a]
diophantineMultiMin z xs =
   do as <- diophantineMulti z xs
      return $ unfoldr
         (\as' -> case as' of
           ((x0,a0):(x1,a1):aRest) ->
              let (b0,b1) = minimizeFirstOperand (x0,x1) (gcd x0 x1) (a0,a1)
              in  Just (b0, (x1,b1):aRest)
           (_,a):[] -> Just (a,[])
           [] -> Nothing) $
         zip xs as

{-
diophantineMultiMin :: C a => a -> [a] -> Maybe [a]
diophantineMultiMin z xs =
   do as <- diophantineMulti z xs
      return $
         case as of
           (c:cs'@(_:_)) ->
              let (cs,cLast) = splitLast cs'
                  (d,as') = mapAccumL (\a b -> swap $ minimizeFirstOperand (gcd a b) (a,b)) c cs
                  (d',cLast') = minimizeFirstOperand (gcd d cLast) d cLast
              in  as' ++ [d',cLast']
           _ -> as
-}

{- |
Not efficient enough, because GCD\/LCM is computed twice.
-}
chineseRemainder :: C a => (a,a) -> (a,a) -> Maybe (a,a)
chineseRemainder (m0,a0) (m1,a1) =
   liftM (\(k,_) -> let m = lcm m0 m1 in (m, mod (a0-k*m0) m)) $
   diophantineMin (a0-a1) m0 m1
{-
a0-k*m0 == a1+l*m1
a0-a1 == k*m0+l*m1
-}

{- |
For @Just (b,n) = chineseRemainder [(a0,m0), (a1,m1), ..., (an,mn)]@
and all @x@ with @x = b mod n@ the congruences
@x=a0 mod m0, x=a1 mod m1, ..., x=an mod mn@
are fulfilled.
-}
chineseRemainderMulti :: C a => [(a,a)] -> Maybe (a,a)
chineseRemainderMulti congs =
   case congs of
      [] -> Nothing
      (c:cs) -> foldM chineseRemainder c cs



{- * Instances for atomic types -}


{-
There is the binary GCD algorithm,
that is specialised for integers in binary representation.
It does not need a division.
However, since we have an optimized division,
the standard implementation is probably faster.

TODO: Can Integer make use of the GMP GCD routine?
-}

instance C Integer where
    -- possibly more efficient than the default method
    gcd = euclid mod

instance C Int where
    gcd = euclid mod

instance C Int8 where
    gcd = euclid mod

instance C Int16 where
    gcd = euclid mod

instance C Int32 where
    gcd = euclid mod

instance C Int64 where
    gcd = euclid mod


propGCDIdentity     :: (Eq a, C a) => a -> Bool
propGCDAssociative :: (Eq a, C a) => a -> a -> a -> Bool
propGCDCommutative :: (Eq a, C a) => a -> a -> Bool
propGCDDiophantine :: (Eq a, C a) => a -> a -> Bool
propExtendedGCDMulti :: (Eq a, C a) => [a] -> Bool
propDiophantineGen :: (Eq a, C a) =>
   (a -> a -> a -> Maybe (a,a)) -> a -> a -> a -> a -> Bool
propDiophantine    :: (Eq a, C a) => a -> a -> a -> a -> Bool
propDiophantineMin :: (Eq a, C a) => a -> a -> a -> a -> Bool
propDiophantineMultiGen :: (Eq a, C a) =>
   (a -> [a] -> Maybe [a]) -> [(a,a)] -> Bool
propDiophantineMulti    :: (Eq a, C a) => [(a,a)] -> Bool
propDiophantineMultiMin :: (Eq a, C a) => [(a,a)] -> Bool
propDivisibleGCD   :: C a => a -> a -> Bool
propDivisibleLCM   :: C a => a -> a -> Bool
propGCD_LCM        :: (Eq a, C a) => a -> a -> Bool
propGCDHomogeneous :: (Eq a, C a) => a -> a -> a -> Bool
propMaximalDivisor :: C a => a -> a -> a -> Property
propChineseRemainder :: (Eq a, C a) => a -> a -> [a] -> Property

propMaximalDivisor x y z =
   divides z x && divides z y ==> divides z (gcd x y)
propGCDDiophantine x y =
   let (g,(a,b)) = extendedGCD x y
   in  g == gcd x y  &&  g == a*x+b*y
propExtendedGCDMulti xs =
   let (g,as) = extendedGCDMulti xs
   in  g == scalarProduct as xs  &&
       (isZero g || all (divides g) xs)
propDiophantineGen dio a b x y =
   let z = a*x+b*y
   in  maybe False (\(a',b') -> z == a'*x+b'*y) (dio z x y)
propDiophantine    = propDiophantineGen diophantine
propDiophantineMin = propDiophantineGen diophantineMin
propDiophantineMultiGen dio axs =
   let (as,xs) = unzip axs
       z = scalarProduct as xs
   in  maybe False (\as' -> z == scalarProduct as' xs) (dio z xs)
propDiophantineMulti    = propDiophantineMultiGen diophantineMulti
propDiophantineMultiMin = propDiophantineMultiGen diophantineMultiMin
propDivisibleGCD x y  =  divides (gcd x y) x
propDivisibleLCM x y  =  divides x (lcm x y)

propGCDIdentity     =  Laws.identity gcd zero . stdAssociate
propGCDCommutative  =  Laws.commutative gcd
propGCDAssociative  =  Laws.associative gcd
propGCDHomogeneous  =  Laws.leftDistributive (*) gcd . stdAssociate
propGCD_LCM x y     =  gcd x y * lcm x y == x * y
propChineseRemainder k x ms =
   not (null ms) && all (not . isZero) ms ==>
   -- cf. Useful.functionToGraph
   let congs = zip ms (map (mod x) ms)
   in  maybe False
          (\(mGlob,y) ->
             let yk = y+mGlob*k
             in  all (\(m,a) -> Integral.sameResidueClass m a yk) congs)
          (chineseRemainderMulti congs)