diff options
Diffstat (limited to 'packages/base/src/Numeric')
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Util/Modular.hs | 251 |
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 | {- | | ||
15 | Module : Numeric.LinearAlgebra.Util.Modular | ||
16 | Copyright : (c) Alberto Ruiz 2015 | ||
17 | License : BSD3 | ||
18 | Stability : experimental | ||
19 | |||
20 | Proof of concept of statically checked modular arithmetic. | ||
21 | |||
22 | -} | ||
23 | |||
24 | module Numeric.LinearAlgebra.Util.Modular( | ||
25 | Mod, F | ||
26 | ) where | ||
27 | |||
28 | import Data.Packed.Numeric | ||
29 | import Numeric.LinearAlgebra.Util(Indexable(..),gaussElim) | ||
30 | import GHC.TypeLits | ||
31 | import Data.Proxy(Proxy) | ||
32 | import Foreign.ForeignPtr(castForeignPtr) | ||
33 | import Data.Vector.Storable(unsafeToForeignPtr, unsafeFromForeignPtr) | ||
34 | import Foreign.Storable | ||
35 | import Data.Ratio | ||
36 | import Data.Packed.Internal.Matrix hiding (mat,size) | ||
37 | import Data.Packed.Internal.Numeric | ||
38 | |||
39 | |||
40 | -- | Wrapper with a phantom integer for statically checked modular arithmetic. | ||
41 | newtype Mod (n :: Nat) t = Mod {unMod:: t} | ||
42 | deriving (Storable) | ||
43 | |||
44 | instance KnownNat m => Enum (F m) | ||
45 | where | ||
46 | toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m)) | ||
47 | fromEnum = fromIntegral . unMod | ||
48 | |||
49 | instance KnownNat m => Eq (F m) | ||
50 | where | ||
51 | a == b = (unMod a) == (unMod b) | ||
52 | |||
53 | instance KnownNat m => Ord (F m) | ||
54 | where | ||
55 | compare a b = compare (unMod a) (unMod b) | ||
56 | |||
57 | instance KnownNat m => Real (F m) | ||
58 | where | ||
59 | toRational x = toInteger x % 1 | ||
60 | |||
61 | instance 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 | ||
69 | instance 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 | |||
79 | l2 :: forall m a b c. (KnownNat m) => (Int -> a -> b -> c) -> Mod m a -> Mod m b -> Mod m c | ||
80 | l2 f (Mod u) (Mod v) = Mod (f m' u v) | ||
81 | where | ||
82 | m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int | ||
83 | |||
84 | l1 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> Mod m a -> Mod m b | ||
85 | l1 f (Mod u) = Mod (f m' u) | ||
86 | where | ||
87 | m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int | ||
88 | |||
89 | l0 :: forall m a b . (KnownNat m) => (Int -> a -> b) -> a -> Mod m b | ||
90 | l0 f u = Mod (f m' u) | ||
91 | where | ||
92 | m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int | ||
93 | |||
94 | |||
95 | instance Show (F n) | ||
96 | where | ||
97 | show = show . unMod | ||
98 | |||
99 | instance 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 | ||
110 | type F n = Mod n I | ||
111 | |||
112 | type V n = Vector (F n) | ||
113 | type M n = Matrix (F n) | ||
114 | |||
115 | |||
116 | instance 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 | |||
127 | instance 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 | |||
165 | instance Indexable (Vector (F m)) (F m) | ||
166 | where | ||
167 | (!) = (@>) | ||
168 | |||
169 | |||
170 | type instance RealOf (F n) = I | ||
171 | |||
172 | |||
173 | instance KnownNat m => Product (F m) where | ||
174 | norm2 = undefined | ||
175 | absSum = undefined | ||
176 | norm1 = undefined | ||
177 | normInf = undefined | ||
178 | multiply = lift2 multiply | ||
179 | |||
180 | |||
181 | instance KnownNat m => Numeric (F m) | ||
182 | |||
183 | i2f :: Vector I -> Vector (F n) | ||
184 | i2f v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n) | ||
185 | where (fp,i,n) = unsafeToForeignPtr v | ||
186 | |||
187 | f2i :: Vector (F n) -> Vector I | ||
188 | f2i v = unsafeFromForeignPtr (castForeignPtr fp) (i) (n) | ||
189 | where (fp,i,n) = unsafeToForeignPtr v | ||
190 | |||
191 | f2iM :: Matrix (F n) -> Matrix I | ||
192 | f2iM = liftMatrix f2i | ||
193 | |||
194 | i2fM :: Matrix I -> Matrix (F n) | ||
195 | i2fM = liftMatrix i2f | ||
196 | |||
197 | |||
198 | lift1 f a = fromInt (f (toInt a)) | ||
199 | lift2 f a b = fromInt (f (toInt a) (toInt b)) | ||
200 | |||
201 | instance 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 | |||
214 | instance (KnownNat m) => Testable (M m) | ||
215 | where | ||
216 | checkT _ = test | ||
217 | |||
218 | test = (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 | |||