module Math.Sym.Internal
(
Lehmercode
, Perm0
, unrankLehmercode
, fromLehmercode
, randomLehmercode
, lehmercodes
, size
, toList
, fromList
, act
, unrankPerm
, randomPerm
, sym
, idperm
, revIdperm
, sti
, st
, ordiso
, simple
, copies
, avoiders
, reverse
, complement
, inverse
, rotate
, asc
, des
, exc
, fp
, inv
, maj
, peak
, vall
, dasc
, ddes
, lmin
, lmax
, rmin
, rmax
, head
, last
, lir
, ldr
, rir
, rdr
, comp
, ep
, dim
, stackSort
, bubbleSort
, del
, onesCUInt
, nextCUInt
, nextIntegral
) where
import Prelude hiding (reverse, head, last)
import qualified Prelude (head)
import System.Random (getStdRandom, randomR)
import Control.Monad (forM_, liftM)
import Control.Monad.ST (runST)
import Data.List (group)
import Data.Bits (Bits, shiftR, (.|.), (.&.), popCount)
import qualified Data.Vector.Storable as SV
( Vector, toList, fromList, length, (!), thaw
, unsafeFreeze, unsafeWith, enumFromN, enumFromStepN
, head, last, filter, maximum, minimum, null, reverse, map
)
import qualified Data.Vector.Storable.Mutable as MV
( unsafeNew, unsafeWrite, unsafeWith, swap, replicate
)
import Foreign (Ptr, castPtr)
import System.IO.Unsafe (unsafePerformIO)
import Foreign.C.Types (CLong(..), CInt(..), CUInt(..))
import Foreign.Marshal.Utils (toBool)
type Lehmercode = SV.Vector Int
type Perm0 = SV.Vector Int
unrankLehmercode :: Int -> Integer -> Lehmercode
unrankLehmercode n rank = runST $ do
v <- MV.unsafeNew n
iter v n rank (toInteger n)
SV.unsafeFreeze v
where
iter _ 0 _ _ = return ()
iter v i r m = do
let (r',j) = quotRem r m
MV.unsafeWrite v (ni) (fromIntegral j)
iter v (i1) r' (m1)
fromLehmercode :: Lehmercode -> Perm0
fromLehmercode code = runST $ do
let n = SV.length code
v <- MV.unsafeNew n
forM_ [0..n1] $ \i -> MV.unsafeWrite v i i
forM_ [0..n1] $ \i -> MV.swap v i (i + (SV.!) code i)
SV.unsafeFreeze v
randomLehmercode :: Int -> IO Lehmercode
randomLehmercode n = unrankLehmercode n `liftM` getStdRandom (randomR (0, factorial n 1))
lehmercodes :: Int -> [Lehmercode]
lehmercodes n = map (unrankLehmercode n) [0 .. factorial n 1]
size :: Perm0 -> Int
size = SV.length
toList :: Perm0 -> [Int]
toList = SV.toList
fromList :: [Int] -> Perm0
fromList = SV.fromList
act :: Perm0 -> Perm0 -> Perm0
act u v = runST $ do
let n = SV.length u
w <- MV.unsafeNew n
forM_ [0..n1] $ \i -> MV.unsafeWrite w i ((SV.!) v ((SV.!) u i))
SV.unsafeFreeze w
factorial :: Integral a => a -> Integer
factorial = product . enumFromTo 1 . toInteger
unrankPerm :: Int -> Integer -> Perm0
unrankPerm n = fromLehmercode . unrankLehmercode n
randomPerm :: Int -> IO Perm0
randomPerm n = fromLehmercode `liftM` randomLehmercode n
sym :: Int -> [Perm0]
sym n = map (unrankPerm n) [0 .. factorial n 1]
idperm :: Int -> Perm0
idperm = SV.enumFromN 0
revIdperm :: Int -> Perm0
revIdperm n = SV.enumFromStepN (n1) (1) n
sti :: Perm0 -> Perm0
sti w = runST $ do
let a = if SV.null w then 0 else SV.minimum w
let b = if SV.null w then 0 else SV.maximum w
let n = SV.length w
v <- MV.replicate (1 + b a) (1)
forM_ [0..n1] $ \i -> MV.unsafeWrite v ((SV.!) w i a) i
SV.filter (>=0) `liftM` SV.unsafeFreeze v
st :: Perm0 -> Perm0
st = inverse . sti
foreign import ccall unsafe "ordiso.h ordiso" c_ordiso
:: Ptr CLong -> Ptr CLong -> Ptr CLong -> CLong -> CInt
ordiso :: Perm0 -> Perm0 -> SV.Vector Int -> Bool
ordiso u v m =
let k = fromIntegral (SV.length u)
in unsafePerformIO $
SV.unsafeWith u $ \u' ->
SV.unsafeWith v $ \v' ->
SV.unsafeWith m $ \m' ->
return . toBool $ c_ordiso (castPtr u') (castPtr v') (castPtr m') k
foreign import ccall unsafe "simple.h simple" c_simple
:: Ptr CLong -> CLong -> CInt
simple :: Perm0 -> Bool
simple w =
let n = fromIntegral (SV.length w)
in unsafePerformIO $
SV.unsafeWith w $ \w' ->
return . toBool $ c_simple (castPtr w') n
copies :: (Int -> Int -> [SV.Vector Int]) -> Perm0 -> Perm0 -> [SV.Vector Int]
copies subsets p w = filter (ordiso p w) $ subsets n k
where
n = SV.length w
k = SV.length p
avoiders1 :: (Int -> Int -> [SV.Vector Int]) -> (a -> Perm0) -> Perm0 -> [a] -> [a]
avoiders1 subsets f p ws =
let ws0 = map f ws
ws2 = zip ws0 ws
in case group (map SV.length ws0) of
[] -> []
[_] -> let k = SV.length p
n = SV.length (Prelude.head ws0)
in [ v | (v0,v) <- ws2, not $ any (ordiso p v0) (subsets n k) ]
_ -> [ v | (v0,v) <- ws2, null $ copies subsets p v0 ]
avoiders :: (Int -> Int -> [SV.Vector Int]) -> (a -> Perm0) -> [Perm0] -> [a] -> [a]
avoiders _ _ [] ws = ws
avoiders subsets f (p:ps) ws = avoiders subsets f ps $ avoiders1 subsets f p ws
reverse :: Perm0 -> Perm0
reverse = SV.reverse
complement :: Perm0 -> Perm0
complement w = SV.map (\x -> SV.length w x 1) w
inverse :: Perm0 -> Perm0
inverse w = runST $ do
let n = SV.length w
v <- MV.unsafeNew n
forM_ [0..n1] $ \i -> MV.unsafeWrite v ((SV.!) w i) i
SV.unsafeFreeze v
rotate :: Perm0 -> Perm0
rotate w = runST $ do
let n = SV.length w
v <- MV.unsafeNew n
forM_ [0..n1] $ \i -> MV.unsafeWrite v ((SV.!) w (n1i)) i
SV.unsafeFreeze v
foreign import ccall unsafe "stat.h asc" c_asc
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h des" c_des
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h exc" c_exc
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h fp" c_fp
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h inv" c_inv
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h maj" c_maj
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h peak" c_peak
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h vall" c_vall
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h dasc" c_dasc
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h ddes" c_ddes
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h lmin" c_lmin
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h lmax" c_lmax
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h lir" c_lir
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h ldr" c_ldr
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h comp" c_comp
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h ep" c_ep
:: Ptr CLong -> CLong -> CLong
foreign import ccall unsafe "stat.h dim" c_dim
:: Ptr CLong -> CLong -> CLong
stat :: (Ptr CLong -> CLong -> CLong) -> Perm0 -> Int
stat f w = unsafePerformIO $
SV.unsafeWith w $ \ptr ->
return . fromIntegral $ f (castPtr ptr) (fromIntegral (SV.length w))
head :: Perm0 -> Int
head = SV.head
last :: Perm0 -> Int
last = SV.last
rmin :: Perm0 -> Int
rmin = lmin . SV.reverse
rmax :: Perm0 -> Int
rmax = lmax . SV.reverse
rir :: Perm0 -> Int
rir = ldr . SV.reverse
rdr :: Perm0 -> Int
rdr = lir . SV.reverse
asc :: Perm0 -> Int
asc = stat c_asc
des :: Perm0 -> Int
des = stat c_des
inv :: Perm0 -> Int
inv = stat c_inv
maj :: Perm0 -> Int
maj = stat c_maj
peak :: Perm0 -> Int
peak = stat c_peak
vall :: Perm0 -> Int
vall = stat c_vall
dasc :: Perm0 -> Int
dasc = stat c_dasc
ddes :: Perm0 -> Int
ddes = stat c_ddes
lmin :: Perm0 -> Int
lmin = stat c_lmin
lmax :: Perm0 -> Int
lmax = stat c_lmax
lir :: Perm0 -> Int
lir = stat c_lir
ldr :: Perm0 -> Int
ldr = stat c_ldr
exc :: Perm0 -> Int
exc = stat c_exc
fp :: Perm0 -> Int
fp = stat c_fp
comp :: Perm0 -> Int
comp = stat c_comp
ep :: Perm0 -> Int
ep = stat c_ep
dim :: Perm0 -> Int
dim = stat c_dim
foreign import ccall unsafe "sortop.h stacksort" c_stacksort
:: Ptr CLong -> CLong -> IO ()
foreign import ccall unsafe "sortop.h bubblesort" c_bubblesort
:: Ptr CLong -> CLong -> IO ()
sortop :: (Ptr CLong -> CLong -> IO ()) -> Perm0 -> Perm0
sortop f w = unsafePerformIO $ do
v <- SV.thaw w
MV.unsafeWith v $ \ptr -> do
f (castPtr ptr) (fromIntegral (SV.length w))
SV.unsafeFreeze v
stackSort :: Perm0 -> Perm0
stackSort = sortop c_stacksort
bubbleSort :: Perm0 -> Perm0
bubbleSort = sortop c_bubblesort
del :: Int -> Perm0 -> Perm0
del i u = runST $ do
let n = SV.length u
let j = (SV.!) u i
v <- MV.unsafeNew (n1)
forM_ [0..i1] $ \k -> do
let m = (SV.!) u k
MV.unsafeWrite v k (if m < j then m else m1)
forM_ [i+1..n1] $ \k -> do
let m = (SV.!) u k
MV.unsafeWrite v (k1) (if m < j then m else m1)
SV.unsafeFreeze v
foreign import ccall unsafe "bit.h next" c_next :: CUInt -> CUInt
nextCUInt :: CUInt -> CUInt
nextCUInt = c_next
foreign import ccall unsafe "bit.h ones" c_ones :: Ptr CUInt -> CUInt -> IO ()
onesCUInt :: CUInt -> SV.Vector Int
onesCUInt m = SV.map fromIntegral . unsafePerformIO $ do
v <- MV.unsafeNew (popCount m)
MV.unsafeWith v $ \ptr -> do
c_ones ptr m
SV.unsafeFreeze v
nextIntegral :: (Integral a, Bits a) => a -> a
nextIntegral a =
let b = (a .|. (a 1)) + 1
in b .|. ((((b .&. (b)) `div` (a .&. (a))) `shiftR` 1) 1)