added experimental module for modular operation

5 files changed, 245 insertions, 7 deletions

diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal
index 69ec8fa..6614633 100644
--- a/packages/base/hmatrix.cabal
+++ b/packages/base/hmatrix.cabal
@@ -75,6 +75,7 @@ library
                         Numeric.Sparse
                         Numeric.LinearAlgebra.Static.Internal
                         Numeric.LinearAlgebra.Util
+                        Numeric.LinearAlgebra.Util.Modular
 
     C-sources:          src/C/lapack-aux.c
                         src/C/vector-aux.c
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
index 44d18b3..a03159d 100644
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ b/packages/base/src/Data/Packed/Internal/Numeric.hs
@@ -24,7 +24,7 @@ module Data.Packed.Internal.Numeric (
     scalar, conj, scale, arctan2, cmap, cmod,
     atIndex, minIndex, maxIndex, minElement, maxElement,
     sumElements, prodElements,
-    step, cond, find, assoc, accum,
+    step, cond, find, assoc, accum, findV, assocV, accumV,
     Transposable(..), Linear(..), Testable(..),
     -- * Matrix product and related functions
     Product(..), udot,
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 4ba0c98..c7eb491 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -141,6 +141,7 @@ module Numeric.LinearAlgebra (
     Complexable, RealElement,
     RealOf, ComplexOf, SingleOf, DoubleOf,
     IndexOf,
+    Mod,
     Field,
 --    Normed,
     Transposable,
@@ -156,6 +157,7 @@ import Data.Packed.Numeric hiding ((<>), mul)
 import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth)
 import qualified Numeric.LinearAlgebra.Algorithms as A
 import Numeric.LinearAlgebra.Util
+import Numeric.LinearAlgebra.Util.Modular
 import Numeric.LinearAlgebra.Random
 import Numeric.Sparse((!#>))
 import Numeric.LinearAlgebra.Util.CG
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs
index 8bacb09..9935c15 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Data.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs
@@ -14,9 +14,10 @@ Basic data processing.
 module Numeric.LinearAlgebra.Data(
 
     -- * Vector
-    -- | 1D arrays are storable vectors from the vector package.
+    -- | 1D arrays are storable vectors from the vector package. There is no distinction
+    --   between row and column vectors.
 
-    vector, (|>),
+    fromList, toList, vector, (|>),
 
     -- * Matrix
 
@@ -26,12 +27,12 @@ module Numeric.LinearAlgebra.Data(
 
     size, rows, cols,
 
-    -- * conversion from\/to lists
-    fromList, toList,
+    -- * Conversion from\/to lists
+
     fromLists, toLists,
     row, col,
 
-    -- * conversions vector\/matrix
+    -- * Conversions vector\/matrix
 
     flatten, reshape, asRow, asColumn,
     fromRows, toRows, fromColumns, toColumns,
@@ -86,7 +87,7 @@ module Numeric.LinearAlgebra.Data(
     separable,
     fromArray2D,
     module Data.Complex,
-    I,
+    I,F,
     Vector, Matrix, GMatrix, nRows, nCols
 
 ) where
@@ -97,5 +98,6 @@ import Data.Packed.Numeric
 import Numeric.LinearAlgebra.Util hiding ((&),(#))
 import Data.Complex
 import Numeric.Sparse
+import Numeric.LinearAlgebra.Util.Modular
 
 
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs
new file mode 100644
index 0000000..18b635d
--- /dev/null
+++ b/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs
@@ -0,0 +1,233 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE TypeFamilies #-}
+
+
+{- |
+Module      :  Numeric.LinearAlgebra.Util.Modular
+Copyright   :  (c) Alberto Ruiz 2015
+License     :  BSD3
+Stability   :  experimental
+
+Proof of concept of statically checked modular arithmetic.
+
+-}
+
+module Numeric.LinearAlgebra.Util.Modular(
+    Mod, F
+) where
+
+import Data.Packed.Numeric
+import Numeric.LinearAlgebra.Util(Indexable(..),size)
+import GHC.TypeLits
+import Data.Proxy(Proxy)
+import Foreign.ForeignPtr(castForeignPtr)
+import Data.Vector.Storable(unsafeToForeignPtr, unsafeFromForeignPtr)
+import Foreign.Storable
+import Data.Ratio
+import Data.Packed.Internal.Matrix hiding (mat,size)
+import Data.Packed.Internal.Numeric
+
+
+-- | Wrapper with a phantom integer for statically checked modular arithmetic.
+newtype Mod (n :: Nat) t = Mod {unMod:: t}
+  deriving (Storable)
+
+instance KnownNat m => Enum (F m)
+  where
+    toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m))
+    fromEnum = fromIntegral . unMod
+
+instance KnownNat m => Eq (F m)
+  where
+    a == b = (unMod a) == (unMod b)
+
+instance KnownNat m => Ord (F m)
+  where
+    compare a b = compare (unMod a) (unMod b)
+
+instance KnownNat m => Real (F m)
+  where
+    toRational x = toInteger x % 1
+
+instance KnownNat m => Integral (F m)
+  where
+    toInteger = toInteger . unMod
+    quotRem a b = (Mod q, Mod r)
+      where
+         (q,r) = quotRem (unMod a) (unMod b)
+
+-- | this instance is only valid for prime m (not checked)
+instance KnownNat m => Fractional (F m)
+  where
+    recip x = x^(m'-2)
+      where
+        m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+    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 f (Mod u) (Mod v) = Mod (f m' u v)
+  where
+    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+
+l1 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> Mod m a -> Mod m b
+l1 f (Mod u) = Mod (f m' u)
+  where
+    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+
+l0 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> a -> Mod m b
+l0 f u = Mod (f m' u)
+  where
+    m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+
+
+instance Show (F n)
+  where
+    show = show . unMod
+
+instance forall n . KnownNat n => Num (F n)
+  where
+    (+) = l2 (\m a b -> (a + b) `mod` (fromIntegral m))
+    (*) = l2 (\m a b -> (a * b) `mod` (fromIntegral m))
+    (-) = l2 (\m a b -> (a - b) `mod` (fromIntegral m))
+    abs = l1 (const abs)
+    signum = l1 (const signum)
+    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)
+  where
+    transdata n v m = i2f (transdata n (f2i v) m)
+    constantD x n = i2f (constantD (unMod x) n)
+    extractR m mi is mj js = i2fM (extractR (f2iM m) mi is mj js)
+    sortI = sortI . f2i
+    sortV = i2f . sortV . f2i
+    compareV u v = compareV (f2i u) (f2i v)
+    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)
+  where
+    conj' = id
+    size' = size
+    scale' s x = fromInt (scale (unMod s) (toInt x))
+    addConstant c x = fromInt (addConstant (unMod c) (toInt x))
+    add a b = fromInt (add (toInt a) (toInt b))
+    sub a b = fromInt (sub (toInt a) (toInt b))
+    mul a b = fromInt (mul (toInt a) (toInt b))
+    equal u v = equal (toInt u) (toInt 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 (toInt x) k)
+    minIndex'     = minIndex . toInt
+    maxIndex'     = maxIndex . toInt
+    minElement'   = Mod . minElement . toInt
+    maxElement'   = Mod . maxElement . toInt
+    sumElements'  = fromIntegral . sumElements . toInt
+    prodElements' = fromIntegral . sumElements . toInt
+    step'         = i2f . step . toInt
+    find' = findV
+    assoc' = assocV
+    accum' = accumV
+    cond' x y l e g = cselect (ccompare x y) l e g
+    ccompare' a b = ccompare (toInt a) (toInt b)
+    cselect' c l e g = i2f $ cselect c (toInt l) (toInt e) (toInt g)
+    scaleRecip s x = scale' s (cmap recip x)
+    divide x y = mul x (cmap recip y)
+    arctan2' = undefined
+    cmod' m = fromInt' . cmod' (unMod m) . toInt'
+    fromInt' v = i2f $ cmod' (fromIntegral m') (fromInt' v)
+      where
+        m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
+    toInt'   = f2i
+
+
+instance Indexable (Vector (F m)) (F m)
+  where
+    (!) = (@>)
+
+
+type instance RealOf (F n) = I
+
+
+instance KnownNat m => Product (F m) where
+    norm2      = undefined
+    absSum     = undefined
+    norm1      = undefined
+    normInf    = undefined
+    multiply   = lift2 multiply
+
+
+instance KnownNat m => Numeric (F m)
+
+i2f :: Vector I -> Vector (F n)
+i2f v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n)
+    where (fp,i,n) = unsafeToForeignPtr v
+
+f2i :: Vector (F n) -> Vector I
+f2i v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n)
+    where (fp,i,n) = unsafeToForeignPtr v
+
+f2iM :: Matrix (F n) -> Matrix I
+f2iM = liftMatrix f2i
+
+i2fM :: Matrix I -> Matrix (F n)
+i2fM = liftMatrix i2f
+
+
+lift1 f a   = fromInt (f (toInt a))
+lift2 f a b = fromInt (f (toInt a) (toInt b))
+
+instance forall m . KnownNat m => Num (V m)
+  where
+    (+) = lift2 (+)
+    (*) = lift2 (*)
+    (-) = lift2 (-)
+    abs = lift1 abs
+    signum = lift1 signum
+    negate = lift1 negate
+    fromInteger x = fromInt (fromInteger x)
+
+
+--------------------------------------------------------------------------------
+
+instance (KnownNat m) => Testable (M m)
+  where
+    checkT _ = test
+
+test = (ok, info)
+  where
+    v = fromList [3,-5,75] :: V 11
+    m = (3><3) [1..]   :: M 11
+    info = do
+        print v
+        print m
+        print (tr m)
+        print $ v+v
+        print $ m+m
+        print $ m <> m
+        print $ m #> v
+
+
+    ok = and
+      [ toInt (m #> v) == cmod 11 (toInt m #> toInt v )
+      ]
+
+