{-# LANGUAGE NoImplicitPrelude #-}
module Algebra.PrincipalIdealDomain (
{- * Class -}
C,
extendedGCD,
gcd,
lcm,
{- * 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, safeDiv, 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 PreludeBase
import Prelude (Integer, Int)
import qualified Prelude as P
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 = safeDiv x (gcd x y) * y -- avoid big temporary results
-- lcm x y = safeDiv (x * y) (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
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 = safeDiv x g
yl = safeDiv 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
-- cf. MathObj.Permutation.Table
swap :: (a,b) -> (b,a)
swap (x,y) = (y,x)
-}
{- |
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 exists a GCD variant,
that is specialised for integers and 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
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)