diff --git a/docs/NOTES b/docs/NOTES --- a/docs/NOTES +++ b/docs/NOTES @@ -83,6 +83,17 @@ Juergen Bokowski DMV-Nachrichten 2004/3 + Gerhard Navratil + Partial derivatives, Haskell-Cafe 2006-05-08 + + Paul Johnson + interval arithmetic + http://sourceforge.net/projects/ranged-sets + + David Amos + algebra + http://polyomino.f2s.com/ + * RealFloat Defines the properties of a Floating type, thus should be named 'Floating'. diff --git a/numeric-prelude.cabal b/numeric-prelude.cabal --- a/numeric-prelude.cabal +++ b/numeric-prelude.cabal @@ -1,5 +1,5 @@ Name: numeric-prelude -Version: 0.1.1 +Version: 0.1.2 License: GPL License-File: LICENSE Author: Dylan Thurston , Henning Thielemann , Mikael Johansson @@ -185,6 +185,7 @@ MathObj.DiscreteMap MathObj.LaurentPolynomial MathObj.Matrix + MathObj.Monoid MathObj.PartialFraction MathObj.Permutation MathObj.Permutation.CycleList @@ -204,6 +205,7 @@ Number.DimensionTerm.SI Number.FixedPoint Number.FixedPoint.Check + Number.GaloisField2p32m5 Number.NonNegative Number.NonNegativeChunky Number.PartiallyTranscendental @@ -246,11 +248,14 @@ GHC-Options: -Wall Other-modules: Test.NumericPrelude.Utility + Test.Number.GaloisField2p32m5 Test.MathObj.PartialFraction + Test.MathObj.Matrix Test.MathObj.Polynomial Test.MathObj.PowerSeries Test.MathObj.Gaussian.Variance Test.MathObj.Gaussian.Bell + Test.MathObj.Gaussian.Polynomial Main-Is: Test/Run.hs If flag(buildTests) Build-Depends: HUnit >=1 && <2 diff --git a/src/Algebra/Monoid.hs b/src/Algebra/Monoid.hs --- a/src/Algebra/Monoid.hs +++ b/src/Algebra/Monoid.hs @@ -5,16 +5,59 @@ Stability : provisional Portability : -Abstract concept of a Monoid. Will be used in order to generate -type classes for generic algebras. An algebra is a vector space -that also is a monoid. +Abstract concept of a Monoid. +Will be used in order to generate type classes for generic algebras. +An algebra is a vector space that also is a monoid. +Should we use the Monoid class from base library +despite its unfortunate method name @mappend@? -} module Algebra.Monoid where -{- | We expect a monoid to adher to associativity and the identity -behaving decently. Nothing more, really. -} +import qualified Algebra.Additive as Additive +import qualified Algebra.Ring as Ring +import Data.Monoid as Mn + +{- | +We expect a monoid to adher to associativity and +the identity behaving decently. +Nothing more, really. +-} class C a where idt :: a (<*>) :: a -> a -> a + +instance C Mn.All where + idt = mempty + (<*>) = mappend + +instance C Any where + idt = mempty + (<*>) = mappend + +instance C a => C (Dual a) where + idt = Mn.Dual idt + (Mn.Dual x) <*> (Mn.Dual y) = Mn.Dual (y <*> x) + +instance C (Endo a) where + idt = mempty + (<*>) = mappend + +instance C (First a) where + idt = mempty + (<*>) = mappend + +instance C (Last a) where + idt = mempty + (<*>) = mappend + + +instance Ring.C a => C (Product a) where + idt = Mn.Product Ring.one + (Mn.Product x) <*> (Mn.Product y) = Mn.Product (x Ring.* y) + +instance Additive.C a => C (Sum a) where + idt = Mn.Sum Additive.zero + (Mn.Sum x) <*> (Mn.Sum y) = Mn.Sum (x Additive.+ y) + diff --git a/src/Algebra/PrincipalIdealDomain.hs b/src/Algebra/PrincipalIdealDomain.hs --- a/src/Algebra/PrincipalIdealDomain.hs +++ b/src/Algebra/PrincipalIdealDomain.hs @@ -55,6 +55,8 @@ import Control.Monad (foldM, liftM) import Data.List (mapAccumL, mapAccumR, unfoldr) +import Data.Int (Int, Int8, Int16, Int32, Int64, ) + import PreludeBase import Prelude (Integer, Int) import qualified Prelude as P @@ -137,6 +139,11 @@ 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 = @@ -288,6 +295,18 @@ 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 diff --git a/src/Algebra/Units.hs b/src/Algebra/Units.hs --- a/src/Algebra/Units.hs +++ b/src/Algebra/Units.hs @@ -32,6 +32,8 @@ import Algebra.Additive (negate) import Algebra.ZeroTestable (isZero) +import Data.Int (Int, Int8, Int16, Int32, Int64, ) + import PreludeBase import Prelude (Integer, Int) import qualified Prelude as P @@ -96,6 +98,30 @@ stdUnitInv = intStandardInverse instance C Integer where + isUnit = intQuery + stdAssociate = intAssociate + stdUnit = intStandard + stdUnitInv = intStandardInverse + +instance C Int8 where + isUnit = intQuery + stdAssociate = intAssociate + stdUnit = intStandard + stdUnitInv = intStandardInverse + +instance C Int16 where + isUnit = intQuery + stdAssociate = intAssociate + stdUnit = intStandard + stdUnitInv = intStandardInverse + +instance C Int32 where + isUnit = intQuery + stdAssociate = intAssociate + stdUnit = intStandard + stdUnitInv = intStandardInverse + +instance C Int64 where isUnit = intQuery stdAssociate = intAssociate stdUnit = intStandard diff --git a/src/MathObj/Matrix.hs b/src/MathObj/Matrix.hs --- a/src/MathObj/Matrix.hs +++ b/src/MathObj/Matrix.hs @@ -2,16 +2,45 @@ {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FlexibleInstances #-} {- | -Copyright : (c) Mikael Johansson 2006 -Maintainer : mik@math.uni-jena.de +Copyright : (c) Henning Thielemann 2009, Mikael Johansson 2006 +Maintainer : numericprelude@henning-thielemann.de Stability : provisional Portability : requires multi-parameter type classes Routines and abstractions for Matrices and basic linear algebra over fields or rings. + +We stick to simple Int indices. +Although advanced indices would be nice +e.g. for matrices with sub-matrices, +this is not easily implemented since arrays +do only support a lower and an upper bound +but no additional parameters. + +ToDo: + - Matrix inverse, determinant -} -module MathObj.Matrix where +module MathObj.Matrix ( + T, Dimension, + format, + transpose, + rows, + columns, + fromRows, + fromColumns, + fromList, + dimension, + numRows, + numColumns, + zipWith, + zero, + one, + diagonal, + scale, + random, + randomR, + ) where import qualified Algebra.Module as Module import qualified Algebra.Vector as Vector @@ -20,32 +49,34 @@ import Algebra.Module((*>), ) import Algebra.Ring((*), fromInteger, scalarProduct, ) -import Algebra.Additive((+), (-), zero, subtract, ) +import Algebra.Additive((+), (-), subtract, ) -import Data.Array (Array, listArray, elems, bounds, (!), ixmap, range, ) +import qualified System.Random as Rnd +import Data.Array (Array, array, listArray, accumArray, elems, bounds, (!), ixmap, range, ) import qualified Data.List as List import Control.Monad (liftM2, ) import Control.Exception (assert, ) -import Data.Tuple.HT (swap, ) +import Data.Tuple.HT (swap, mapFst, ) import Data.List.HT (outerProduct, ) -import NumericPrelude (Integer, ) + +import NumericPrelude (Int, ) import PreludeBase hiding (zipWith, ) + {- | -A matrix is a twodimensional array of ring elements, indexed by integers. +A matrix is a twodimensional array, indexed by integers. -} +data T a = + Cons (Array (Dimension, Dimension) a) + deriving (Eq,Ord,Read) -data {-(Ring.C a) =>-} - T a = Cons (Array (Integer, Integer) a) deriving (Eq,Ord,Read) +type Dimension = Int {- | -Transposition of matrices is just transposition in the sense of -Data.List. +Transposition of matrices is just transposition in the sense of Data.List. -} - - transpose :: T a -> T a transpose (Cons m) = let (lower,upper) = bounds m @@ -59,31 +90,63 @@ columns :: T a -> [[a]] columns (Cons m) = let ((lr,lc), (ur,uc)) = bounds m - in outerProduct (curry(m!)) (range (lc,uc)) (range (lr,ur)) + in outerProduct (flip(curry(m!))) (range (lc,uc)) (range (lr,ur)) -fromList :: Integer -> Integer -> [a] -> T a -fromList m n xs = Cons (listArray ((1,1),(m,n)) xs) +fromRows :: Dimension -> Dimension -> [[a]] -> T a +fromRows m n = + Cons . + array (indexBounds m n) . + concat . + List.zipWith (\r -> map (\(c,x) -> ((r,c),x))) allIndices . + map (zip allIndices) -instance (Ring.C a, Show a) => Show (T a) where - show m = List.unlines $ map (\r -> "(" ++ r ++ ")") - $ map (unwords . map show) $ rows m +fromColumns :: Dimension -> Dimension -> [[a]] -> T a +fromColumns m n = + Cons . + array (indexBounds m n) . + concat . + List.zipWith (\r -> map (\(c,x) -> ((c,r),x))) allIndices . + map (zip allIndices) +fromList :: Dimension -> Dimension -> [a] -> T a +fromList m n xs = Cons (listArray (indexBounds m n) xs) -dimension :: T a -> (Integer,Integer) +appPrec :: Int +appPrec = 10 + +instance (Show a) => Show (T a) where + showsPrec p m = + showParen (p >= appPrec) + (showString "Matrix.fromRows " . showsPrec appPrec (rows m)) + +format :: (Show a) => T a -> String +format m = formatS m "" + +formatS :: (Show a) => T a -> ShowS +formatS = + concatS . + map (\r -> showString "(" . concatS r . showString ")\n") . + map (List.intersperse (' ':) . map (showsPrec 11)) . + rows + +concatS :: [ShowS] -> ShowS +concatS = flip (foldr ($)) + +dimension :: T a -> (Dimension,Dimension) dimension (Cons m) = uncurry subtract (bounds m) + (1,1) -numRows :: T a -> Integer +numRows :: T a -> Dimension numRows = fst . dimension -numColumns :: T a -> Integer +numColumns :: T a -> Dimension numColumns = snd . dimension -- These implementations may benefit from a better exception than -- just assertions to validate dimensionalities instance (Additive.C a) => Additive.C (T a) where - (+) = zipWith (+) - (-) = zipWith (-) - zero = zeroMatrix 1 1 + (+) = zipWith (+) + (-) = zipWith (-) + zero = zero 1 1 zipWith :: (a -> b -> c) -> T a -> T b -> T c zipWith op mM@(Cons m) nM@(Cons n) = @@ -93,30 +156,71 @@ in assert (d == dimension nM) $ uncurry fromList d (List.zipWith op em en) -zeroMatrix :: (Additive.C a) => Integer -> Integer -> T a -zeroMatrix m n = fromList m n $ - List.repeat zero +zero :: (Additive.C a) => Dimension -> Dimension -> T a +zero m n = + fromList m n $ + List.repeat Additive.zero -- List.replicate (fromInteger (m*n)) zero +one :: (Ring.C a) => Dimension -> T a +one n = + Cons $ + accumArray (flip const) Additive.zero + (indexBounds n n) + (map (\i -> ((i,i), Ring.one)) (indexRange n)) + +diagonal :: (Additive.C a) => [a] -> T a +diagonal xs = + let n = List.length xs + in Cons $ + accumArray (flip const) Additive.zero + (indexBounds n n) + (zip (map (\i -> (i,i)) allIndices) xs) + +scale :: (Ring.C a) => a -> T a -> T a +scale s = Vector.functorScale s + instance (Ring.C a) => Ring.C (T a) where - mM * nM = assert (numRows mM == numColumns nM) $ - fromList (numColumns mM) (numRows nM) - (liftM2 scalarProduct (rows mM) (columns nM)) - fromInteger n = fromList 1 1 [fromInteger n] + mM * nM = + assert (numColumns mM == numRows nM) $ + fromList (numRows mM) (numColumns nM) + (liftM2 scalarProduct (rows mM) (columns nM)) + fromInteger n = fromList 1 1 [fromInteger n] instance Functor T where fmap f (Cons m) = Cons (fmap f m) instance Vector.C T where - zero = zero + zero = Additive.zero (<+>) = (+) - (*>) = Vector.functorScale + (*>) = scale instance Module.C a b => Module.C a (T b) where x *> m = fmap (x*>) m -{- | -What more do we need from our matrix class? We have addition, + +random :: (Rnd.RandomGen g, Rnd.Random a) => + Dimension -> Dimension -> g -> (T a, g) +random = + randomAux Rnd.random + +randomR :: (Rnd.RandomGen g, Rnd.Random a) => + Dimension -> Dimension -> (a,a) -> g -> (T a, g) +randomR m n rng = + randomAux (Rnd.randomR rng) m n + +{- +could be made nicer with the State monad, +but I like to keep dependencies minimal +-} +randomAux :: (Rnd.RandomGen g, Rnd.Random a) => + (g -> (a,g)) -> Dimension -> Dimension -> g -> (T a, g) +randomAux rnd m n g0 = + mapFst (fromList m n) $ swap $ + List.mapAccumL (\g _i -> swap $ rnd g) g0 (indexRange (m*n)) + +{- +What more do we need from our matrix type? We have addition, subtraction and multiplication, and thus composition of generic free-module-maps. We're going to want to solve linear equations with or without fields underneath, so we're going to want an implementation @@ -125,6 +229,7 @@ with the Gaussian algorithm or some other goodish method. -} +{- {- | We'll want generic linear equation solving, returning one solution, any solution really, or nothing. Basically, this is asking for the @@ -144,7 +249,23 @@ (numRows a == numRows y && -- they match numColumns y == 1) -- and y is a column vector Nothing +-} {- Cf. /usr/lib/hugs/demos/Matrix.hs -} + + +-- these functions control whether we use 0 or 1 based indices + +indexRange :: Dimension -> [Dimension] +indexRange n = [0..(n-1)] + +indexBounds :: + Dimension -> Dimension -> + ((Dimension,Dimension), (Dimension,Dimension)) +indexBounds m n = + ((0,0), (m-1,n-1)) + +allIndices :: [Dimension] +allIndices = [0..] diff --git a/src/MathObj/Monoid.hs b/src/MathObj/Monoid.hs new file mode 100644 --- /dev/null +++ b/src/MathObj/Monoid.hs @@ -0,0 +1,56 @@ +{-# LANGUAGE NoImplicitPrelude #-} +module MathObj.Monoid where + +import qualified Algebra.PrincipalIdealDomain as PID + +import Algebra.PrincipalIdealDomain (gcd, lcm, ) +import Algebra.Additive (zero, ) +import Algebra.Monoid (C, idt, (<*>), ) + +import PreludeBase + +{- | +It is only a monoid for non-negative numbers. + +> idt <*> GCD (-2) = GCD 2 + +Thus, use this Monoid only for non-negative numbers! +-} +newtype GCD a = GCD {runGCD :: a} + deriving (Show, Eq) + +instance PID.C a => C (GCD a) where + idt = GCD zero + (GCD x) <*> (GCD y) = GCD (gcd x y) + + +newtype LCM a = LCM {runLCM :: a} + deriving (Show, Eq) + +instance PID.C a => C (LCM a) where + idt = LCM zero + (LCM x) <*> (LCM y) = LCM (lcm x y) + + +{- | +@Nothing@ is the largest element. +-} +newtype Min a = Min {runMin :: Maybe a} + deriving (Show, Eq) + +instance Ord a => C (Min a) where + idt = Min Nothing + (Min x) <*> (Min y) = Min $ + maybe y (\x' -> maybe x (Just . min x') y) x + + +{- | +@Nothing@ is the smallest element. +-} +newtype Max a = Max {runMax :: Maybe a} + deriving (Show, Eq) + +instance Ord a => C (Max a) where + idt = Max Nothing + (Max x) <*> (Max y) = Max $ + maybe y (\x' -> maybe x (Just . max x') y) x diff --git a/src/Number/GaloisField2p32m5.hs b/src/Number/GaloisField2p32m5.hs new file mode 100644 --- /dev/null +++ b/src/Number/GaloisField2p32m5.hs @@ -0,0 +1,93 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{- | +This number type is intended for tests of functions over fields, +where the field elements need constant space. +This way we can provide a Storable instance. +For 'Rational' this would not be possible. + +However, be aware that sums of non-zero elements may yield zero. +Thus division is not always safe, where it is for rational numbers. +-} +module Number.GaloisField2p32m5 where + +import qualified Number.ResidueClass as RC +import qualified Algebra.Module as Module +import qualified Algebra.Field as Field +import qualified Algebra.Ring as Ring +import qualified Algebra.Additive as Additive + +import Data.Int (Int64, ) +import Data.Word (Word32, Word64, ) + +import qualified Foreign.Storable.Newtype as SN +import qualified Foreign.Storable as St + +import Test.QuickCheck (Arbitrary(..), ) + +import PreludeBase +import NumericPrelude + + +newtype T = Cons {decons :: Word32} + deriving Eq + +{-# INLINE appPrec #-} +appPrec :: Int +appPrec = 10 + +instance Show T where + showsPrec p (Cons x) = + showsPrec p x +{- + showParen (p >= appPrec) + (showString "GF2p32m5.Cons " . shows x) +-} + +instance Arbitrary T where + arbitrary = fmap (Cons . fromInteger . flip mod base) arbitrary + coarbitrary = undefined + +instance St.Storable T where + sizeOf = SN.sizeOf decons + alignment = SN.alignment decons + peek = SN.peek Cons + poke = SN.poke decons + + +base :: Ring.C a => a +base = 2^32-5 + + +{-# INLINE lift2 #-} +lift2 :: (Word64 -> Word64 -> Word64) -> (T -> T -> T) +lift2 f (Cons x) (Cons y) = + Cons (fromIntegral (mod (f (fromIntegral x) (fromIntegral y)) base)) + +{-# INLINE lift2Integer #-} +lift2Integer :: (Int64 -> Int64 -> Int64) -> (T -> T -> T) +lift2Integer f (Cons x) (Cons y) = + Cons (fromIntegral (mod (f (fromIntegral x) (fromIntegral y)) base)) + + +instance Additive.C T where + zero = Cons zero + (+) = lift2 (+) +-- (-) = lift2 (-) + x-y = x + negate y + negate n@(Cons x) = + if x==0 + then n + else Cons (base-x) + +instance Ring.C T where + one = Cons one + (*) = lift2 (*) + fromInteger = + Cons . fromInteger . flip mod base + +instance Field.C T where + (/) = lift2Integer (RC.divide base) + +instance Module.C T T where + (*>) = (*) diff --git a/test/Test/MathObj/Gaussian/Polynomial.hs b/test/Test/MathObj/Gaussian/Polynomial.hs new file mode 100644 --- /dev/null +++ b/test/Test/MathObj/Gaussian/Polynomial.hs @@ -0,0 +1,144 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} +module Test.MathObj.Gaussian.Polynomial where + +import qualified MathObj.Gaussian.Polynomial as G +import qualified MathObj.Gaussian.Bell as B + +import qualified MathObj.Polynomial as Poly + +-- import qualified Algebra.Ring as Ring + +import qualified Algebra.Laws as Laws + +import qualified Number.Complex as Complex + +import Test.NumericPrelude.Utility (testUnit) +import Test.QuickCheck (Testable, quickCheck, (==>)) +import qualified Test.HUnit as HUnit + +import qualified Number.NonNegative as NonNeg +import Data.Function.HT (nest, ) +import Data.Tuple.HT (mapSnd, ) + +-- import Debug.Trace (trace, ) + +import PreludeBase as P +import NumericPrelude as NP + + +simple :: + (Testable t) => + (G.T Rational -> t) -> IO () +simple f = + quickCheck (\x -> f (x :: G.T Rational)) + +tests :: HUnit.Test +tests = + HUnit.TestLabel "polynomial" $ + HUnit.TestList $ + map testUnit $ + testList + +testList :: [(String, IO ())] +testList = +{- + ("convolution, dirac", + simple $ Laws.identity (+) zero) : +-} + ("convolution, commutative", + simple $ Laws.commutative G.convolve) : +-- simple $ \x -> Laws.commutative G.convolve (trace (show x) x)) : + ("convolution, associative", + simple $ Laws.associative G.convolve) : + ("multiplication, one", + simple $ Laws.identity G.multiply G.constant) : + ("multiplication, commutative", + simple $ Laws.commutative G.multiply) : + ("multiplication, associative", + simple $ Laws.associative G.multiply) : + ("convolution, multplication, fourier", + simple $ \x y -> + G.fourier (G.convolve x y) + == G.multiply (G.fourier x) (G.fourier y)) : + ("fourier reverse", + simple $ \x -> nest 2 G.fourier x == G.reverse x) : + ("reverse identity", + simple $ \x -> nest 2 G.reverse x == x) : + ("fourier eigenfunction differential", + quickCheck $ \m -> + m <= 15 ==> + let n = NonNeg.toNumber m + x = G.eigenfunctionDifferential n :: G.T Rational + k = Complex.conjugate Complex.imaginaryUnit ^ fromIntegral n + in G.fourier x == G.scaleComplex k x) : + ("fourier eigenfunction iterative", + quickCheck $ \m -> + m <= 15 ==> + let n = NonNeg.toNumber m + x = G.eigenfunctionIterative n :: G.T Rational + k = Complex.conjugate Complex.imaginaryUnit ^ fromIntegral n + in G.fourier x == G.scaleComplex k x) : +{- this does not hold, both functions compute different eigenbases + ("fourier eigenfunction diff vs. iterative", + quickCheck $ \n -> + n <= 15 ==> + G.eigenfunctionDifferential n == + (G.eigenfunctionIterative n :: G.T Rational)) : +-} + ("translate additive", + simple $ \x a b -> + G.translate a (G.translate b x) == G.translate (a+b) x) : + ("translateComplex additive", + simple $ \x a b -> + G.translateComplex a (G.translateComplex b x) == G.translateComplex (a+b) x) : + ("translateComplex real", + simple $ \x a -> + G.translateComplex (Complex.fromReal a) x == G.translate a x) : + ("modulate additive", + simple $ \x a b -> + G.modulate a (G.modulate b x) == G.modulate (a+b) x) : + ("commute translate modulate", + simple $ \x a b -> + G.modulate b (G.translate a x) + == G.turn (a*b) (G.translate a (G.modulate b x))) : + ("fourier translate", + simple $ \x a -> + G.fourier (G.translate a x) + == G.modulate a (G.fourier x)) : + ("dilate multiplicative", + simple $ \x a b -> a>0 && b>0 ==> + G.dilate a (G.dilate b x) == G.dilate (a*b) x) : + ("dilate by shrink", + simple $ \x a -> a>0 ==> + G.shrink a x == G.dilate (recip a) x) : + ("fourier dilate", + simple $ \x a -> a>0 ==> + G.fourier (G.dilate a x) == G.amplify a (G.shrink a (G.fourier x))) : + ("integrate differentiate", + simple $ \x -> + G.integrate (G.differentiate x) == (zero, x)) : + ("fourier differentiate", + simple $ \x -> + G.fourier (G.differentiate x) == + let y = G.fourier x + in y{G.polynomial = + Poly.fromCoeffs [0, 0 Complex.+: 2] * G.polynomial y}) : + ("approximate by bells, translate", + simple $ \x unit d -> unit/=0 ==> + G.approximateByBells unit (G.translateComplex d x) == + map (mapSnd (B.translateComplex d)) (G.approximateByBells unit x)) : + ("approximate by bells, dilate", + simple $ \x unit d -> unit/=0 && d/=0 ==> + G.approximateByBells unit (G.dilate d x) == + map (mapSnd (B.dilate d)) (G.approximateByBells (unit/d) x)) : + ("approximate by bells, shrink", + simple $ \x unit d -> unit/=0 && d/=0 ==> + G.approximateByBells unit (G.shrink d x) == + map (mapSnd (B.shrink d)) (G.approximateByBells (unit*d) x)) : + ("approximate by bells, different implementations", + quickCheck $ \unit d s p -> unit/=0 ==> + G.approximateByBellsAtOnce unit d s (p::Poly.T (Complex.T Rational)) == + G.approximateByBellsByTranslation unit d s p) : + [] diff --git a/test/Test/MathObj/Matrix.hs b/test/Test/MathObj/Matrix.hs new file mode 100644 --- /dev/null +++ b/test/Test/MathObj/Matrix.hs @@ -0,0 +1,96 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} +module Test.MathObj.Matrix where + +import qualified MathObj.Matrix as Matrix + +import qualified Algebra.Ring as Ring + +import qualified Algebra.Laws as Laws + +import qualified Number.NonNegative as NonNeg + +import qualified System.Random as Random + +import Test.NumericPrelude.Utility (testUnit, ) +import Test.QuickCheck (quickCheck, ) +import qualified Test.HUnit as HUnit + + +import PreludeBase as P +import NumericPrelude as NP + + +type Seed = Int +type Dimension = NonNeg.Int + +random :: Dimension -> Dimension -> Seed -> Matrix.T Integer +random mn nn seed = + fst $ + Matrix.random (NonNeg.toNumber mn) (NonNeg.toNumber nn) $ + Random.mkStdGen seed + + +tests :: HUnit.Test +tests = + HUnit.TestLabel "matrix" $ + HUnit.TestList $ + map testUnit $ + ("dimension", + quickCheck (\m n a -> + (NonNeg.toNumber m, NonNeg.toNumber n) == Matrix.dimension (random m n a))) : + ("to and from rows", + quickCheck (\m n a' -> + let a = random m n a' + in a == Matrix.fromRows (NonNeg.toNumber m) (NonNeg.toNumber n) (Matrix.rows a))) : + ("to and from columns", + quickCheck (\m n a' -> + let a = random m n a' + in a == Matrix.fromColumns (NonNeg.toNumber m) (NonNeg.toNumber n) (Matrix.columns a))) : + ("transpose, rows, columns", + quickCheck (\m n a' -> + let a = random m n a' + in Matrix.rows a == Matrix.columns (Matrix.transpose a))) : + ("transpose, columns, rows", + quickCheck (\m n a' -> + let a = random m n a' + in Matrix.columns a == Matrix.rows (Matrix.transpose a))) : + ("addition, zero", + quickCheck (\m n a -> + Laws.identity (+) (Matrix.zero (NonNeg.toNumber m) (NonNeg.toNumber n)) (random m n a))) : + ("addition, commutative", + quickCheck (\m n a b -> + Laws.commutative (+) (random m n a) (random m n b))) : + ("addition, associative", + quickCheck (\m n a b c -> + Laws.associative (+) (random m n a) (random m n b) (random m n c))) : + ("addition, transpose", + quickCheck (\m n a b -> + Laws.homomorphism Matrix.transpose (+) (+) (random m n a) (random m n b))) : + ("one, diagonal", + quickCheck (\n' -> + let n = NonNeg.toNumber n' + in Matrix.one n == (Matrix.diagonal $ replicate n Ring.one :: Matrix.T Integer))) : + ("multiplication, one left", + quickCheck (\m n a -> + Laws.leftIdentity (*) (Matrix.one (NonNeg.toNumber m)) (random m n a))) : + ("multiplication, one right", + quickCheck (\m n a -> + Laws.rightIdentity (*) (Matrix.one (NonNeg.toNumber n)) (random m n a))) : + ("multiplication, associative", + quickCheck (\k l m n a b c -> + Laws.associative (*) (random k l a) (random l m b) (random m n c))) : + ("multiplication and addition, distributive left", + quickCheck (\l m n a b c -> + Laws.leftDistributive (*) (+) (random n l a) (random m n b) (random m n c))) : + ("multiplication and addition, distributive right", + quickCheck (\l m n a b c -> + Laws.rightDistributive (*) (+) (random l m a) (random m n b) (random m n c))) : + ("multiplication, transpose", + quickCheck (\l m n a b -> + Laws.homomorphism Matrix.transpose (*) (flip (*)) (random l m a) (random m n b))) : +{- + ("division", quickCheck (\x -> Integral.propInverse (x :: Poly.T Rational))) : +-} + [] diff --git a/test/Test/MathObj/Polynomial.hs b/test/Test/MathObj/Polynomial.hs --- a/test/Test/MathObj/Polynomial.hs +++ b/test/Test/MathObj/Polynomial.hs @@ -14,7 +14,7 @@ import qualified Data.List as List import Test.NumericPrelude.Utility (testUnit) -import Test.QuickCheck (Property, quickCheck, (==>)) +import Test.QuickCheck (Property, quickCheck, (==>), Testable, ) import qualified Test.HUnit as HUnit @@ -32,6 +32,13 @@ mul xs ys = Poly.equal (Poly.mul xs ys) (Poly.mulShear xs ys) +test :: Testable a => (Poly.T Integer -> a) -> IO () +test = quickCheck + +testRat :: Testable a => (Poly.T Rational -> a) -> IO () +testRat = quickCheck + + tests :: HUnit.Test tests = HUnit.TestLabel "polynomial" $ @@ -39,12 +46,12 @@ map testUnit $ ("tensor product", quickCheck (tensorProductTranspose :: [Integer] -> [Integer] -> Property)) : ("mul speed", quickCheck (mul :: [Integer] -> [Integer] -> Bool)) : - ("addition, zero", quickCheck (\x -> Laws.identity (+) zero (x :: Poly.T Integer))) : - ("addition, commutative", quickCheck (\x -> Laws.commutative (+) (x :: Poly.T Integer))) : - ("addition, associative", quickCheck (\x -> Laws.associative (+) (x :: Poly.T Integer))) : - ("multiplication, one", quickCheck (\x -> Laws.identity (*) one (x :: Poly.T Integer))) : - ("multiplication, commutative", quickCheck (\x -> Laws.commutative (*) (x :: Poly.T Integer))) : - ("multiplication, associative", quickCheck (\x -> Laws.associative (*) (x :: Poly.T Integer))) : - ("multiplication and addition, distributive", quickCheck (\x -> Laws.leftDistributive (*) (+) (x :: Poly.T Integer))) : - ("division", quickCheck (\x -> Integral.propInverse (x :: Poly.T Rational))) : + ("addition, zero", test (Laws.identity (+) zero)) : + ("addition, commutative", test (Laws.commutative (+))) : + ("addition, associative", test (Laws.associative (+))) : + ("multiplication, one", test (Laws.identity (*) one)) : + ("multiplication, commutative", test (Laws.commutative (*))) : + ("multiplication, associative", test (Laws.associative (*))) : + ("multiplication and addition, distributive", test (Laws.leftDistributive (*) (+))) : + ("division", testRat (Integral.propInverse)) : [] diff --git a/test/Test/Number/GaloisField2p32m5.hs b/test/Test/Number/GaloisField2p32m5.hs new file mode 100644 --- /dev/null +++ b/test/Test/Number/GaloisField2p32m5.hs @@ -0,0 +1,37 @@ +{-# LANGUAGE NoImplicitPrelude #-} +module Test.Number.GaloisField2p32m5 where + +import qualified Number.GaloisField2p32m5 as GF + +import qualified Algebra.Laws as Laws + +import Test.NumericPrelude.Utility (testUnit) +import Test.QuickCheck (Testable, quickCheck, (==>)) +import qualified Test.HUnit as HUnit + + +import PreludeBase as P +import NumericPrelude as NP + + +test :: Testable a => (GF.T -> a) -> IO () +test = quickCheck + + +tests :: HUnit.Test +tests = + HUnit.TestLabel "galois field 2^32-5" $ + HUnit.TestList $ + map testUnit $ + ("addition, zero", test (Laws.identity (+) zero)) : + ("addition, commutative", test (Laws.commutative (+))) : + ("addition, associative", test (Laws.associative (+))) : + ("addition, negate", test (Laws.inverse (+) negate zero)) : + ("addition, subtract", test (\x -> Laws.inverse (+) (x-) x)) : + ("multiplication, one", test (Laws.identity (*) one)) : + ("multiplication, commutative", test (Laws.commutative (*))) : + ("multiplication, associative", test (Laws.associative (*))) : + ("multiplication, recip", test (\y -> y /= 0 ==> Laws.inverse (*) recip one y)) : + ("multiplication, division", test (\y x -> y /= 0 ==> Laws.inverse (*) (x/) x y)) : + ("multiplication and addition, distributive", test (Laws.leftDistributive (*) (+))) : + [] diff --git a/test/Test/Run.hs b/test/Test/Run.hs --- a/test/Test/Run.hs +++ b/test/Test/Run.hs @@ -4,8 +4,10 @@ import qualified Test.MathObj.Gaussian.Variance as GaussVariance import qualified Test.MathObj.Gaussian.Bell as GaussBell import qualified Test.MathObj.PartialFraction as PartialFraction +import qualified Test.MathObj.Matrix as Matrix import qualified Test.MathObj.Polynomial as Polynomial import qualified Test.MathObj.PowerSeries as PowerSeries +import qualified Test.Number.GaloisField2p32m5 as GF import qualified Test.HUnit.Text as HUnitText import qualified Test.HUnit as HUnit @@ -16,7 +18,9 @@ GaussBell.tests : GaussPoly.tests : PartialFraction.tests : + Matrix.tests : Polynomial.tests : PowerSeries.tests : + GF.tests : []) return ()