{-# OPTIONS -fno-implicit-prelude #-} 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 NumericPrelude.Condition (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)