$ ghc --make Vector.hs -fth
Chasing modules from: Vector.hs
[1 of 2] Skipping Fu.Prepose ( Fu/Prepose.hs, Fu/Prepose.o )
[2 of 2] Compiling Fu.Vector ( Vector.hs, Vector.o )
ghc-6.5.20051208: panic! (the `impossible' happened, GHC version 6.5.20051208):
Failed binder lookup: b{tv a1hL}
Please report it as a compiler bug to glasgow-haskell-bugs@haskell.org,
or http://sourceforge.net/projects/ghc/.
{-# OPTIONS -fglasgow-exts -fallow-undecidable-instances #-}
module Fu.Vector where
import System.IO.Unsafe
import Control.Exception
import Data.Ix
import Data.Array.IO
import qualified Data.Array as Array
import Data.Array hiding ((!))
import Fu.Prepose hiding (__)
import Prelude hiding (sum)
import Language.Haskell.TH
import Language.Haskell.TH.Syntax
import qualified Data.List
__ = __
data L n = L Int deriving (Show, Eq, Ord)
unL l@(L x) = check l x
checkBounds :: ReflectNum n => L n -> a -> a
checkBounds (l::L n) y =
if l < (L 0) || l > maxBound then
error ("Value "++show l++" out of bounds for L "++(show $ (maxBound ::
L n)))
else y
-- need Ord too?
class (Bounded a, Ix a, Eq a, Show a) => IxB a
instance (Bounded a, Ix a, Eq a, Show a) => IxB a
--instance (IxB a, IxB b) => IxB (a,b)
--instance IxB ()
instance (Bounded a, Bounded b) => Bounded (Either a b) where
minBound = Left minBound
maxBound = Right maxBound
instance (IxB a, IxB b) => Ix (Either a b) where
range (Left x, Left y) = map Left $ range (x,y)
range (Right x, Right y) = map Right $ range (x,y)
range (Left x, Right y) =
(map Left $ range (x,maxBound))++
(map Right $ range (minBound,y))
range (Right _, Left _) = error "range endpoints out of order"
index (Left x, Left y) (Left z) = index (x,y) z
index (Left x, Left y) (Right z) = error "out of bounds"
index (Right x, Right y) (Right z) = index (x,y) z
index (Right x, Right y) (Left z) = error "out of bounds"
index (Left x, Right y) (Left z) = index (x,maxBound) z
index (Left x, Right y) (Right z) =
rangeSize (x,maxBound) + index (minBound,y) z
index (Right _, Left _) _ = error "range endpoints out of order"
inRange (x,y) z = (x<=z) && (z<=y)
--instance (IxB a, IxB b) => IxB (Either a b)
domain :: IxB a => [a]
domain = range (minBound, maxBound)
instance ReflectNum n => Bounded (L n) where
minBound = L 0
maxBound = L $ reflectNum (__ :: n) - 1
instance ReflectNum n => Check (L n) where
check l = checkBounds l
class Check a where
check :: a -> (b -> b)
instance ReflectNum n => Ix (L n) where
index (a, b) (c) = index (unL a, unL b) (unL c)
range (a, b) = map L $ range (unL a, unL b)
inRange (a, b) c = inRange (unL a, unL b) (unL c)
rangeSize (a, b) = rangeSize (unL a, unL b)
-- needed?
type One = Succ Zero
-- needed?
-- singleton domain
type S = L One
-- needed?
-- empty domain
type Z = L Zero
-- like Array but uses IxB instead of Ix
newtype Vector k i = Vector (Array.Array i k) deriving (Eq, Show)
type Matrix k a b = Vector k (a,b)
boundaries :: Bounded a => (a,a)
boundaries = (minBound, maxBound)
vector :: IxB a => (a -> k) -> Vector k a
vector f = Vector
(array boundaries (map (\i -> (i, f i)) domain))
(!) :: IxB a => Vector k a -> a -> k
(!) (Vector ax) i = (Array.!) ax i
slice :: (IxB a, IxB b) => Vector k a -> (b -> a) -> Vector k b
slice v f = vector ((v!) . f)
updateArray ax i f = do
e <- readArray ax i
writeArray ax i (f e)
margin :: (IxB a, IxB b, Num k) => Vector k a -> (a -> b) -> Vector k b
margin (Vector v) f = unsafePerformIO $ (do
(ax::IOArray a k) <- newArray (minBound, maxBound) 0
mapM_ (\ (i,x) -> updateArray ax (f i) (+x)) (assocs v)
ac <- getAssocs ax
return $ Vector $ array boundaries ac
)
-- terminal map
terminal :: a -> ()
terminal x = ()
sum :: (IxB a, Num k) => Vector k a -> k
sum v = margin v (const ()) ! ()
diag :: a -> (a,a)
diag x = (x,x)
trace m = sum $ slice m diag
-- initial map
initial :: Z -> a
initial _ = undefined
(*>) :: (Num k, IxB a, IxB b, IxB c) => Matrix k a b -> Matrix k b c -> Matrix
k a c
(*>) m n =
vector (\ (i,k) -> Data.List.sum $ map (\j -> m ! (i,j) + n ! (j,k)) domain)
mapv :: (IxB a) => (k -> l) -> Vector k a -> Vector l a
mapv f v = vector (\i -> f (v!i))
mapv2 :: (IxB a) => (k -> l -> m) -> Vector k a -> Vector l a -> Vector m a
mapv2 f u v = vector (\i -> f (u!i) (v!i))
-- element-wise operations
instance (Num k, IxB a) => Num (Vector k a) where
(+) = mapv2 (+)
(-) = mapv2 (-)
(*) = mapv2 (*)
negate = mapv negate
signum = mapv signum
abs = mapv abs
fromInteger n = (fromInteger n) .* ones
-- should probably be able to do this in constant time
vec :: (Num k, IxB a, IxB b) => Matrix k a b -> Matrix k (a,b) S
vec m = slice m (\ ((i,j),_) -> (i,j))
-- better name?
inv :: (RealFrac k, IxB a) => Matrix k a a -> Matrix k a a
inv = undefined
-- need least squares operator: />?
det :: (Num k, IxB a) => Matrix k a a -> k
det = undefined
-- should reorder type args so this can be a Functor
mapm :: (IxB a, IxB b) => (k -> l) -> Matrix k a b -> Matrix l a b
-- in fact should it be a Monad? could be similar but types would need to change
mapm = mapv
vlen :: IxB a => Vector k a -> Int
vlen (v::Vector k a) = rangeSize $ (boundaries :: (a,a))
unif :: (Fractional k, IxB a) => Vector k a
unif = v
where
p = 1/(fromIntegral $ vlen v)
v = vector (const p)
ones :: (Num k, IxB a) => Vector k a
ones = vector (const 1)
zeros :: (Num k, IxB a) => Vector k a
zeros = vector (const 0)
(.*) :: (Num k, IxB a) => k -> Vector k a -> Vector k a
(.*) x v = vector (\i -> x*v!i)
normalize :: Vector k a -> Vector k a
normalize = undefined
concat :: (IxB a, IxB b) => Vector k a -> Vector k b -> Vector k (Either a b)
concat u v = vector $ \i ->
case i of
Left j -> u!j
Right k -> v!k
kronecker :: (IxB a, IxB b, Num k) => Vector k a -> Vector k b -> Vector k (a,b)
kronecker u v = vector (\(i,j) -> u!i * v!j)
-- range
(|>) :: (IxB a, IxB b) => Vector k b -> Vector b a -> Vector k a
(|>) u ix = vector (\i -> u ! (ix ! i))
(||>) :: (IxB a, IxB b, IxB c, IxB d) => Matrix k c d -> (Vector c a, Vector d
b) -> Matrix k a b
(||>) m (ix0, ix1) = vector (\ (i,j) -> m ! (ix0 ! i, ix1 ! j))
-- modify just the 'true' entries of the argument
subsetModify :: Vector k b -> (b -> Bool) -> (forall a . Vector k a -> Vector k
a) -> Vector k b
subsetModify = undefined
filterModify :: Vector k b -> (k -> b -> Bool) -> (forall a . Vector k a ->
Vector k a) -> Vector k b
filterModify = undefined
-- need to inject values into a vector...
-- permutations, size is n!
data Perm n = Perm (Vector (L n) (L n))
--det :: (Num k, IxB a) => Matrix a a k -> k
listVec :: [k] -> (forall a . IxB a => Vector k a -> w) -> w
listVec (l::[k]) f =
let n = length l in
reifyIntegral n (\ (_::rn) ->
f (listToVec l::Vector k (L rn)))
vecQ :: Lift k => [k] -> Q Exp
vecQ (l::[k]) =
let n = length l in
reifyIntegral n (\ (x::rn) ->
let r = T x
rt = [| r |] in
[| let (_::b) = $(rt) in listToVec l::Vector k (L b) |])
listToVec :: IxB a => [k] -> Vector k a
listToVec l::Vector k a = (v::Vector k a)
where
v = assert (length l == vlen v) $
Vector $ listArray (boundaries :: (a,a)) l
--listMat :: [[k]] -> (forall a b . (IxB a, IxB b) => Matrix k (a,b) -> w) -> w
main = undefined
{-# OPTIONS -fglasgow-exts -fallow-undecidable-instances #-}
-- Big comment:
-- Note that the terms "reify" and "reflect" may seem to be used in a way
-- which is the reverse of what one might intuit! Here, the
-- representation of a number in type-level binary is considered "real".
-- Thus to turn a value into a type-level term is to "reify" it; to go
-- the other way is to "reflect".
module Fu.Prepose where
import System.IO.Unsafe (unsafePerformIO)
import Control.Exception (handle, handleJust, errorCalls)
import Foreign.Marshal.Utils (with, new)
import Foreign.Marshal.Array (peekArray, pokeArray)
import Foreign.Marshal.Alloc (alloca)
import Foreign.Ptr (Ptr, castPtr)
import Foreign.Storable (Storable(sizeOf, peek))
import Foreign.C.Types (CChar)
import Foreign.StablePtr (StablePtr, newStablePtr,
deRefStablePtr, freeStablePtr)
import Data.Bits (Bits(..))
import Prelude hiding (getLine)
import Language.Haskell.TH hiding (reify)
import Language.Haskell.TH.Syntax hiding (reify)
import Debug.Trace
newtype Modulus s a = Modulus a deriving (Eq, Show)
newtype M s a = M a deriving (Eq, Show)
unM :: M s a -> a
unM (M a) = a
__ = __
class Modular s a | s -> a where modulus :: s -> a
normalize :: (Modular s a, Integral a) => a -> M s a
normalize a :: M s a = M (mod a (modulus (__ :: s)))
instance (Modular s a, Integral a) => Num (M s a) where
M a + M b = normalize (a + b)
M a - M b = normalize (a - b)
M a * M b = normalize (a * b)
negate (M a) = normalize (negate a)
fromInteger i = normalize (fromInteger i)
signum = error "Modular numbers are not signed"
abs = error "Modular numbers are not signed"
aaa (M a :: M s1 a) (M b :: M s2 a) = normalize (a + b) :: M s1 a
where _ = [__::s1, __ :: s2] -- make sure input moduli are equal
withModulus :: a -> (forall s. Modular s a => s -> w) -> w
data Zero; data Twice s; data Succ s; data Pred s
data Positive s; data Negative s; data TwiceSucc s
----------------------------------------------------------------
class GetType s where
getType :: s -> Type
instance GetType (Zero) where
getType _ = ConT $ mkName "Zero"
instance GetType s => GetType (Twice s) where
getType _ = AppT (ConT $ mkName "Twice") (getType (__ :: s))
instance GetType s => GetType (Succ s) where
getType _ = AppT (ConT $ mkName "Succ") (getType (__ :: s))
instance GetType s => GetType (Pred s) where
getType _ = AppT (ConT $ mkName "Pred") (getType (__ :: s))
instance GetType s => GetType (Positive s) where
getType _ = AppT (ConT $ mkName "Positive") (getType (__ :: s))
instance GetType s => GetType (Negative s) where
getType _ = AppT (ConT $ mkName "Negative") (getType (__ :: s))
instance GetType s => GetType (TwiceSucc s) where
getType _ = AppT (ConT $ mkName "TwiceSucc") (getType (__ :: s))
data T s = T s
instance GetType s => Lift (T s) where
lift _ = return $ SigE undE (getType (__ :: s))
undE = VarE $ mkName "undefined"
----------------------------------------------------------------
class GetType s => ReflectUnsigned s where reflectUnsigned :: Num a => s -> a
instance ReflectUnsigned Zero where
reflectUnsigned _ = 0
instance ReflectUnsigned s => ReflectUnsigned (Twice s) where
reflectUnsigned _ = reflectUnsigned (__ :: s) * 2
instance ReflectUnsigned s => ReflectUnsigned (TwiceSucc s) where
reflectUnsigned _ = reflectUnsigned (__ :: s) * 2 + 1
instance ReflectUnsigned s => ReflectNum (Positive s) where
reflectNum _ = reflectUnsigned (__ :: s)
instance ReflectUnsigned s => ReflectNum (Negative s) where
reflectNum _ = -1 - reflectUnsigned (__ :: s)
class GetType s => ReflectNum s where reflectNum :: Num a => s -> a
instance ReflectNum Zero where
reflectNum _ = 0
instance ReflectNum s => ReflectNum (Twice s) where
reflectNum _ = reflectNum (__ :: s) * 2
instance ReflectNum s => ReflectNum (Succ s) where
reflectNum _ = reflectNum (__ :: s) + 1
instance ReflectNum s => ReflectNum (Pred s) where
reflectNum _ = reflectNum (__ :: s) - 1
reifyIntegral :: Integral a =>
a -> (forall s. ReflectNum s => s -> w) -> w
reifyIntegral i k = --(trace $ "reifyIntegral "++show i) $
case quotRem i 2 of
(0, 0) -> k (__ :: Zero)
(0, 1) -> k (__ :: Succ Zero)
(0,- 1) -> k (__ :: Pred Zero)
(j, 0) -> reifyIntegral j (\ (_ :: s) -> k (__ :: Twice s))
(j, 1) -> reifyIntegral j (\ (_ :: s) -> k (__ :: Succ (Twice s)))
(j,- 1) -> reifyIntegral j (\ (_ :: s) -> k (__ :: Pred (Twice s)))
data ModulusNum s a
instance (ReflectNum s, Num a) =>
Modular (ModulusNum s a) a where
modulus _ = reflectNum (__ :: s)
withIntegralModulus :: Integral a =>
a -> (forall s. Modular s a => s -> w) -> w
withIntegralModulus i k =
reifyIntegral i (\(_ :: s) -> k (__ :: ModulusNum s a))
data Nil; data Cons s ss
class ReflectNums ss where reflectNums :: Num a => ss -> [a]
instance ReflectNums Nil where
reflectNums _ = []
instance (ReflectNum s, ReflectNums ss) =>
ReflectNums (Cons s ss) where
reflectNums _ = reflectNum (__ :: s) : reflectNums (__ :: ss)
reifyIntegrals :: Integral a =>
[a] -> (forall ss. ReflectNums ss => ss -> w) -> w
reifyIntegrals [] k = k (__ :: Nil)
reifyIntegrals (i:ii) k = reifyIntegral i (\(_ :: s) ->
reifyIntegrals ii (\(_ :: ss) ->
k (__ :: Cons s ss)))
type Byte = CChar
data Store s a
class ReflectStorable s where
reflectStorable :: Storable a => s a -> a
instance ReflectNums s => ReflectStorable (Store s) where
reflectStorable _ = unsafePerformIO
$ alloca
$ \p -> do pokeArray (castPtr p) bytes
peek p
where bytes = reflectNums (__ :: s) :: [Byte]
reifyStorable :: Storable a =>
a -> (forall s. ReflectStorable s => s a -> w) -> w
reifyStorable a k =
reifyIntegrals (bytes :: [Byte]) (\(_ :: s) -> k (__ :: Store s a))
where bytes = unsafePerformIO
$ with a (peekArray (sizeOf a) . castPtr)
class Reflect s a | s -> a where reflect :: s -> a
data Stable (s :: * -> *) a
instance ReflectStorable s => Reflect (Stable s a) a where
reflect = unsafePerformIO
$ do a <- deRefStablePtr p
return (const a)
where p = reflectStorable (__ :: s p)
reify :: a -> (forall s. Reflect s a => s -> w) -> w
reify (a :: a) k = unsafePerformIO
$ do p <- newStablePtr a
reifyStorable p (\(_ :: s p) ->
k' (__ :: Stable s a))
where k' s = return (k s)
data ModulusAny s
instance Reflect s a => Modular (ModulusAny s) a where
modulus _ = reflect (__ :: s)
withModulus a k = reify a (\(_ :: s) -> k (__ :: ModulusAny s))
withIntegralModulus' :: Integral a =>
a -> (forall s. Modular s a => M s w) -> w
withIntegralModulus' (i::a) k :: w =
reifyIntegral i ( \(_ :: t) ->
unM (k :: M (ModulusNum t a) w))
test4' :: (Modular s a, Integral a) => M s a
test4' = 3*3 + 5*5
test4 = withIntegralModulus' 4 test4'
data Even p q u v a = E a a deriving (Eq, Show)
normalizeEven :: ( ReflectNum p, ReflectNum q, Integral a,
Bits a) => a -> a -> Even p q u v a
normalizeEven a b :: Even p q u v a =
E (a .&. (shiftL 1 (reflectNum (__ :: p)) - 1)) -- $a \bmod 2^p$
(mod b (reflectNum (__ :: q))) -- $b \bmod q$
instance ( ReflectNum p, ReflectNum q,
ReflectNum u, ReflectNum v,
Integral a, Bits a) => Num (Even p q u v a) where
E a1 b1 + E a2 b2 = normalizeEven (a1 + a2) (b1 + b2)
{-"\raisebox{0pt}[2.5ex][0pt]{$\vdots$}"-}
E a1 b1 - E a2 b2 = normalizeEven (a1 - a2) (b1 - b2)
E a1 b1 * E a2 b2 = normalizeEven (a1 * a2) (b1 * b2)
negate (E a b) = normalizeEven (-a) (-b)
fromInteger i = normalizeEven (fromInteger i)
(fromInteger i)
signum = error "Modular numbers are not signed"
abs = error "Modular numbers are not signed"
gcd' x 0 = []
gcd' x y = (x `div` y) : (gcd' y (x `mod` y))
-- chinese remainder theorem
crt a b =
let (c,d) = snd $ foldl (\ ((a,b),(a',b')) p -> ((-(a*p+a'),b*p+b'),(-a,b)))
((1,0),(0,1)) (gcd' a b)
in (-c*b,d*a)
withIntegralModulus'' :: (Integral a, Bits a) =>
a -> (forall s. Num (s a) => s a) -> a
withIntegralModulus'' (i::a) k = case factor 0 i of
(0, i) -> withIntegralModulus' i k -- odd case
(p, q) -> let (u, v) = crt (2^p) q in -- even case: $i =
2^p q$
-- trace ("(u,v)="++show(u,v)) $
reifyIntegral p (\(_::p ) ->
reifyIntegral q (\(_::q ) ->
reifyIntegral u (\(_::u ) ->
reifyIntegral v (\(_::v ) ->
unEven (k :: Even p q u v a)))))
factor :: (Num p, Integral q) => p -> q -> (p, q)
factor p i = case quotRem i 2 of
(0, 0) -> (0, 0) -- just zero
(j, 0) -> factor (p+1) j -- accumulate powers of two
_ -> (p, i) -- not even
unEven ::( ReflectNum p, ReflectNum q, ReflectNum u,
ReflectNum v, Integral a, Bits a) => Even p q u v a -> a
unEven (E a b :: Even p q u v a) = -- trace "unEven" $
mod (a * (reflectNum (__ :: u)) + b * (reflectNum (__ :: v)))
(shiftL (reflectNum (__ :: q)) (reflectNum (__ :: p)))
test5 :: Num (s a) => s a
test5 = 1000*1000 + 513*513
test5' = withIntegralModulus'' 1279 test5 :: Integer
test5'' = withIntegralModulus'' 1280 test5 :: Integer
_______________________________________________
Glasgow-haskell-bugs mailing list
Glasgow-haskell-bugs@haskell.org
http://www.haskell.org/mailman/listinfo/glasgow-haskell-bugs
