constructive-algebra 0.2.0 → 0.3.0

raw patch · 18 files changed

+668/−325 lines, 18 filesdep ~basePVP ok

version bump matches the API change (PVP)

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• README.hs +0/−80
• constructive-algebra.cabal +19/−16
• src/Algebra/Matrix.hs +197/−16
• src/Algebra/Q.hs +14/−11
• src/Algebra/Structures/BezoutDomain.hs +81/−66
• src/Algebra/Structures/Coherent.hs +38/−41
• src/Algebra/Structures/CommutativeRing.hs +2/−0
• src/Algebra/Structures/EuclideanDomain.hs +25/−21
• src/Algebra/Structures/ExplicitUnits.hs +34/−0
• src/Algebra/Structures/Field.hs +1/−1
• src/Algebra/Structures/FieldOfFractions.hs +16/−13
• src/Algebra/Structures/GCDDomain.hs +23/−4
• src/Algebra/Structures/Module.hs +8/−6
• src/Algebra/Structures/PruferDomain.hs +8/−3
• src/Algebra/Structures/Ring.hs +6/−1
• src/Algebra/UPoly.hs +164/−32
• src/Algebra/Z.hs +21/−3
• src/Algebra/Zn.hs +11/−11

@@ -1,80 +0,0 @@----------------------------------------------------------------------------------- | Constructive Algebra Library --- --- Anders Mortberg    <mortberg@student.chalmers.se>--- Bassel Mannaa      <mannaa@student.chalmers.se>------ Abstract:--- This is a library written as part of our master theses. It focuses mainly--- on the theory of commutative rings from a constructive point of view. -------------------------------------------------------------------------------------module README where-------------------------------------------------------------------------------------- Structures---- Rings with basic operations. -import Algebra.Structures.Ring---- Commutative rings.-import Algebra.Structures.CommutativeRing---- Integral domains.-import Algebra.Structures.IntegralDomain---- Fields.-import Algebra.Structures.Field---- Strongly discrete rings - Rings with decidable ideal membership.-import Algebra.Structures.StronglyDiscrete---- EuclideanDomains - Integral domains with decidable division and and Euclidean--- function. Contains lots of functions that are possible at the level of --- Euclidean domain like the Euclidean algorithm and extended Euclidean --- algorithm.-import Algebra.Structures.EuclideanDomain---- Bezout domains - Non-Noetherian analogues of principal ideal domains. All --- finitely generated ideals are principal.-import Algebra.Structures.BezoutDomain---- GCD domains - Non-Noetherian analogues of unique factorization domains. --- All pairs of nonzero elements have a greatest common divisor.-import Algebra.Structures.GCDDomain---- Field of fractions of a GCD domain.-import Algebra.Structures.FieldOfFractions---- Coherent rings. That is rings in which it is possible to solve homogenous--- linear equations. -import Algebra.Structures.Coherent------------------------------------------------------------------------------------- Special constructions.---- Finitely generated ideals over commutative rings. -import Algebra.Ideal---- Simple matrix library-import Algebra.Matrix---- Principle localization matrices-import Algebra.PLM------------------------------------------------------------------------------------- Instances.---- The integers.-import Algebra.Z---- The rational numbers as the field of fractions of Z. -import Algebra.Q------------------------------------------------------------------------------------- The end.

@@ -7,7 +7,7 @@ -- The package version. See the Haskell package versioning policy -- (http://www.haskell.org/haskellwiki/Package_versioning_policy) for -- standards guiding when and how versions should be incremented.-Version:             0.2.0+Version:             0.3.0  Synopsis:            A library of constructive algebra. Description:         @@ -19,7 +19,8 @@         generated. For example, instead of principal ideal domains one gets          Bezout domains which are integral domains in which all finitely          generated ideals are principal (and not necessarily that all ideals are-        principal).+        principal). This give a good framework for implementing many +        interesting algorithms.  License:             BSD3 License-file:        LICENSE@@ -40,7 +41,7 @@  -- Extra files to be distributed with the package, such as examples or -- a README.-Extra-source-files:  README.hs, examples/Z_Examples.hs+Extra-source-files:  examples/Z_Examples.hs  -- Constraint on the version of Cabal needed to build this package. Cabal-version:       >=1.2@@ -48,34 +49,36 @@  Library   -- Modules exported by the library.-  Exposed-modules:     Algebra.Structures.Group, -                       Algebra.Structures.Ring,+  Exposed-modules:     Algebra.Structures.BezoutDomain,+                       Algebra.Structures.Coherent,                        Algebra.Structures.CommutativeRing,-                       Algebra.Structures.IntegralDomain, +                       Algebra.Structures.EuclideanDomain,+                       Algebra.Structures.ExplicitUnits,                        Algebra.Structures.Field,+                       Algebra.Structures.FieldOfFractions,+                       Algebra.Structures.GCDDomain, +                       Algebra.Structures.Group, +                       Algebra.Structures.IntegralDomain,                         Algebra.Structures.Module,-                       Algebra.Structures.BezoutDomain,                        Algebra.Structures.PruferDomain,-                       Algebra.Structures.EuclideanDomain,+                       Algebra.Structures.Ring,                        Algebra.Structures.StronglyDiscrete,-                       Algebra.Structures.FieldOfFractions,-                       Algebra.Structures.GCDDomain, -                       Algebra.Structures.Coherent,                        Algebra.TypeChar.Char,+                       Algebra.EllipticCurve,                        Algebra.FieldOfRationalFunctions,                        Algebra.Ideal,                        Algebra.Matrix,                        Algebra.PLM,+                       Algebra.Q,+--                       Algebra.SmithNormalForm,  <- Should be added to separate package                        Algebra.UPoly,-                       Algebra.EllipticCurve,+                       Algebra.Z,                                          Algebra.ZSqrt5,-                       Algebra.Zn,-                       Algebra.Z,-                       Algebra.Q+                       Algebra.Zn                            -- Packages needed in order to build this package.-  Build-depends:       base >= 3 && <= 4, QuickCheck >= 2, type-level >= 0.2+  Build-depends:       base >= 3 && <= 4.3.1.0, QuickCheck >= 2, type-level >= 0.2      -- Modules not exported by this package.   -- Other-modules:       

@@ -4,16 +4,23 @@   ( Vector(Vec)   , unVec, lengthVec   , Matrix(M), matrix-  , matrixToVector, vectorToMatrix, unMVec, unM -  , identity, mulM, addM, transpose, isSquareMatrix, dimension+  , matrixToVector, vectorToMatrix, unMVec, unM, (!!!)+  , identity, propLeftIdentity, propRightIdentity+  , mulM, addM, transpose, isSquareMatrix, dimension+  , scale, swap, pivot+  , addRow, subRow, addCol, subCol+  , findPivot, forwardElim, gaussElim, gaussElimCorrect   ) where  import qualified Data.List as L+import Data.Function (on) import Control.Monad (liftM)+import Control.Arrow hiding ((<+>)) import Test.QuickCheck -import Algebra.Structures.IntegralDomain+import Algebra.Structures.Field +import Debug.Trace  ------------------------------------------------------------------------------- -- | Row vectors@@ -23,15 +30,19 @@ instance Show r => Show (Vector r) where   show (Vec vs) = show vs -instance (Ring r, Arbitrary r, Eq r) => Arbitrary (Vector r) where+-- Generate vector of length 1-10+instance Arbitrary r => Arbitrary (Vector r) where   arbitrary = do n <- choose (1,10) :: Gen Int                  liftM Vec $ gen n     where     gen 0 = return []     gen n = do x <- arbitrary                xs <- gen (n-1)-               if x == zero then return (one:xs) else return (x:xs)+               return (x:xs) +instance Functor Vector where +  fmap f = Vec . map f . unVec+ {- instance Ring r => Ring (Vector r) where   (Vec xs) <+> (Vec ys) | length xs == length ys = Vec (zipWith (<+>) xs ys)@@ -62,7 +73,8 @@     [] -> "[]"      xs -> init xs ++ "\n" -instance (Eq r, Arbitrary r, Ring r) => Arbitrary (Matrix r) where+-- Generate matrices with at most 10 rows+instance Arbitrary r => Arbitrary (Matrix r) where   arbitrary = do n <- choose (1,10) :: Gen Int                  m <- choose (1,10) :: Gen Int                  xs <- sequence [ liftM Vec (gen n) | _ <- [1..m]]@@ -71,8 +83,10 @@     gen 0 = return []     gen n = do x <- arbitrary                xs <- gen (n-1)-               if x == zero then return (one:xs) else return (x:xs)+               return (x:xs) +instance Functor Matrix where+  fmap f = M . map (fmap f) . unM  -- | Construct a mxn matrix. matrix :: [[r]] -> Matrix r@@ -96,11 +110,18 @@ matrixToVector m | fst (dimension m) == 1 = head (unM m)                  | otherwise              = error "matrixToVector: Bad dimension" +(!!!) :: Matrix a -> (Int,Int) -> a+m !!! (r,c) | r >= 0 && r < rows && c >= 0 && c < cols = unMVec m !! r !! c+            | otherwise = error "!!!: Out of bounds"+  where+  (rows,cols) = dimension m + -- | Compute the dimension of a matrix. dimension :: Matrix r -> (Int, Int) dimension (M xs) | null xs   = (0,0)                  | otherwise = (length xs, length (unVec (head xs))) + isSquareMatrix :: Matrix r -> Bool isSquareMatrix (M xs) = all (== length xs) (map lengthVec xs) @@ -110,7 +131,7 @@  -- | Matrix addition. addM :: Ring r => Matrix r -> Matrix r -> Matrix r-addM (M xs) (M ys) +addM (M xs) (M ys)   | dimension (M xs) == dimension (M ys) = m   | otherwise = error "Bad dimensions in matrix addition"   where@@ -118,33 +139,33 @@  -- | Matrix multiplication. mulM :: Ring r => Matrix r -> Matrix r -> Matrix r-mulM (M xs) (M ys) +mulM (M xs) (M ys)   | snd (dimension (M xs)) == fst (dimension (M ys)) = m   | otherwise = error "Bad dimensions in matrix multiplication"     where-    m = matrix [ [ foldr1 (<+>) (zipWith (<*>) x y) -                 | y <- L.transpose (map unVec ys) ]+    m = matrix [ [ mulVec x y | y <- L.transpose (map unVec ys) ]                | x <- map unVec xs ] -+mulVec xs ys | length xs == length ys = foldr (<+>) zero $ zipWith (<*>) xs ys+             | otherwise = error "mulVec: Bad dimension"  {- -- In order to do this the size of the matrix need to be encoded in the type--- There is also a problem with the fact that it is not possible to add or +-- There is also a problem with the fact that it is not possible to add or -- multiply matrices with bad dimensions, so the generation of matrices has to be better... instance Ring r => Ring (Matrix r) where   (<+>) = add   (<*>) = mul   neg (Vec xs d) = Vec [ map neg x | x <- xs ] d-  zero  = undefined +  zero  = undefined -}  -- | Construct a nxn identity matrix. identity :: IntegralDomain r => Int -> Matrix r identity n = matrix (xs 0)   where-  xs x | x == n    = [] -       | otherwise = (replicate x zero ++ [one] ++ +  xs x | x == n    = []+       | otherwise = (replicate x zero ++ [one] ++                       replicate (n-x-1) zero) : xs (x+1)  -- Specification of identity.@@ -155,3 +176,163 @@ propRightIdentity :: (IntegralDomain r, Eq r) => Matrix r -> Bool propRightIdentity a = a == a `mulM` identity m   where m = snd (dimension a)+++-------------------------------------------------------------------------------+-- Operations on matrices.++-- | Scale a row in a matrix.+scale :: CommutativeRing a => Matrix a -> Int -> a -> Matrix a+scale m r s+  | 0 <= r && r < rows = matrix $ take r m' ++ map (s <*>) (m' !! r) : drop (r+1) m'+  | otherwise = error "scale: Index out of bounds"+  where+  (rows,_) = dimension m+  m'       = unMVec m++-- Scaling does not affect dimension+propScaleDimension :: (Arbitrary r, CommutativeRing r) => Matrix r -> Int -> r -> Bool+propScaleDimension m r s = d == dimension (scale m (mod r rows) s)+  where d@(rows,_) = dimension m++-- | Swap two rows of a matrix.+swap :: Matrix a -> Int -> Int -> Matrix a+swap m i j+  | 0 <= i && i <= r && 0 <= j && j <= r = matrix $ swap' m' i j+  | otherwise = error "swap: Index out of bounds"+  where+  (r,_) = dimension m+  m'    = unMVec m++  swap' xs 0 0     = xs+  swap' (x:xs) 0 j = (x:xs) !! j : take (j-1) xs ++ x : drop j xs+  swap' xs i 0     = swap' xs 0 i+  swap' (x:xs) i j = x : swap' xs (i-1) (j-1)++-- Swapping does not affect dimension+propSwapDimension :: Matrix () -> Int -> Int -> Bool+propSwapDimension m i j = d == dimension (swap m (mod i r) (mod j r))+  where d@(r,_) = dimension m++-- Swap is itselfs identity.+propSwapIdentity :: Matrix () -> Int -> Int -> Bool+propSwapIdentity m i j = m == swap (swap m i' j') i' j'+  where+  d@(r,_) = dimension m+  i'      = mod i r+  j'      = mod j r+++-- Add the row-vector to the specified row of the matrix.+addRow :: CommutativeRing a => Matrix a -> Vector a -> Int -> Matrix a+addRow m row@(Vec xs) x+  | 0 <= x && x < r = matrix $ take x m' +++                               zipWith (<+>) (m' !! x) xs :+                               drop (x+1) m'+  | c /= length xs  = error "addRow: Bad length of row"+  | otherwise       = error "addRow: Bad row number"+    where+    (r,c) = dimension m+    m'    = unMVec m++propAddRowDimension :: (CommutativeRing a, Arbitrary a)+                    => Matrix a -> Vector a -> Int -> Property+propAddRowDimension m row@(Vec xs) r =+  length xs == c ==> d == dimension (addRow m row (mod r r'))+  where d@(r',c) = dimension m++addCol :: CommutativeRing a => Matrix a -> Vector a -> Int -> Matrix a+addCol m c x = transpose $ addRow (transpose m) c x++subRow, subCol :: CommutativeRing a => Matrix a -> Vector a -> Int -> Matrix a+subRow m (Vec xs) x = addRow m (Vec (map neg xs)) x+subCol m (Vec xs) x = addCol m (Vec (map neg xs)) x++-- Multiply the pivot row and add it to the target row.+pivot :: CommutativeRing a => Matrix a -> a -> Int -> Int -> Matrix a+pivot m s p t = addRow m (fmap (s <*>) (unM m !! p)) t++-- Find first non-zero number below the pivot and return its value and row number+-- given that it exists+findPivot :: (CommutativeRing a, Eq a) => Matrix a -> (Int,Int) -> Maybe (a,Int)+findPivot m (r,c) = safeHead $ filter ((/= zero) . fst) $ drop (r+1) $ zip (head $ drop c $ unMVec $ transpose m) [0..]+  where+  m' = unMVec m++  safeHead []     = Nothing+  safeHead (x:xs) = Just x++fE :: (Field a, Eq a) => Matrix a -> Matrix a+fE (M [])         = M []+fE (M (Vec []:_)) = M []+fE m     = case L.findIndices (/= zero) (map head xs) of+  (i:is) -> case fE (cancelOut m [ (i,map head xs !! i) | i <- is ] (i,map head xs !! i)) of+    ys -> matrix (xs !! i : map (zero :) (unMVec ys))+  []     -> case fE (matrix (map tail xs)) of+    ys -> matrix (map (zero:) (unMVec ys))+  where+  cancelOut :: (Field a, Eq a) => Matrix a -> [(Int,a)] -> (Int,a) -> Matrix a+  cancelOut m [] (i,_)    = let xs = unMVec m in matrix $ map tail (L.delete (xs !! i) xs)+  cancelOut m ((t,x):xs) (i,p) = cancelOut (pivot m (neg (x </> p)) i t) xs (i,p)++  xs = unMVec m+++-- | Compute row echelon form of a system Ax=b.+forwardElim :: (Field a, Eq a) => (Matrix a,Vector a) -> (Matrix a,Vector a)+forwardElim (m,v) = fE m' (0,0)+  where+  -- fE takes the matrix to eliminate and the current row and column+  fE :: (Field a, Eq a) => Matrix a -> (Int,Int) -> (Matrix a,Vector a)+  fE (M []) _  = error "forwardElim: Empty input matrix"+  fE m rc@(r,c)+      -- The algorithm is done when it reaches the last column or row.+    | c == mc || r == mr =+      -- Decompose the matrix into A and b again+      (matrix *** Vec) $ unzip $ map (init &&& last) $ unMVec m++    | m !!! rc == zero   = case findPivot m rc of+      -- If the pivot element is zero swap the pivot row with the first row+      -- with a nonzero element in the pivot column.+      Just (_,r') -> fE (swap m r r') rc+      -- If all elements in the pivot column is zero the move right.+      Nothing     -> fE m (r,c+1)++    | m !!! rc /= one    =+      -- Make the pivot element 1.+      fE (scale m r (inv (m !!! rc))) rc++    | otherwise          = case findPivot m rc of+      -- Make the first nonzero element in the pivot row 0.+      Just (v,r') -> fE (pivot m (neg v) r r') (r,c)+      -- If all elements in the pivot column is zero then move down and right.+      Nothing     -> fE m (r+1,c+1)++  (mr,mc) = dimension m++  -- Combine A and b to a matrix where the last column is b+  m' = matrix $ [ r ++ [x] | (r,x) <- zip (unMVec m) (unVec v) ]+++-- | Perform "jordan"-step in Gauss-Jordan elimination. That is make every+-- element above the diagonal zero. In other words compute the reduced+-- echelon form of a matrix given that the input is in row echelon form.+jordan :: (Field a, Eq a) => (Matrix a, Vector a) -> (Matrix a, Vector a)+jordan (m, Vec ys) = case L.unzip (jordan' (zip (unMVec m) ys) (r-1)) of+  (a,b) -> (matrix a, Vec b)+  where+  (r,_) = dimension m++  jordan' [] _ = []+  jordan' xs c =+    jordan' [ (take c x ++ zero : drop (c+1) x, v <-> x !! c <*> snd (last xs))+            | (x,v) <- init xs ] (c-1) ++ [last xs]+++-- | Gauss-Jordan elimination: Given A and B solve Ax=B.+gaussElim :: (Field a, Eq a, Show a) => (Matrix a, Vector a) -> (Matrix a, Vector a)+gaussElim = jordan . forwardElim++gaussElimCorrect :: (Field a, Eq a, Arbitrary a, Show a) => (Matrix a, Vector a) -> Property+gaussElimCorrect m@(a,b) = fst (dimension a) == lengthVec b && isSquareMatrix a ==>+  matrixToVector (transpose (a `mulM` transpose (M [snd (gaussElim m)]))) == b

@@ -1,16 +1,13 @@ {-# LANGUAGE TypeSynonymInstances #-} -- | Representation of rational numbers as the field of fractions of Z.-module Algebra.Q -  ( Q-  , toQ, toZ-  ) where+module Algebra.Q ( Q, toQ, toZ ) where +import Data.Ratio (numerator, denominator) import Test.QuickCheck--- import qualified Math.Algebra.Field.Base as A (Q(..)) --- import Data.Ratio   import Algebra.Structures.Field-import Algebra.Structures.FieldOfFractions+import Algebra.Structures.FieldOfFractions hiding (numerator, denominator)+import Algebra.Structures.ExplicitUnits import Algebra.Z  -------------------------------------------------------------------------------@@ -26,12 +23,18 @@   fromInteger      = toQ  instance Fractional Q where-  (/) = (</>)-  fromRational = undefined---   fromRational (a :% b) = reduce $ F (a,b)+  (/)            = (</>)+  fromRational x = +    reduce $ F (fromIntegral (numerator x), fromIntegral (denominator x))  toQ :: Z -> Q toQ = toFieldOfFractions  toZ :: Q -> Z-toZ = fromFieldOfFractions   +toZ = fromFieldOfFractions++propFieldQ :: Q -> Q -> Q -> Property+propFieldQ = propField++instance ExplicitUnits Q where+  unit a = if a == 0 then Nothing else Just (inv a)

@@ -1,18 +1,16 @@ {-# LANGUAGE FlexibleInstances, UndecidableInstances #-}--- | Representation of Bezout domains. That is non-Noetherian analogues of +-- | Representation of Bezout domains. That is non-Noetherian analogues of -- principal ideal domains. This means that all finitely generated ideals are -- principal. -- module Algebra.Structures.BezoutDomain-  ( BezoutDomain(..)-  , propToPrincipal, propIsSameIdeal, propBezoutDomain-  , dividesB, gcdB-  , intersectionB, intersectionBWitness-  , solveB-  , crt+  ( BezoutDomain(..), propBezoutDomain+  , toPrincipal, propToPrincipal, propIsSameIdeal+  , gcdB, intersectionB, intersectionBWitness+  , solveB, crt   ) where -import Test.QuickCheck +import Test.QuickCheck  import Algebra.Structures.IntegralDomain import Algebra.Structures.Coherent@@ -23,21 +21,49 @@   ---------------------------------------------------------------------------------- | Bezout domains--- --- Compute a principal ideal from another ideal. Also give witness that the--- principal ideal is equal to the first ideal.------ toPrincipal \<a_1,...,a_n> = (\<a>,u_i,v_i)---   where------   sum (u_i * a_i) = a------   a_i = v_i * a---+{- | Bezout domains++Has a Bezout function which given a and b give g, a1, b1, x and y such that:++ - g = gcd(a,b)++ - a = g * a1 and b = g * b1++ - g = a * x + b * y++-} class IntegralDomain a => BezoutDomain a where-  toPrincipal :: Ideal a -> (Ideal a,[a],[a])+  bezout :: a -> a -> (a,a,a,a,a) +propBezoutDomain :: (BezoutDomain a, Eq a) => a -> a -> Property+propBezoutDomain a b = +  let (g,a1,b1,x,y) = bezout a b+  in if a == g <*> a1 && b == g <*> b1 && a <*> x <+> b <*> y == g+        then propIntegralDomain a b b+        else whenFail (print "propBezoutDomain") False+++{- | Compute a principal ideal from another ideal. Also give witness that the+principal ideal is equal to the first ideal.++toPrincipal \<a_1,...,a_n> = (\<a>,u_i,v_i)+  where++  sum (u_i * a_i) = a++  a_i = v_i * a+-}+toPrincipal :: BezoutDomain a => Ideal a -> (Ideal a,[a],[a])+toPrincipal (Id [])    = error "toPrincipal: Empty input"+toPrincipal (Id [a])   = (Id [a],[one],[one])+toPrincipal (Id [a,b]) = +  let (g,a1,b1,x,y) = bezout a b+  in (Id [g],[x,y],[a1,b1])+toPrincipal (Id (a:xs)) = +  let (Id [g],us,vs) = toPrincipal (Id xs)+      (g',a1,b1,x,y) = bezout a g+  in (Id [g'],x : map (y <*>) us,a1 : map (b1 <*>) vs)+ -- | Test that the generated ideal is principal. propToPrincipal :: (BezoutDomain a, Eq a) => Ideal a -> Bool propToPrincipal = isPrincipal . (\(a,_,_) -> a) . toPrincipal@@ -45,24 +71,12 @@ -- | Test that the generated ideal generate the same elements as the given. propIsSameIdeal :: (BezoutDomain a, Eq a) => Ideal a -> Bool propIsSameIdeal (Id as) =-  let (Id [a], us, vs) = toPrincipal (Id as) -  in a == foldr1 (<+>) (zipWith (<*>) as us) +  let (Id [a], us, vs) = toPrincipal (Id as)+  in a == foldr1 (<+>) (zipWith (<*>) as us)   && and [ ai == a <*> vi | (ai,vi) <- zip as vs ]   && length us == l_as && length vs == l_as   where l_as = length as -propBezoutDomain :: (BezoutDomain a, Eq a) => Ideal a -> a -> a -> a -> Property-propBezoutDomain id@(Id xs) a b c = zero `notElem` xs ==> -  if propToPrincipal id-     then if propIsSameIdeal id-             then propIntegralDomain a b c -             else whenFail (print "propIsSameIdeal") False-     else whenFail (print "propToPrincipal") False--dividesB :: (BezoutDomain a, Eq a) => a -> a -> Bool-dividesB a b = a == x || a == neg x-    where (Id [x],_,_) = toPrincipal (Id [a,b])- -- TODO: Add error cases... gcdB :: BezoutDomain a => a -> a -> a gcdB a b = g@@ -72,24 +86,28 @@ -- Euclidean domain -> Bezout domain  instance (EuclideanDomain a, Eq a) => BezoutDomain a where-  toPrincipal (Id [x]) = (Id [x], [one], [one])-  toPrincipal (Id xs)  = (Id [a], as, [ quotient ai a | ai <- xs ])-    where-    a  = genEuclidAlg xs-    as = genExtendedEuclidAlg xs+  bezout a b+    | b == zero        = (a,one,zero,one,zero)+    | a == zero        = (b,zero,one,zero,one)+    | norm a <= norm b = let (q,r)           = quotientRemainder b a+                             (g,a1,r1,u',v') = bezout a r+                         in (g,a1,r1<+>(q<*>a1),u'<->(q<*>v'),v')+    | otherwise        = let (q,r)           = quotientRemainder a b+                             (g,b1,r1,u',v') = bezout b r+                         in (g,r1<+>(q<*>b1),b1,v',u'<->(v'<*>q))   ------------------------------------------------------------------------------- -- | Intersection of ideals with witness.--- --- If one of the ideals is the zero ideal then the intersection is the zero +--+-- If one of the ideals is the zero ideal then the intersection is the zero -- ideal.--- -intersectionBWitness :: (BezoutDomain a, Eq a) -              => Ideal a -              -> Ideal a ++intersectionBWitness :: (BezoutDomain a, Eq a)+              => Ideal a+              -> Ideal a               -> (Ideal a, [[a]], [[a]])-intersectionBWitness (Id xs) (Id ys) +intersectionBWitness (Id xs) (Id ys)   | xs' == [] = zeroIdealWitnesses xs ys   | ys' == [] = zeroIdealWitnesses xs ys   | otherwise = (Id [l], [handleZero xs as], [handleZero ys bs])@@ -97,7 +115,7 @@   xs'            = filter (/= zero) xs   ys'            = filter (/= zero) ys -  (Id [a],us1,vs1) = toPrincipal (Id xs') +  (Id [a],us1,vs1) = toPrincipal (Id xs')   (Id [b],us2,vs2) = toPrincipal (Id ys')    (Id [g],[u1,u2],[v1,v2]) = toPrincipal (Id [a,b])@@ -105,14 +123,14 @@   l  = g <*> v1 <*> v2   as = map (v2 <*>) us1   bs = map (v1 <*>) us2-  + -- Handle the zeroes specially. If the first element in xs is a zero--- then the witness should be zero otherwise use the computed witness. +-- then the witness should be zero otherwise use the computed witness. handleZero :: (Ring a, Eq a) => [a] -> [a] -> [a]-handleZero xs [] +handleZero xs []   | all (==zero) xs = xs   | otherwise       = error "intersectionB: This should be impossible"-handleZero (x:xs) (a:as) +handleZero (x:xs) (a:as)   | x == zero = zero : handleZero xs (a:as)   | otherwise = a    : handleZero xs as handleZero [] _  = error "intersectionB: This should be impossible"@@ -134,13 +152,13 @@  ------------------------------------------------------------------------------- -- | Strongly discreteness for Bezout domains--- +-- -- Given x, compute as such that x = sum (a_i * x_i) -- instance (BezoutDomain a, Eq a) => StronglyDiscrete a where   member x (Id xs) | x == zero = Just (replicate (length xs) zero)-                   | otherwise = if a == g -                                    then Just witness +                   | otherwise = if a == g+                                    then Just witness                                     else Nothing     where     -- (<g>, as, bs)   = <x1,...,xn>@@ -148,30 +166,27 @@     -- x_i             = b_i * g     (Id [g], as, bs) = toPrincipal (Id (filter (/= zero) xs))     (Id [a], _,[q1,q2]) = toPrincipal (Id [x,g])-    +     -- x = qg = q (sum (ai * xi)) = sum (q * ai * xi)     witness = handleZero xs (map (q1 <*>) as)   ------------------------------------------------------------------------------- -- | Chinese remainder theorem--- +-- -- Given a_1,...,a_n and m_1,...,m_n such that gcd(m_i,m_j) = 1. -- Let m = m_1*...*m_n compute a such that:--- --- (1) a = a_i (mod m_i) --- --- (2) If b is such that---      ---        b = a_i (mod m_i)--- ---     then a = b (mod m) --+-- (1) a = a_i (mod m_i)+--+-- (2) If b is such that b = a_i (mod m_i) then a = b (mod m)+-- -- The function return (a,m).+ crt :: (BezoutDomain a, Eq a) => [a] -> [a] -> (a,a) crt as ms   | length as /= length ms = error "crt: Input lists need to have same length"-  | not (and [ gcdB m1 m2 == one | m1 <- ms, m2 <- ms, m1 /= m2 ]) = +  | not (and [ gcdB m1 m2 == one | m1 <- ms, m2 <- ms, m1 /= m2 ]) =       error "crt: All ms need to be relatively prime"   | otherwise = crt' as ms   where

@@ -1,9 +1,9 @@ -- | Representation of coherent rings. Traditionally a ring is coherent if every--- finitely generated ideal is finitely presented. This means that it is +-- finitely generated ideal is finitely presented. This means that it is -- possible to solve homogenous linear equations in them. module Algebra.Structures.Coherent   ( Coherent(solve)-  , propCoherent+  , propCoherent, isSolution   , solveMxN, propSolveMxN   , solveWithIntersection   , solveGeneralEquation, propSolveGeneralEquation@@ -25,7 +25,7 @@ -- --   MX=0   \<-\>  \exists Y. X=LY ----- that is, iff we can generate the solutions of any linear homogeous system +-- that is, iff we can generate the solutions of any linear homogeous system -- of equations. -- -- The main point here is that ML=0, it is not clear how to represent the@@ -34,6 +34,7 @@ class IntegralDomain a => Coherent a where   solve :: Vector a -> Matrix a +-- | Test that the second matrix is a solution to the first. isSolution :: (CommutativeRing a, Eq a) => Matrix a -> Matrix a -> Bool isSolution m sol = all (==zero) (concat (unMVec (m `mulM` sol))) @@ -45,14 +46,14 @@ solveMxN :: (Coherent a, Eq a) => Matrix a -> Matrix a solveMxN (M (l:ls)) = solveMxN' (solve l) ls   where-  -- Inductively solve all subsystems. If the computed solution is in fact a -  -- solution to the next set of equations then don't do anything. +  -- Inductively solve all subsystems. If the computed solution is in fact a+  -- solution to the next set of equations then don't do anything.   -- This solves the problems with having many identical rows in the system,   -- like [[1,1],[1,1]].   solveMxN' :: (Coherent a, Eq a) => Matrix a -> [Vector a] -> Matrix a   solveMxN' m []      = m-  solveMxN' m1 (x:xs) = if isSolution (vectorToMatrix x) m1 -                           then solveMxN' m1 xs +  solveMxN' m1 (x:xs) = if isSolution (vectorToMatrix x) m1+                           then solveMxN' m1 xs                            else solveMxN' (m1 `mulM` m2) xs     where m2 = solve (matrixToVector (mulM (vectorToMatrix x) m1)) @@ -64,38 +65,34 @@  ------------------------------------------------------------------------------- -- | Intersection computable -> Coherence.--- --- Proof that if there is an algorithm to compute a f.g. set of generators for +--+-- Proof that if there is an algorithm to compute a f.g. set of generators for -- the intersection of two f.g. ideals then the ring is coherent. ----- Takes the vector to solve, \[x1,...,xn\], and a function (int) that computes --- the intersection of two ideals. +-- Takes the vector to solve, \[x1,...,xn\], and a function (int) that computes+-- the intersection of two ideals. ----- If---       ---     \[ x_1, ..., x_n \] \`int\` \[ y_1, ..., y_m \] = \[ z_1, ..., z_l \]+-- If \[ x_1, ..., x_n \] \`int\` \[ y_1, ..., y_m \] = \[ z_1, ..., z_l \] -- -- then int should give witnesses us and vs such that: -----     z_k = n_k1 * x_1 + ... + u_kn * x_n------         = u_k1 * y_1 + ... + n_km * y_m+--     z_k = n_k1 * x_1 + ... + u_kn * x_n = u_k1 * y_1 + ... + n_km * y_m -- solveWithIntersection :: (IntegralDomain a, Eq a)-                      => Vector a +                      => Vector a                       -> (Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]]))                       -> Matrix a-solveWithIntersection (Vec xs) int = transpose $ matrix $ solveInt xs +solveWithIntersection (Vec xs) int = transpose $ matrix $ solveInt xs   where   solveInt []     = error "solveInt: Can't solve an empty system"   solveInt [x]    = [[zero]] -- Base case, could be [x,y] also...-                             -- That wouldn't give the trivial solution... ---  solveInt [x,y]  | x == zero || y == zero = [[zero,zero]]---                  | otherwise = ---    let (Id ts,us,vs) = (Id [x]) `int` (Id [neg y])---    in [ u ++ v | (u,v) <- zip us vs ]+                             -- That wouldn't give the trivial solution...+  solveInt [x,y]  | x == zero || y == zero = [[zero,zero]]+                  | otherwise =+    let (Id ts,us,vs) = (Id [x]) `int` (Id [neg y])+    in [ u ++ v | (u,v) <- zip us vs ]   solveInt (x:xs)-    | x == zero             = map (zero:) $ solveInt xs +    | x == zero             = (one : replicate (length xs) zero) : (map (zero:) $ solveInt xs)     | isSameIdeal int as bs = s ++ m'     | otherwise             = error "solveInt: This does not compute the intersection"       where@@ -105,7 +102,7 @@       -- Compute the intersection of <x1> and <-x2,...,-xn>       (Id ts,us,vs) = as `int` bs       s             = [ u ++ v | (u,v) <- zip us vs ]-      +       -- Solve <0,x2,...,xn> recursively       m             = solveInt xs       m'            = map (zero:) m@@ -114,11 +111,11 @@ ------------------------------------------------------------------------------- -- | Strongly discrete coherent rings. ----- If the ring is strongly discrete and coherent then we can solve arbitrary --- equations of the type AX=b. +-- If the ring is strongly discrete and coherent then we can solve arbitrary+-- equations of the type AX=b. -- solveGeneralEquation :: (Coherent a, StronglyDiscrete a) => Vector a -> a -> Maybe (Matrix a)-solveGeneralEquation v@(Vec xs) b = +solveGeneralEquation v@(Vec xs) b =   let sol = solve v   in case b `member` (Id xs) of     Just as -> Just $ transpose (M (replicate (length (head (unMVec sol))) (Vec as)))@@ -126,9 +123,9 @@     Nothing -> Nothing  -propSolveGeneralEquation :: (Coherent a, StronglyDiscrete a, Eq a) -                         => Vector a -                         -> a +propSolveGeneralEquation :: (Coherent a, StronglyDiscrete a, Eq a)+                         => Vector a+                         -> a                          -> Bool propSolveGeneralEquation v b = case solveGeneralEquation v b of   Just sol -> all (==b) $ concat $ unMVec $ vectorToMatrix v `mulM` sol@@ -140,24 +137,24 @@  -- | Solves general linear systems of the kind AX = B. ----- A is given as a matrix and B is given as a row vector (it should be column +-- A is given as a matrix and B is given as a row vector (it should be column -- vector). ---solveGeneral :: (Coherent a, StronglyDiscrete a, Eq a) +solveGeneral :: (Coherent a, StronglyDiscrete a, Eq a)              => Matrix a   -- M              -> Vector a   -- B              -> Maybe (Matrix a, Matrix a)  -- (L,X0)-solveGeneral (M (l:ls)) (Vec (a:as)) = +solveGeneral (M (l:ls)) (Vec (a:as)) =   case solveGeneral' (solveGeneralEquation l a) ls as [(l,a)] of     Just x0 -> Just (solveMxN (M (l:ls)), x0)     Nothing -> Nothing   where-  -- Compute a new solution inductively and check that the new solution +  -- Compute a new solution inductively and check that the new solution   -- satisfies all the previous equations.   solveGeneral' Nothing _ _ _              = Nothing-  solveGeneral' (Just m) [] [] old         = Just m -  solveGeneral' (Just m) (l:ls) (a:as) old = -    if isSolutionB l m a +  solveGeneral' (Just m) [] [] old         = Just m+  solveGeneral' (Just m) (l:ls) (a:as) old =+    if isSolutionB l m a        then solveGeneral' (Just m) ls as old        else case solveGeneralEquation (matrixToVector (vectorToMatrix l `mulM` m)) a of          Just m' -> let m'' = m `mulM` m'@@ -165,10 +162,10 @@                           then solveGeneral' (Just m'') ls as ((l,a):old)                           else Nothing          Nothing -> Nothing-  solveGeneral' _ _ _ _ = error "solveGeneral: Bad input"      +  solveGeneral' _ _ _ _ = error "solveGeneral: Bad input"  -- It would be great to only generate solvable systems...--- propSolveGeneral :: (Coherent a, StronglyDiscrete a, Eq a) => Matrix a -> Vector a -> Property +-- propSolveGeneral :: (Coherent a, StronglyDiscrete a, Eq a) => Matrix a -> Vector a -> Property propSolveGeneral m b = length (unM m) == length (unVec b) ==> case solveGeneral m b of   Just (l,x) -> all (==b) (unM (transpose (m `mulM` x))) &&                 isSolution m l

@@ -1,3 +1,5 @@+-- | Structure for commutative rings. +-- module Algebra.Structures.CommutativeRing   ( module Algebra.Structures.Ring   , CommutativeRing(..)

@@ -1,10 +1,10 @@--- | Representation of Euclidean domains. That is integral domains with an +-- | Representation of Euclidean domains. That is integral domains with an -- Euclidean functions and decidable division. ---module Algebra.Structures.EuclideanDomain +module Algebra.Structures.EuclideanDomain   ( EuclideanDomain(..)-  , propD, propQuotRem, propEuclideanDomain-  , modulo, quotient, divides +  , propNorm, propQuotRem, propEuclideanDomain+  , modulo, quotient, divides   , euclidAlg, genEuclidAlg   , lcmE, genLcmE   , extendedEuclidAlg, genExtendedEuclidAlg@@ -20,25 +20,27 @@ ------------------------------------------------------------------------------- -- | 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.+-- Given a and b compute (q,r) such that a = bq + r and r = 0 || norm r < norm b.+-- Where norm is the Euclidean function.  class IntegralDomain a => EuclideanDomain a where-  d :: a -> Integer+  norm :: 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)+propNorm :: (EuclideanDomain a, Eq a) => a -> a -> Bool+propNorm a b =+  a == zero || b == zero || +  (norm a <= norm (a <*> b) && norm a >= 0 && norm 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 +propQuotRem a b = +  b == zero || (a == b <*> q <+> r && (r == zero || norm r < norm b))+    where (q,r) = quotientRemainder a b  propEuclideanDomain :: (EuclideanDomain a, Eq a) => a -> a -> a -> Property propEuclideanDomain a b c =-  if propD a b +  if propNorm a b      then if propQuotRem a b              then propIntegralDomain a b c              else whenFail (print "propQuotRem") False@@ -59,7 +61,7 @@  -- | 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"+euclidAlg a b | a == zero && b == zero = zero               | b == zero = a               | otherwise = euclidAlg b (a `modulo` b) @@ -75,18 +77,20 @@ genLcmE :: (EuclideanDomain a, Eq a) => [a] -> a genLcmE xs = quotient (foldr1 (<*>) xs) (genEuclidAlg xs) --- | The extended Euclidean algorithm. --- +-- | 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)+extendedEuclidAlg a b | a == zero = (zero,one)+                      | b == zero = (one,zero)+--                      | modulo a b == zero = (one,zero)                       | 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 ==> +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.@@ -99,7 +103,7 @@   let (x,y) = extendedEuclidAlg (genEuclidAlg (init xs)) (last xs)   in map (x<*>) (genExtendedEuclidAlg (init xs)) ++ [y] --- Specification of generalized extended Euclidean algorithm. +-- Specification of generalized extended Euclidean algorithm. propGenExtEuclidAlg :: (EuclideanDomain a, Eq a) => [a] -> Property-propGenExtEuclidAlg xs = all (/= zero) xs && length xs >= 2 ==> +propGenExtEuclidAlg xs = all (/= zero) xs && length xs >= 2 ==>   foldr (<+>) zero (zipWith (<*>) (genExtendedEuclidAlg xs) xs) == genEuclidAlg xs

@@ -0,0 +1,34 @@+-- | Structure of rings with explicit units. +module Algebra.Structures.ExplicitUnits+  ( ExplicitUnits(..)+  , propUnit, isUnit, (%|), (~=)+  ) where ++import Algebra.Structures.IntegralDomain+import Algebra.Structures.GCDDomain++infix 5 %|+infix 4 ~=++-- | A ring has explicit units if there is a function that can test if an+-- element is invertible and if this is the case give the inverse. +class IntegralDomain a => ExplicitUnits a where+  unit :: a -> Maybe a++propUnit :: (ExplicitUnits a, Eq a) => a -> Bool+propUnit a = case unit a of+  Just a' -> a <*> a' == one+  Nothing -> True++-- | An element is a unit if it is invertible. +isUnit :: ExplicitUnits a => a -> Bool+isUnit = maybe False (const True) . unit++-- | Decidable units is sufficient to decide divisibility in GCD domains. +(%|) :: (ExplicitUnits a, GCDDomain a) => a -> a -> Bool+a %| b = let (g,x,y) = gcd' a b+         in isUnit x++-- | Test for associatedness, i.e. a ~ b iff a | b /\\ b | a.+(~=) :: (ExplicitUnits a, GCDDomain a) => a -> a -> Bool+a ~= b = a %| b && b %| a 

@@ -1,3 +1,4 @@+-- | Structure for fields. module Algebra.Structures.Field   ( module Algebra.Structures.IntegralDomain   , Field(inv)@@ -9,7 +10,6 @@  import Algebra.Structures.Ring import Algebra.Structures.IntegralDomain-  infixl 7 </> 

@@ -2,6 +2,7 @@ -- domain is that we only want to work over reduced quotients. module Algebra.Structures.FieldOfFractions   ( FieldOfFractions(..)+  , numerator, denominator   , toFieldOfFractions, fromFieldOfFractions   , reduce, propReduce   ) where@@ -17,7 +18,11 @@  newtype GCDDomain a => FieldOfFractions a = F (a,a) +numerator, denominator :: GCDDomain a => FieldOfFractions a -> a+numerator (F (x,_))   = x+denominator (F (_,x)) = x + -------------------------------------------------------------------------------- -- Instances @@ -33,13 +38,13 @@     b <- arbitrary     if b == zero         then return $ F (a,one)-       else return $ F (a,b)+       else return $ reduce $ F (a,b)  instance (GCDDomain a, Eq a) => Eq (FieldOfFractions a) where-  f == g = a <*> d == b <*> c-    where-    F (a,b) = reduce f-    F (c,d) = reduce g+  (F (a,b)) == (F (c,d)) = a <*> d == b <*> c+--    where+--    F (a,b) = reduce f+--    F (c,d) = reduce g  instance (GCDDomain a, Eq a) => Ring (FieldOfFractions a) where   (F (a,b)) <+> (F (c,d)) = reduce (F (a <*> d <+> c <*> b,b <*> d))@@ -52,7 +57,7 @@ instance (GCDDomain a, Eq a) => IntegralDomain (FieldOfFractions a)  instance (GCDDomain a, Eq a) => Field (FieldOfFractions a) where-  inv (F (a,b)) | b /= zero && a /= zero = reduce $ F (b,a)+  inv (F (a,b)) | a /= zero && b /= zero = reduce $ F (b,a)                 | otherwise = error "FieldOfFraction: Division by zero"  @@ -66,17 +71,15 @@ -- | Extract a value from the field of fractions. This is only possible if the -- divisor is one. fromFieldOfFractions :: (GCDDomain a, Eq a) => FieldOfFractions a -> a-fromFieldOfFractions (F (a,b)) -  | b == one  = a-  | otherwise = error "FieldOfFractions: Can't extract value"+fromFieldOfFractions x | b' == one = a'+                       | otherwise = error "fromFieldOfFractions: Division by zero"+  where F (a',b') = reduce x   -- | Reduce an element. reduce :: (GCDDomain a, Eq a) => FieldOfFractions a -> FieldOfFractions a-reduce (F (a,b)) | b == zero = error "FieldOfFractions: Division by zero"+reduce (F (a,b)) | b == zero = error "reduce: Division by zero"                  | a == zero = F (zero,one)-                 | otherwise = if g == one-                                  then F (a,b)-                                  else F (x,y)+                 | otherwise = F (x,y)   where   (g,x,y) = gcd' a b 

@@ -8,6 +8,7 @@ module Algebra.Structures.GCDDomain    ( GCDDomain(gcd')   , propGCD, propGCDDomain+  , ggcd   ) where  import Test.QuickCheck@@ -16,6 +17,7 @@ import Algebra.Structures.BezoutDomain import Algebra.Ideal +-- infix 4 ~~  ------------------------------------------------------------------------------- -- | GCD domains@@ -41,10 +43,27 @@                          then propIntegralDomain a b c                          else whenFail (print "propGCD") False --- This can be used to compute gcd of a list of non-zero elements--- genGCD :: ?--- genGCD = ?+fst3 :: (a,b,c) -> a+fst3 (x,_,_) = x +-- Generalized greatest common divisor, computes the gcd of a list of elements.+ggcd :: GCDDomain a => [a] -> a+ggcd []     = error "ggcd: Can't compute ggcd of the empty list"+ggcd (x:xs) = foldr (\a b -> fst3 (gcd' a b)) x xs+ instance BezoutDomain a => GCDDomain a where   gcd' a b = (g,x,y)-    where (Id [g],_,[x,y]) = toPrincipal (Id [a,b])+   where (Id [g],_,[x,y]) = toPrincipal (Id [a,b])++{-+class IntegralDomain a => DecidableUnits a where+  unit :: a -> Maybe a -- Just x = the inverse++-- Divisibility is decidable if A is a gcd domain with decidable units+divides :: (GCDDomain a, DecidableUnits a) => a -> a -> Bool+divides a b = unit g+  where (g,_,_) = gcd' a b++(~~) :: (GCDDomain a, DecidableUnits a) => a -> a -> Bool+(~~) = divides+-}

@@ -1,24 +1,26 @@ {-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, FlexibleInstances #-}+-- | R-modules.  module Algebra.Structures.Module -  ( Module(..), (<*)+  ( Module((*>)), (<*)   , propScalarMul, propScalarAdd, propScalarAssoc, propModule   ) where  import Algebra.Structures.Group as G-import Algebra.Structures.CommutativeRing as R -- hiding ((<*),(*>))+import Algebra.Structures.CommutativeRing as R import Algebra.Z import Algebra.Zn  import Test.QuickCheck  infixl 7 *>-infixl 7 <*+infixr 7 <*  -- Consider only the commutative case, it would be possible to implement left -- and right modules instead. --- A module over a commutative ring r.+-- | Module over a commutative ring r. class (CommutativeRing r, AbelianGroup m) => Module r m where+  -- | Scalar multiplication.   (*>) :: r -> m -> m  propScalarMul :: (Module r m, Eq m) => r -> m -> m -> Bool@@ -38,11 +40,11 @@     (_,False,_)      -> whenFail (print "propScalarAdd") False     (_,_,False)      -> whenFail (print "propScalarAssoc") False --- Since the ring is commutative we can turn this around.+-- | Since the ring is commutative we can turn the scalar multiplication around. (<*) :: Module r m => m -> r -> m (<*) = flip (*>) --- Z-module+-- | Z-module structure. instance AbelianGroup m => Module Z m where   n *> x | n > 0  = sumGroup (replicate (fromInteger n) x)          | n == 0 = G.zero

@@ -25,7 +25,12 @@   ---------------------------------------------------------------------------------- | Prufer domain+-- | Given a and b it computes u, v and we such that:+-- +--  (1) au = bv+--+--  (2) b(1-u) = aw+-- class IntegralDomain a => PruferDomain a where   --         a    b     u v w   calcUVW :: a -> a -> (a,a,a)@@ -59,7 +64,7 @@   ----------------------------------------------------------------------------------- | Bezout domain -> Prüfer domain+-- | Bezout domain -> Prufer domain -- {- Prufer: forall a b exists u v w t.  u+t = 1 &  ua = vb & wa = tb@@ -179,7 +184,7 @@  -- | Compute the intersection of I and J by: --       ---       (I ∩ J)(I + J) = IJ  => (I ∩ J)(I + J)(I + J)' = IJ(I + J)'+--       (I \\cap J)(I + J) = IJ  => (I \\cap J)(I + J)(I + J)' = IJ(I + J)' -- intersectionPDWitness :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]]) intersectionPDWitness (Id is) (Id js) = (int,wis,wjs)

@@ -6,6 +6,7 @@   , propRing   , (<->), (<^>) -- , (*>), (<*)   , sumRing, productRing+--  , (~~)   ) where  import Test.QuickCheck@@ -17,7 +18,7 @@ -- infixl 7 <* infixl 6 <+> infixl 6 <->-+--infix  4 ~~  ------------------------------------------------------------------------------- -- | Definition of rings.@@ -111,6 +112,10 @@ x <^> y = if y < 0               then error "<^>: Input should be positive"              else x <*> x <^> (y-1)++-- | Check if a == b or -a == b or a == -b or -a == -b+-- (~~) :: (Ring a, Eq a) => a -> a -> Bool+-- x ~~ y = x == y || neg x == y || x == neg y || neg x == neg y  {- -- | Multiply from left with an integer; n *> x means x + x + ... + x, n times.

@@ -1,11 +1,14 @@ {-# LANGUAGE ScopedTypeVariables, FlexibleContexts #-} -- | Univariate polynomials parametrised by the variable name.-module Algebra.UPoly +module Algebra.UPoly   ( UPoly(..)   , deg   , Qx, x   , toUPoly, monomial   , lt, deriv+  , cont, isPrimitive+  , toPrimitive, propToPrimitive, gaussLemma+  , gcdUPolyWitness   , sqfr, sqfrDec   ) where @@ -13,17 +16,22 @@ import Test.QuickCheck import Control.Monad (liftM) -import Algebra.TypeChar.Char hiding (Q)+import Algebra.TypeChar.Char hiding (Z,Q) import Algebra.Structures.Field-import Algebra.Structures.BezoutDomain+import Algebra.Structures.FieldOfFractions import Algebra.Structures.EuclideanDomain+import Algebra.Structures.ExplicitUnits+import Algebra.Structures.BezoutDomain+import Algebra.Structures.GCDDomain import Algebra.Structures.PruferDomain-import Algebra.Structures.StronglyDiscrete+-- import Algebra.Structures.StronglyDiscrete+ import Algebra.Ideal+import Algebra.Z import Algebra.Q --- | Polynomials over a commutative ring, indexed by a phantom type x that --- denote the name of the variable that the polynomial is over. For example +-- | Polynomials over a commutative ring, indexed by a phantom type x that+-- denote the name of the variable that the polynomial is over. For example -- UPoly Q X_ is Q[x] and UPoly Q T_ is Q[t]. newtype CommutativeRing r => UPoly r x = UP [r]   deriving (Eq,Ord)@@ -44,7 +52,11 @@ toUPoly :: (CommutativeRing r, Eq r) => [r] -> UPoly r x toUPoly = UP . reverse . dropWhile (==zero) . reverse --- | Take an element of the ring and the degree of the desired monomial, for +isZeroPoly :: CommutativeRing r => UPoly r x -> Bool+isZeroPoly (UP []) = True+isZeroPoly _       = False++-- | Take an element of the ring and the degree of the desired monomial, for -- example: monomial 3 7 = 3x^7 monomial :: CommutativeRing r => r -> Integer -> UPoly r x monomial a i = UP $ replicate (fromInteger i) zero ++ [a]@@ -67,7 +79,7 @@  instance (CommutativeRing r, Eq r, Show r, Show x) => Show (UPoly r x) where   show (UP []) = "0"-  show (UP ps) = init $ fixSign $ concat +  show (UP ps) = init $ fixSign $ concat                   [ show' (show (undefined :: x)) p n                   | (p,n) <- zip ps [0..]                   , p /= zero ]@@ -75,10 +87,10 @@     show' :: (CommutativeRing r, Show r) => String -> r -> Integer -> String     show' x p 0 = show p ++ "+"     show' x p 1 = if p == one then x ++ "+" else show p ++ x ++ "+"-    show' x p n = if p == one  -                     then x ++ "^" ++ show n ++ "+" +    show' x p n = if p == one+                     then x ++ "^" ++ show n ++ "+"                      else show p ++ x ++ "^" ++ show n ++ "+"-    +     fixSign []  = []     fixSign [x] = [x]     fixSign ('+':'-':xs) = '-' : fixSign xs@@ -89,9 +101,9 @@  -- Addition of polynomials. addUP :: (CommutativeRing r, Eq r) => UPoly r x -> UPoly r x -> UPoly r x-addUP (UP ps) (UP qs) | length ps >= length qs = add' ps qs +addUP (UP ps) (UP qs) | length ps >= length qs = add' ps qs                       | otherwise              = add' qs ps-  where add' a b = toUPoly $ zipWith (<+>) a b ++ drop (length b) a +  where add' a b = toUPoly $ zipWith (<+>) a b ++ drop (length b) a  -- Multiplication of polynomials. mulUP :: (CommutativeRing r, Eq r) => UPoly r x -> UPoly r x -> UPoly r x@@ -99,7 +111,7 @@   where   m ps qs r | r > length ps + length qs - 2 = []             | otherwise = c r 0 (length ps-1) (length qs-1) : m ps qs (r+1)-  +   c (-1) _ _ _ = zero   c r k m n | r > m || k > n = c (r-1) (k+1) m n             | otherwise      = ps !! r <*> qs !! k <+> c (r-1) (k+1) m n@@ -115,7 +127,7 @@   (+)    = (<+>)   (-)    = (<->)   (*)    = (<*>)-  abs    = fromInteger . d+  abs    = fromInteger . norm   signum = undefined -- Is it possible to define this?   fromInteger x = UP [fromInteger x] @@ -124,31 +136,151 @@  -- Polynomial rings are Euclidean. instance (Field k, Eq k) => EuclideanDomain (UPoly k x) where-  d (UP ps)             = fromIntegral (length ps) - 1+  norm (UP ps)             = fromIntegral (length ps) - 1   quotientRemainder f g = qr zero f     where     -- This is the division algorithm in k[x]. Page 39 in Cox.-    qr q r | d g <= d r = qr (q <+> monomial (lt r </> lt g) (d r - d g))-                            (r <-> monomial (lt r </> lt g) (d r - d g) <*> g)+    qr q r | norm g <= norm r = +              qr (q <+> monomial (lt r </> lt g) (norm r - norm g))+                 (r <-> monomial (lt r </> lt g) (norm r - norm g) <*> g)            | otherwise = (q,r)  instance (Field k, Eq k) => PruferDomain (UPoly k x) where   calcUVW = calcUVW_B +instance (ExplicitUnits a, Eq a) => ExplicitUnits (UPoly a x) where+  unit (UP [a]) = case unit a of+    Just a' -> Just (UP [a'])+    Nothing -> Nothing+  unit _        = Nothing+ -- Now that we know that the polynomial ring k[x] is a Bezout domain it is -- possible to implement membership in an ideal of k[x]. f is a member of the -- ideal <f1,...,fn> if the rest is zero when dividing f with the principal -- ideal <h>. -- instance (Field k, Eq k, Show x) => StronglyDiscrete (UPoly k x) where---  member p ps = modulo p h == zero +--  member p ps = modulo p h == zero --    where Id [h] = (\(a,_,_) -> a) $ toPrincipal ps ++-------------------------------------------------------------------------------+-- Proof that if A is a GCD domain then A[x] is a GCD domain following+-- section 4.4 in A course in constructive algebra.++-- Some test polynomials:++test1 :: UPoly Z X_+test1 = toUPoly [1,2,3]++test2 :: UPoly Z X_+test2 = toUPoly [2,4,6,8,10]++test3 :: UPoly Q X_+test3 = toUPoly [inv 2, inv 3, inv 4]+++-- | Compute the content of a polynomial, i.e. the gcd of the coefficients.+cont :: (GCDDomain a, Eq a) => UPoly a x -> a+cont (UP xs) = case filter (/= zero) xs of+  []  -> error "cont: Can't compute the content of the zero polynomial"+  xs' -> ggcd xs'++-- *Algebra.UPoly> cont test1+-- 1+-- *Algebra.UPoly> cont test2+-- 2+++-- | If all coefficients are relatively prime then the polynomial is primitive.+isPrimitive :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly a x -> Bool+isPrimitive = isUnit . cont++-- *Algebra.UPoly> isPrimitive test1+-- True+-- *Algebra.UPoly> isPrimitive test2+-- False+++-- | Lemma 4.2: Given a polynomial p in K[x] where K=Quot(A) we can find c in K+-- and q primitive in A[x] such that p = cq.+toPrimitive :: (GCDDomain a, Eq a)+            => UPoly (FieldOfFractions a) x+            -> (FieldOfFractions a, UPoly a x)+toPrimitive p@(UP xs) = (c,q)+  where+  c0' = toFieldOfFractions $ productRing $ map denominator xs+  c0  = inv c0'+  g0  = map (fromFieldOfFractions . (c0' <*>)) xs+  cg0 = toFieldOfFractions $ cont $ UP g0+  c   = c0 <*> cg0+  q   = toUPoly $ map (\x -> fromFieldOfFractions (toFieldOfFractions x </> cg0)) g0++-- *Algebra.UPoly> toPrimitive test3+-- (1/12,6+4x+3x^2)++-- Specification of toPrimitive.+propToPrimitive :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly (FieldOfFractions a) x -> Property+propToPrimitive p =+  not (isZeroPoly p) ==>+    p == toUPoly (map (\x -> c <*> toFieldOfFractions x) q) && isPrimitive (UP q)+  where (c,UP q) = toPrimitive p++-- *Algebra.UPoly> quickCheck (propToPrimitive :: UPoly Q X_ -> Property)+-- +++ OK, passed 100 tests.+++-- | Gauss lemma says that if p and q are polynomials over a GCD domain then+-- cont(pq) = cont(p) * cont(q).+gaussLemma :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly a x -> UPoly a x -> Property+gaussLemma p q =+  not (isZeroPoly p) && not (isZeroPoly q) ==>+    cont (p <*> q) ~= cont p <*> cont q++-- *Algebra.UPoly> quickCheck (gaussLemma :: UPoly Z X_ -> UPoly Z X_ -> Property)+-- +++ OK, passed 100 tests.++liftUPoly :: (GCDDomain a, Eq a) => UPoly a x -> UPoly (FieldOfFractions a) x+liftUPoly (UP xs) = toUPoly $ map toFieldOfFractions xs++-- | Proof that if A is a GCD domain then A[x] also is a GCD domain. This also+-- computes witnesses that the computed GCD divides the given polynomials.+gcdUPolyWitness :: (GCDDomain a, Eq a)+                => UPoly a x+                -> UPoly a x+                -> (UPoly a x, UPoly a x, UPoly a x)+gcdUPolyWitness p q = (constUPoly d <*> h, constUPoly x <*> a, constUPoly y <*> b)+  where+  (h',a',b') = gcd' (liftUPoly p) (liftUPoly q)++  (_,h) = toPrimitive h'+  (_,a) = toPrimitive a'+  (_,b) = toPrimitive b'++  (d,x,y) = gcd' (cont p) (cont q)++test4, test5 :: UPoly Z X_+test4 = toUPoly [6,7,1]+test5 = toUPoly [-6,-5,1]++-- *Algebra.UPoly> gcdUPolyWitness test4 test5+-- (1+x,6+x,-6+x)+-- *Algebra.UPoly> gcdUPolyWitness test4 test4+-- (6+7x+x^2,1,1)+-- *Algebra.UPoly> gcdUPolyWitness test1 test2+-- (1,1+2x+3x^2,2+4x+6x^2+8x^3+10x^4)++-- This does not work:+-- instance (GCDDomain a, Eq a) => GCDDomain (UPoly a x) where+--   gcd' = gcdUPolyWitness+++------------------------------------------------------------------------------- -- Square free decomposition. -- Teo Mora; Solving Polynomial Equations Systems I. pg 69-70 -- Works only for Char 0 --TODO: Add check for char --square free associate of f--- | Square free decomposition of a polynomial. +-- | Square free decomposition of a polynomial. sqfr :: (Num k, Field k) => UPoly k x -> UPoly k x sqfr f = f `quotient` euclidAlg f f'   where f' = deriv f@@ -156,24 +288,24 @@ -- | Distinct power factorization, aka square free decomposition sqfrDec :: (Num k, Field k) => UPoly k x -> [UPoly k x] sqfrDec f = help p q-  where +  where   p = euclidAlg f (deriv f)-  q = f `quotient` p -  -  help p q | d q < 1    = []+  q = f `quotient` p++  help p q | norm q < 1    = []            | otherwise  = t : help (p `quotient` s) s-    where +    where     s = euclidAlg p q     t = q `quotient` s  -- | Pseudo-division of polynomials.--- +-- -- Given s(x) and p(x) compute c, q(x) and r(x) such that:---   +-- --   cs(x) = p(x)q(x)+r(x), deg r < deg p.-pseudoDivide :: (CommutativeRing a, Eq a) +pseudoDivide :: (CommutativeRing a, Eq a)              => UPoly a x -> UPoly a x -> (a, UPoly a x, UPoly a x)-pseudoDivide s p +pseudoDivide s p   | m < n     = (one,zero,s)   | otherwise = pD (a' <*> s' <-> b' <*> xmn <*> p') 1 (b' <*> xmn) s2   where@@ -189,7 +321,7 @@   p'  = p <-> monomial a n   s2  = a' <*> s' <-> b' <*> xmn <-> p' -  pD s k out1 out2 +  pD s k out1 out2     | deg s < n = (a <^> k,out1,out2)     | otherwise = pD s3 (k+1) (b2xm2n <+> a' <*> out1) s3     where@@ -200,6 +332,6 @@     s3 = (a' <*> s) <-> (b2xm2n <*> p)  -propPD :: Qx -> Qx -> Property+propPD :: (CommutativeRing a, Eq a) => UPoly a x -> UPoly a x -> Property propPD s p = deg s > 1 && deg p > 1 ==> constUPoly c <*> s == p <*> q <+> r-  where (c,q,r) = pseudoDivide s p +  where (c,q,r) = pseudoDivide s p

@@ -1,15 +1,17 @@ {-# LANGUAGE TypeSynonymInstances #-} module Algebra.Z    ( Z-  , Ring(..)+  , IntegralDomain(..)   ) where  import Test.QuickCheck  import Algebra.Structures.IntegralDomain import Algebra.Structures.EuclideanDomain+import Algebra.Structures.ExplicitUnits import Algebra.Structures.StronglyDiscrete import Algebra.Structures.BezoutDomain+import Algebra.Structures.GCDDomain import Algebra.Structures.PruferDomain import Algebra.Structures.Coherent import Algebra.Ideal@@ -34,17 +36,28 @@ propIntegralDomainZ :: Z -> Z -> Z -> Property propIntegralDomainZ = propIntegralDomain +instance ExplicitUnits Z where+  unit 1    = Just 1+  unit (-1) = Just (-1)+  unit _    = Nothing+ instance EuclideanDomain Z where-  d = abs+  norm = abs   quotientRemainder = quotRem  propEuclideanDomainZ :: Z -> Z -> Z -> Property propEuclideanDomainZ = propEuclideanDomain  -- Euclidean domain -> Bezout domain-propBezoutDomainZ :: Ideal Z -> Z -> Z -> Z -> Property+propBezoutDomainZ :: Z -> Z -> Property propBezoutDomainZ = propBezoutDomain +propToPrincipalZ :: Ideal Z -> Bool+propToPrincipalZ = propToPrincipal++propIsSameIdealZ :: Ideal Z -> Bool+propIsSameIdealZ = propIsSameIdeal+ -- Bezout domain -> Strongly discrete propStronglyDiscreteZ :: Z -> Ideal Z -> Bool propStronglyDiscreteZ = propStronglyDiscrete@@ -65,6 +78,11 @@ -- Not working perfectly... propSolveGeneralZ :: Matrix Z -> Vector Z -> Property propSolveGeneralZ = propSolveGeneral+++-- GCD Domain+propGCDDomainZ :: Z -> Z -> Z -> Property+propGCDDomainZ = propGCDDomain  -- PLM propPLMZ :: Ideal Z -> Bool

@@ -2,7 +2,7 @@              FunctionalDependencies, FlexibleContexts,    UndecidableInstances,              FlexibleInstances #-} -- | Integers modulo n parametrised by the n. This also has type-level primality--- testing used for instantiating integral domain and field type classes. The +-- testing used for instantiating integral domain and field type classes. The -- primality testing is very slow, but it seem to be working fine for relatively -- small numbers. module Algebra.Zn (Zn(..), Z3) where@@ -28,12 +28,12 @@   signum (Zn x) = Zn $ signum x   negate (Zn x) = Zn $ negate x `mod` toNum (undefined :: n)   fromInteger x = Zn $ fromInteger $ x `mod` toNum (undefined :: n)-  + instance Nat n => Arbitrary (Zn n) where   arbitrary = liftM Zn (choose (0,toNum (undefined :: n) - 1))- + instance Nat n => Ring (Zn n) where-  (<+>) = (+) +  (<+>) = (+)   zero  = Zn 0   one   = Zn 1   neg   = negate@@ -61,13 +61,13 @@  type Z17 = Zn D17 -test2 :: Z17-test2 = inv 13+-- test2 :: Z17+-- test2 = inv 13  -- Test that all elements of Z17 get correct inverses-test3 :: Prelude.Bool-test3 = all (==1) [ inv x * x | x <- xs ]-  where xs :: [Z17] = map fromInteger [1..16]+-- test3 :: Prelude.Bool+-- test3 = all (==1) [ inv x * x | x <- xs ]+--   where xs :: [Z17] = map fromInteger [1..16]  ----------------------------------------------------------------------- -- Lots of crazy type-level stuff:@@ -78,7 +78,7 @@ class Sqrt' x y r sqrt | x y r -> sqrt instance Sub y D2 y' => Sqrt' x y LT y' instance Pred y y'   => Sqrt' x y EQ y'-instance (ExpBase y D2 square, Succ y y', Trich x square r, +instance (ExpBase y D2 square, Succ y y', Trich x square r,           Sqrt' x y' r sqrt) => Sqrt' x y GT sqrt  sqrt :: Sqrt x sqrt => x -> sqrt@@ -89,7 +89,7 @@  class Data.TypeLevel.Bool b => Prime' x y r b | x y r -> b instance Prime' x D1 EQ True-instance (Pred y z, Trich z D1 r1, Mod x y rest, IsZero rest b1, +instance (Pred y z, Trich z D1 r1, Mod x y rest, IsZero rest b1,           Not b1 b', Prime' x z r1 b2, And b' b2 b3) => Prime' x y GT b3  prime :: Prime x b => x -> b