From 7c446e77a52dc0a979c0fe570f1e7b6127d9ff42 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 8 Jun 2015 12:07:15 +0200 Subject: instances for Mod m Z and Mod m I --- packages/base/src/Internal/C/lapack-aux.c | 2 +- packages/base/src/Internal/C/vector-aux.c | 8 +- packages/base/src/Internal/Modular.hs | 136 ++++++++++++++++------ packages/base/src/Internal/Numeric.hs | 4 +- packages/base/src/Numeric/LinearAlgebra.hs | 2 - packages/base/src/Numeric/LinearAlgebra/Data.hs | 2 +- packages/tests/src/Numeric/LinearAlgebra/Tests.hs | 2 +- 7 files changed, 108 insertions(+), 48 deletions(-) (limited to 'packages') diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c index 1601bef..dcce1c5 100644 --- a/packages/base/src/Internal/C/lapack-aux.c +++ b/packages/base/src/Internal/C/lapack-aux.c @@ -1308,7 +1308,7 @@ int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) { } OK } int multiplyI(int m, KOIMAT(a), KOIMAT(b), OIMAT(r)) MULT_IMP -int multiplyL(int32_t m, KOLMAT(a), KOLMAT(b), OLMAT(r)) MULT_IMP +int multiplyL(int64_t m, KOLMAT(a), KOLMAT(b), OLMAT(r)) MULT_IMP ////////////////// sparse matrix-product /////////////////////////////////////// diff --git a/packages/base/src/Internal/C/vector-aux.c b/packages/base/src/Internal/C/vector-aux.c index 580aa1c..c161556 100644 --- a/packages/base/src/Internal/C/vector-aux.c +++ b/packages/base/src/Internal/C/vector-aux.c @@ -71,7 +71,7 @@ int sumI(int m, KIVEC(x),IVEC(r)) { OK } -int sumL(int32_t m, KLVEC(x),LVEC(r)) { +int sumL(int64_t m, KLVEC(x),LVEC(r)) { REQUIRES(rn==1,BAD_SIZE); int i; int res = 0; @@ -148,7 +148,7 @@ int prodI(int m, KIVEC(x),IVEC(r)) { OK } -int prodL(int32_t m, KLVEC(x),LVEC(r)) { +int prodL(int64_t m, KLVEC(x),LVEC(r)) { REQUIRES(rn==1,BAD_SIZE); int i; int res = 1; @@ -761,8 +761,8 @@ int mapValL(int code, int64_t* pval, KLVEC(x), LVEC(r)) { OPV(1,val/xp[k]) OPV(2,val+xp[k]) OPV(3,val-xp[k]) - OPV(6,mod(val,xp[k])) - OPV(7,mod(xp[k],val)) + OPV(6,mod_l(val,xp[k])) + OPV(7,mod_l(xp[k],val)) default: ERROR(BAD_CODE); } } diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs index 0274607..6b34010 100644 --- a/packages/base/src/Internal/Modular.hs +++ b/packages/base/src/Internal/Modular.hs @@ -7,6 +7,7 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE TypeFamilies #-} @@ -22,7 +23,7 @@ Proof of concept of statically checked modular arithmetic. -} module Internal.Modular( - Mod, F + Mod ) where import Internal.Vector @@ -30,8 +31,8 @@ import Internal.Matrix hiding (mat,size) import Internal.Numeric import Internal.Element import Internal.Container -import Internal.Vectorized (prodI,sumI) -import Internal.LAPACK (multiplyI) +import Internal.Vectorized (prodI,sumI,prodL,sumL) +import Internal.LAPACK (multiplyI, multiplyL) import Internal.Util(Indexable(..),gaussElim) import GHC.TypeLits import Data.Proxy(Proxy) @@ -45,24 +46,24 @@ import Data.Ratio newtype Mod (n :: Nat) t = Mod {unMod:: t} deriving (Storable) -instance KnownNat m => Enum (F m) +instance (Integral t, Enum t, KnownNat m) => Enum (Mod m t) where toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m)) fromEnum = fromIntegral . unMod -instance KnownNat m => Eq (F m) +instance (Eq t, KnownNat m) => Eq (Mod m t) where a == b = (unMod a) == (unMod b) -instance KnownNat m => Ord (F m) +instance (Ord t, KnownNat m) => Ord (Mod m t) where compare a b = compare (unMod a) (unMod b) -instance KnownNat m => Real (F m) +instance (Real t, KnownNat m, Integral (Mod m t)) => Real (Mod m t) where toRational x = toInteger x % 1 -instance KnownNat m => Integral (F m) +instance (Integral t, KnownNat m, Num (Mod m t)) => Integral (Mod m t) where toInteger = toInteger . unMod quotRem a b = (Mod q, Mod r) @@ -70,7 +71,7 @@ instance KnownNat m => Integral (F m) (q,r) = quotRem (unMod a) (unMod b) -- | this instance is only valid for prime m -instance KnownNat m => Fractional (F m) +instance (Show (Mod m t), Num (Mod m t), Eq t, KnownNat m) => Fractional (Mod m t) where recip x | x*r == 1 = r @@ -80,27 +81,27 @@ instance KnownNat m => Fractional (F m) m' = fromIntegral . natVal $ (undefined :: Proxy m) fromRational x = fromInteger (numerator x) / fromInteger (denominator x) -l2 :: forall m a b c. (KnownNat m) => (Int -> a -> b -> c) -> Mod m a -> Mod m b -> Mod m c +l2 :: forall m a b c. (Num c, KnownNat m) => (c -> a -> b -> c) -> Mod m a -> Mod m b -> Mod m c l2 f (Mod u) (Mod v) = Mod (f m' u v) where m' = fromIntegral . natVal $ (undefined :: Proxy m) -l1 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> Mod m a -> Mod m b +l1 :: forall m a b . (Num b, KnownNat m) => (b -> a -> b) -> Mod m a -> Mod m b l1 f (Mod u) = Mod (f m' u) where m' = fromIntegral . natVal $ (undefined :: Proxy m) -l0 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> a -> Mod m b +l0 :: forall m a b . (Num b, KnownNat m) => (b -> a -> b) -> a -> Mod m b l0 f u = Mod (f m' u) where m' = fromIntegral . natVal $ (undefined :: Proxy m) -instance Show (F n) +instance Show t => Show (Mod n t) where show = show . unMod -instance forall n . KnownNat n => Num (F n) +instance forall n t . (Integral t, KnownNat n) => Num (Mod n t) where (+) = l2 (\m a b -> (a + b) `mod` (fromIntegral m)) (*) = l2 (\m a b -> (a * b) `mod` (fromIntegral m)) @@ -110,14 +111,8 @@ instance forall n . KnownNat n => Num (F n) fromInteger = l0 (\m x -> fromInteger x `mod` (fromIntegral m)) --- | Integer modulo n -type F n = Mod n I -type V n = Vector (F n) -type M n = Matrix (F n) - - -instance Element (F n) +instance (Ord t, Element t) => Element (Mod n t) where transdata n v m = i2f (transdata n (f2i v) m) constantD x n = i2f (constantD (unMod x) n) @@ -128,7 +123,8 @@ instance Element (F n) selectV c l e g = i2f (selectV c (f2i l) (f2i e) (f2i g)) remapM i j m = i2fM (remap i j (f2iM m)) -instance forall m . KnownNat m => Container Vector (F m) + +instance forall m . KnownNat m => Container Vector (Mod m I) where conj' = id size' = dim @@ -168,24 +164,77 @@ instance forall m . KnownNat m => Container Vector (F m) fromZ' = vmod . fromZ' toZ' = toZ' . f2i +instance forall m . KnownNat m => Container Vector (Mod m Z) + where + conj' = id + size' = dim + scale' s x = vmod (scale (unMod s) (f2i x)) + addConstant c x = vmod (addConstant (unMod c) (f2i x)) + add a b = vmod (add (f2i a) (f2i b)) + sub a b = vmod (sub (f2i a) (f2i b)) + mul a b = vmod (mul (f2i a) (f2i b)) + equal u v = equal (f2i u) (f2i v) + scalar' x = fromList [x] + konst' x = i2f . konst (unMod x) + build' n f = build n (fromIntegral . f) + cmap' = cmap + atIndex' x k = fromIntegral (atIndex (f2i x) k) + minIndex' = minIndex . f2i + maxIndex' = maxIndex . f2i + minElement' = Mod . minElement . f2i + maxElement' = Mod . maxElement . f2i + sumElements' = fromIntegral . sumL m' . f2i + where + m' = fromIntegral . natVal $ (undefined :: Proxy m) + prodElements' = fromIntegral . prodL m' . f2i + where + m' = fromIntegral . natVal $ (undefined :: Proxy m) + step' = i2f . step . f2i + find' = findV + assoc' = assocV + accum' = accumV + ccompare' a b = ccompare (f2i a) (f2i b) + cselect' c l e g = i2f $ cselect c (f2i l) (f2i e) (f2i g) + scaleRecip s x = scale' s (cmap recip x) + divide x y = mul x (cmap recip y) + arctan2' = undefined + cmod' m = vmod . cmod' (unMod m) . f2i + fromInt' = vmod . fromInt' + toInt' = toInt . f2i + fromZ' = vmod + toZ' = f2i + + -instance Indexable (Vector (F m)) (F m) +instance (Storable t, Indexable (Vector t) t) => Indexable (Vector (Mod m t)) (Mod m t) where (!) = (@>) -type instance RealOf (F n) = I +type instance RealOf (Mod n I) = I +type instance RealOf (Mod n Z) = Z -instance KnownNat m => Product (F m) where +instance KnownNat m => Product (Mod m I) where norm2 = undefined absSum = undefined norm1 = undefined normInf = undefined - multiply = lift2 (multiplyI m') + multiply = lift2m (multiplyI m') + where + m' = fromIntegral . natVal $ (undefined :: Proxy m) + +instance KnownNat m => Product (Mod m Z) where + norm2 = undefined + absSum = undefined + norm1 = undefined + normInf = undefined + multiply = lift2m (multiplyL m') where m' = fromIntegral . natVal $ (undefined :: Proxy m) -instance KnownNat m => Numeric (F m) + +instance KnownNat m => Numeric (Mod m I) +instance KnownNat m => Numeric (Mod m Z) i2f :: Storable t => Vector t -> Vector (Mod n t) i2f v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n) @@ -206,10 +255,12 @@ vmod = i2f . cmod' m' where m' = fromIntegral . natVal $ (undefined :: Proxy m) -lift1 f a = fromInt (f (toInt a)) -lift2 f a b = fromInt (f (toInt a) (toInt b)) +lift1 f a = vmod (f (f2i a)) +lift2 f a b = vmod (f (f2i a) (f2i b)) + +lift2m f a b = liftMatrix vmod (f (f2iM a) (f2iM b)) -instance forall m . KnownNat m => Num (V m) +instance forall m . KnownNat m => Num (Vector (Mod m I)) where (+) = lift2 (+) (*) = lift2 (*) @@ -222,14 +273,14 @@ instance forall m . KnownNat m => Num (V m) -------------------------------------------------------------------------------- -instance (KnownNat m) => Testable (M m) +instance (KnownNat m) => Testable (Matrix (Mod m I)) where checkT _ = test test = (ok, info) where - v = fromList [3,-5,75] :: V 11 - m = (3><3) [1..] :: M 11 + v = fromList [3,-5,75] :: Vector (Mod 11 I) + m = (3><3) [1..] :: Matrix (Mod 11 I) a = (3><3) [1,2 , 3 ,4,5 , 6 @@ -237,13 +288,17 @@ test = (ok, info) b = (3><2) [0..] :: Matrix I - am = fromInt a :: Matrix (F 13) - bm = fromInt b :: Matrix (F 13) + am = fromInt a :: Matrix (Mod 13 I) + bm = fromInt b :: Matrix (Mod 13 I) ad = fromInt a :: Matrix Double bd = fromInt b :: Matrix Double g = (3><3) (repeat (40000)) :: Matrix I - gm = fromInt g :: Matrix (F 100000) + gm = fromInt g :: Matrix (Mod 100000 I) + + lg = (3><3) (repeat (3*10^(9::Int))) :: Matrix Z + lgm = fromZ lg :: Matrix (Mod 10000000000 Z) + info = do print v @@ -262,11 +317,18 @@ test = (ok, info) print gm print $ gm <> gm + print lg + print $ lg <> lg + print lgm + print $ lgm <> lgm + + ok = and [ toInt (m #> v) == cmod 11 (toInt m #> toInt v ) , am <> gaussElim am bm == bm - , prodElements (konst (9:: F 10) (12::Int)) == product (replicate 12 (9:: F 10)) + , prodElements (konst (9:: Mod 10 I) (12::Int)) == product (replicate 12 (9:: Mod 10 I)) , gm <> gm == konst 0 (3,3) + , lgm <> lgm == konst 0 (3,3) ] diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs index 2ef96bf..ca17c23 100644 --- a/packages/base/src/Internal/Numeric.hs +++ b/packages/base/src/Internal/Numeric.hs @@ -398,8 +398,8 @@ arctan2 = arctan2' -- -- >>> cmod 3 (range 5) -- fromList [0,1,2,0,1] -cmod :: (Integral e, Container c e) => Int -> c e -> c e -cmod m = cmod' (fromIntegral m) +cmod :: (Integral e, Container c e) => e -> c e -> c e +cmod = cmod' -- | -- >>>fromInt ((2><2) [0..3]) :: Matrix (Complex Double) diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index ace47e9..56e5053 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -141,7 +141,6 @@ module Numeric.LinearAlgebra ( Complexable, RealElement, RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf, - Mod, Field, -- Normed, Transposable, @@ -159,7 +158,6 @@ import Internal.Numeric hiding (mul) import Internal.Algorithms hiding (linearSolve,Normed,orth) import qualified Internal.Algorithms as A import Internal.Util -import Internal.Modular import Internal.Random import Internal.Sparse((!#>)) import Internal.CG diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index c3e0c1d..1c9bb68 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs @@ -88,7 +88,7 @@ module Numeric.LinearAlgebra.Data( separable, fromArray2D, module Data.Complex, - R,C,I,Z,F, + R,C,I,Z,Mod, Vector, Matrix, GMatrix, nRows, nCols ) where diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs index 2be75de..4cfa5cb 100644 --- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs +++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs @@ -383,7 +383,7 @@ intTest = utest "int ops" (fst $ checkT (undefined :: Matrix I)) -------------------------------------------------------------------------------- -modularTest = utest "modular ops" (fst $ checkT (undefined :: Matrix (F 13))) +modularTest = utest "modular ops" (fst $ checkT (undefined :: Matrix (Mod 13 I))) -------------------------------------------------------------------------------- -- cgit v1.2.3