summaryrefslogtreecommitdiff
path: root/packages/base/src/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Numeric')
-rw-r--r--packages/base/src/Numeric/LinearAlgebra.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Data.hs14
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs233
3 files changed, 243 insertions, 6 deletions
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 (
141 Complexable, RealElement, 141 Complexable, RealElement,
142 RealOf, ComplexOf, SingleOf, DoubleOf, 142 RealOf, ComplexOf, SingleOf, DoubleOf,
143 IndexOf, 143 IndexOf,
144 Mod,
144 Field, 145 Field,
145-- Normed, 146-- Normed,
146 Transposable, 147 Transposable,
@@ -156,6 +157,7 @@ import Data.Packed.Numeric hiding ((<>), mul)
156import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth) 157import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth)
157import qualified Numeric.LinearAlgebra.Algorithms as A 158import qualified Numeric.LinearAlgebra.Algorithms as A
158import Numeric.LinearAlgebra.Util 159import Numeric.LinearAlgebra.Util
160import Numeric.LinearAlgebra.Util.Modular
159import Numeric.LinearAlgebra.Random 161import Numeric.LinearAlgebra.Random
160import Numeric.Sparse((!#>)) 162import Numeric.Sparse((!#>))
161import Numeric.LinearAlgebra.Util.CG 163import 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.
14module Numeric.LinearAlgebra.Data( 14module Numeric.LinearAlgebra.Data(
15 15
16 -- * Vector 16 -- * Vector
17 -- | 1D arrays are storable vectors from the vector package. 17 -- | 1D arrays are storable vectors from the vector package. There is no distinction
18 -- between row and column vectors.
18 19
19 vector, (|>), 20 fromList, toList, vector, (|>),
20 21
21 -- * Matrix 22 -- * Matrix
22 23
@@ -26,12 +27,12 @@ module Numeric.LinearAlgebra.Data(
26 27
27 size, rows, cols, 28 size, rows, cols,
28 29
29 -- * conversion from\/to lists 30 -- * Conversion from\/to lists
30 fromList, toList, 31
31 fromLists, toLists, 32 fromLists, toLists,
32 row, col, 33 row, col,
33 34
34 -- * conversions vector\/matrix 35 -- * Conversions vector\/matrix
35 36
36 flatten, reshape, asRow, asColumn, 37 flatten, reshape, asRow, asColumn,
37 fromRows, toRows, fromColumns, toColumns, 38 fromRows, toRows, fromColumns, toColumns,
@@ -86,7 +87,7 @@ module Numeric.LinearAlgebra.Data(
86 separable, 87 separable,
87 fromArray2D, 88 fromArray2D,
88 module Data.Complex, 89 module Data.Complex,
89 I, 90 I,F,
90 Vector, Matrix, GMatrix, nRows, nCols 91 Vector, Matrix, GMatrix, nRows, nCols
91 92
92) where 93) where
@@ -97,5 +98,6 @@ import Data.Packed.Numeric
97import Numeric.LinearAlgebra.Util hiding ((&),(#)) 98import Numeric.LinearAlgebra.Util hiding ((&),(#))
98import Data.Complex 99import Data.Complex
99import Numeric.Sparse 100import Numeric.Sparse
101import Numeric.LinearAlgebra.Util.Modular
100 102
101 103
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 @@
1{-# LANGUAGE DataKinds #-}
2{-# LANGUAGE KindSignatures #-}
3{-# LANGUAGE GeneralizedNewtypeDeriving #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE FlexibleContexts #-}
7{-# LANGUAGE ScopedTypeVariables #-}
8{-# LANGUAGE Rank2Types #-}
9{-# LANGUAGE FlexibleInstances #-}
10{-# LANGUAGE GADTs #-}
11{-# LANGUAGE TypeFamilies #-}
12
13
14{- |
15Module : Numeric.LinearAlgebra.Util.Modular
16Copyright : (c) Alberto Ruiz 2015
17License : BSD3
18Stability : experimental
19
20Proof of concept of statically checked modular arithmetic.
21
22-}
23
24module Numeric.LinearAlgebra.Util.Modular(
25 Mod, F
26) where
27
28import Data.Packed.Numeric
29import Numeric.LinearAlgebra.Util(Indexable(..),size)
30import GHC.TypeLits
31import Data.Proxy(Proxy)
32import Foreign.ForeignPtr(castForeignPtr)
33import Data.Vector.Storable(unsafeToForeignPtr, unsafeFromForeignPtr)
34import Foreign.Storable
35import Data.Ratio
36import Data.Packed.Internal.Matrix hiding (mat,size)
37import Data.Packed.Internal.Numeric
38
39
40-- | Wrapper with a phantom integer for statically checked modular arithmetic.
41newtype Mod (n :: Nat) t = Mod {unMod:: t}
42 deriving (Storable)
43
44instance KnownNat m => Enum (F m)
45 where
46 toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m))
47 fromEnum = fromIntegral . unMod
48
49instance KnownNat m => Eq (F m)
50 where
51 a == b = (unMod a) == (unMod b)
52
53instance KnownNat m => Ord (F m)
54 where
55 compare a b = compare (unMod a) (unMod b)
56
57instance KnownNat m => Real (F m)
58 where
59 toRational x = toInteger x % 1
60
61instance KnownNat m => Integral (F m)
62 where
63 toInteger = toInteger . unMod
64 quotRem a b = (Mod q, Mod r)
65 where
66 (q,r) = quotRem (unMod a) (unMod b)
67
68-- | this instance is only valid for prime m (not checked)
69instance KnownNat m => Fractional (F m)
70 where
71 recip x = x^(m'-2)
72 where
73 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
74 fromRational x = fromInteger (numerator x) / fromInteger (denominator x)
75
76l2 :: forall m a b c. (KnownNat m) => (Int -> a -> b -> c) -> Mod m a -> Mod m b -> Mod m c
77l2 f (Mod u) (Mod v) = Mod (f m' u v)
78 where
79 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
80
81l1 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> Mod m a -> Mod m b
82l1 f (Mod u) = Mod (f m' u)
83 where
84 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
85
86l0 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> a -> Mod m b
87l0 f u = Mod (f m' u)
88 where
89 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
90
91
92instance Show (F n)
93 where
94 show = show . unMod
95
96instance forall n . KnownNat n => Num (F n)
97 where
98 (+) = l2 (\m a b -> (a + b) `mod` (fromIntegral m))
99 (*) = l2 (\m a b -> (a * b) `mod` (fromIntegral m))
100 (-) = l2 (\m a b -> (a - b) `mod` (fromIntegral m))
101 abs = l1 (const abs)
102 signum = l1 (const signum)
103 fromInteger = l0 (\m x -> fromInteger x `mod` (fromIntegral m))
104
105
106-- | Integer modulo n
107type F n = Mod n I
108
109type V n = Vector (F n)
110type M n = Matrix (F n)
111
112
113instance Element (F n)
114 where
115 transdata n v m = i2f (transdata n (f2i v) m)
116 constantD x n = i2f (constantD (unMod x) n)
117 extractR m mi is mj js = i2fM (extractR (f2iM m) mi is mj js)
118 sortI = sortI . f2i
119 sortV = i2f . sortV . f2i
120 compareV u v = compareV (f2i u) (f2i v)
121 selectV c l e g = i2f (selectV c (f2i l) (f2i e) (f2i g))
122 remapM i j m = i2fM (remap i j (f2iM m))
123
124instance forall m . KnownNat m => Container Vector (F m)
125 where
126 conj' = id
127 size' = size
128 scale' s x = fromInt (scale (unMod s) (toInt x))
129 addConstant c x = fromInt (addConstant (unMod c) (toInt x))
130 add a b = fromInt (add (toInt a) (toInt b))
131 sub a b = fromInt (sub (toInt a) (toInt b))
132 mul a b = fromInt (mul (toInt a) (toInt b))
133 equal u v = equal (toInt u) (toInt v)
134 scalar' x = fromList [x]
135 konst' x = i2f . konst (unMod x)
136 build' n f = build n (fromIntegral . f)
137 cmap' = cmap
138 atIndex' x k = fromIntegral (atIndex (toInt x) k)
139 minIndex' = minIndex . toInt
140 maxIndex' = maxIndex . toInt
141 minElement' = Mod . minElement . toInt
142 maxElement' = Mod . maxElement . toInt
143 sumElements' = fromIntegral . sumElements . toInt
144 prodElements' = fromIntegral . sumElements . toInt
145 step' = i2f . step . toInt
146 find' = findV
147 assoc' = assocV
148 accum' = accumV
149 cond' x y l e g = cselect (ccompare x y) l e g
150 ccompare' a b = ccompare (toInt a) (toInt b)
151 cselect' c l e g = i2f $ cselect c (toInt l) (toInt e) (toInt g)
152 scaleRecip s x = scale' s (cmap recip x)
153 divide x y = mul x (cmap recip y)
154 arctan2' = undefined
155 cmod' m = fromInt' . cmod' (unMod m) . toInt'
156 fromInt' v = i2f $ cmod' (fromIntegral m') (fromInt' v)
157 where
158 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
159 toInt' = f2i
160
161
162instance Indexable (Vector (F m)) (F m)
163 where
164 (!) = (@>)
165
166
167type instance RealOf (F n) = I
168
169
170instance KnownNat m => Product (F m) where
171 norm2 = undefined
172 absSum = undefined
173 norm1 = undefined
174 normInf = undefined
175 multiply = lift2 multiply
176
177
178instance KnownNat m => Numeric (F m)
179
180i2f :: Vector I -> Vector (F n)
181i2f v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n)
182 where (fp,i,n) = unsafeToForeignPtr v
183
184f2i :: Vector (F n) -> Vector I
185f2i v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n)
186 where (fp,i,n) = unsafeToForeignPtr v
187
188f2iM :: Matrix (F n) -> Matrix I
189f2iM = liftMatrix f2i
190
191i2fM :: Matrix I -> Matrix (F n)
192i2fM = liftMatrix i2f
193
194
195lift1 f a = fromInt (f (toInt a))
196lift2 f a b = fromInt (f (toInt a) (toInt b))
197
198instance forall m . KnownNat m => Num (V m)
199 where
200 (+) = lift2 (+)
201 (*) = lift2 (*)
202 (-) = lift2 (-)
203 abs = lift1 abs
204 signum = lift1 signum
205 negate = lift1 negate
206 fromInteger x = fromInt (fromInteger x)
207
208
209--------------------------------------------------------------------------------
210
211instance (KnownNat m) => Testable (M m)
212 where
213 checkT _ = test
214
215test = (ok, info)
216 where
217 v = fromList [3,-5,75] :: V 11
218 m = (3><3) [1..] :: M 11
219 info = do
220 print v
221 print m
222 print (tr m)
223 print $ v+v
224 print $ m+m
225 print $ m <> m
226 print $ m #> v
227
228
229 ok = and
230 [ toInt (m #> v) == cmod 11 (toInt m #> toInt v )
231 ]
232
233