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/Util/Modular.hs251
1 files changed, 0 insertions, 251 deletions
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs
deleted file mode 100644
index ea4a668..0000000
--- a/packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs
+++ /dev/null
@@ -1,251 +0,0 @@
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(..),gaussElim)
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
69instance KnownNat m => Fractional (F m)
70 where
71 recip x
72 | x*r == 1 = r
73 | otherwise = error $ show x ++" does not have a multiplicative inverse mod "++show m'
74 where
75 r = x^(m'-2)
76 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
77 fromRational x = fromInteger (numerator x) / fromInteger (denominator x)
78
79l2 :: forall m a b c. (KnownNat m) => (Int -> a -> b -> c) -> Mod m a -> Mod m b -> Mod m c
80l2 f (Mod u) (Mod v) = Mod (f m' u v)
81 where
82 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
83
84l1 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> Mod m a -> Mod m b
85l1 f (Mod u) = Mod (f m' u)
86 where
87 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
88
89l0 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> a -> Mod m b
90l0 f u = Mod (f m' u)
91 where
92 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
93
94
95instance Show (F n)
96 where
97 show = show . unMod
98
99instance forall n . KnownNat n => Num (F n)
100 where
101 (+) = l2 (\m a b -> (a + b) `mod` (fromIntegral m))
102 (*) = l2 (\m a b -> (a * b) `mod` (fromIntegral m))
103 (-) = l2 (\m a b -> (a - b) `mod` (fromIntegral m))
104 abs = l1 (const abs)
105 signum = l1 (const signum)
106 fromInteger = l0 (\m x -> fromInteger x `mod` (fromIntegral m))
107
108
109-- | Integer modulo n
110type F n = Mod n I
111
112type V n = Vector (F n)
113type M n = Matrix (F n)
114
115
116instance Element (F n)
117 where
118 transdata n v m = i2f (transdata n (f2i v) m)
119 constantD x n = i2f (constantD (unMod x) n)
120 extractR m mi is mj js = i2fM (extractR (f2iM m) mi is mj js)
121 sortI = sortI . f2i
122 sortV = i2f . sortV . f2i
123 compareV u v = compareV (f2i u) (f2i v)
124 selectV c l e g = i2f (selectV c (f2i l) (f2i e) (f2i g))
125 remapM i j m = i2fM (remap i j (f2iM m))
126
127instance forall m . KnownNat m => Container Vector (F m)
128 where
129 conj' = id
130 size' = dim
131 scale' s x = fromInt (scale (unMod s) (toInt x))
132 addConstant c x = fromInt (addConstant (unMod c) (toInt x))
133 add a b = fromInt (add (toInt a) (toInt b))
134 sub a b = fromInt (sub (toInt a) (toInt b))
135 mul a b = fromInt (mul (toInt a) (toInt b))
136 equal u v = equal (toInt u) (toInt v)
137 scalar' x = fromList [x]
138 konst' x = i2f . konst (unMod x)
139 build' n f = build n (fromIntegral . f)
140 cmap' = cmap
141 atIndex' x k = fromIntegral (atIndex (toInt x) k)
142 minIndex' = minIndex . toInt
143 maxIndex' = maxIndex . toInt
144 minElement' = Mod . minElement . toInt
145 maxElement' = Mod . maxElement . toInt
146 sumElements' = fromIntegral . sumElements . toInt
147 prodElements' = fromIntegral . sumElements . toInt
148 step' = i2f . step . toInt
149 find' = findV
150 assoc' = assocV
151 accum' = accumV
152 cond' x y l e g = cselect (ccompare x y) l e g
153 ccompare' a b = ccompare (toInt a) (toInt b)
154 cselect' c l e g = i2f $ cselect c (toInt l) (toInt e) (toInt g)
155 scaleRecip s x = scale' s (cmap recip x)
156 divide x y = mul x (cmap recip y)
157 arctan2' = undefined
158 cmod' m = fromInt' . cmod' (unMod m) . toInt'
159 fromInt' v = i2f $ cmod' (fromIntegral m') (fromInt' v)
160 where
161 m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int
162 toInt' = f2i
163
164
165instance Indexable (Vector (F m)) (F m)
166 where
167 (!) = (@>)
168
169
170type instance RealOf (F n) = I
171
172
173instance KnownNat m => Product (F m) where
174 norm2 = undefined
175 absSum = undefined
176 norm1 = undefined
177 normInf = undefined
178 multiply = lift2 multiply
179
180
181instance KnownNat m => Numeric (F m)
182
183i2f :: Vector I -> Vector (F n)
184i2f v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n)
185 where (fp,i,n) = unsafeToForeignPtr v
186
187f2i :: Vector (F n) -> Vector I
188f2i v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n)
189 where (fp,i,n) = unsafeToForeignPtr v
190
191f2iM :: Matrix (F n) -> Matrix I
192f2iM = liftMatrix f2i
193
194i2fM :: Matrix I -> Matrix (F n)
195i2fM = liftMatrix i2f
196
197
198lift1 f a = fromInt (f (toInt a))
199lift2 f a b = fromInt (f (toInt a) (toInt b))
200
201instance forall m . KnownNat m => Num (V m)
202 where
203 (+) = lift2 (+)
204 (*) = lift2 (*)
205 (-) = lift2 (-)
206 abs = lift1 abs
207 signum = lift1 signum
208 negate = lift1 negate
209 fromInteger x = fromInt (fromInteger x)
210
211
212--------------------------------------------------------------------------------
213
214instance (KnownNat m) => Testable (M m)
215 where
216 checkT _ = test
217
218test = (ok, info)
219 where
220 v = fromList [3,-5,75] :: V 11
221 m = (3><3) [1..] :: M 11
222
223 a = (3><3) [1,2 , 3
224 ,4,5 , 6
225 ,0,10,-3] :: Matrix I
226
227 b = (3><2) [0..] :: Matrix I
228
229 am = fromInt a :: Matrix (F 13)
230 bm = fromInt b :: Matrix (F 13)
231 ad = fromInt a :: Matrix Double
232 bd = fromInt b :: Matrix Double
233
234 info = do
235 print v
236 print m
237 print (tr m)
238 print $ v+v
239 print $ m+m
240 print $ m <> m
241 print $ m #> v
242
243 print $ am <> gaussElim am bm - bm
244 print $ ad <> gaussElim ad bd - bd
245
246 ok = and
247 [ toInt (m #> v) == cmod 11 (toInt m #> toInt v )
248 , am <> gaussElim am bm == bm
249 ]
250
251