diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,5 @@ +# Revision history for variety + +## 0.1.0.0 -- 2025-06-04 + +* First version. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,20 @@ +Copyright (c) 2025 nbos + +Permission is hereby granted, free of charge, to any person obtaining +a copy of this software and associated documentation files (the +"Software"), to deal in the Software without restriction, including +without limitation the rights to use, copy, modify, merge, publish, +distribute, sublicense, and/or sell copies of the Software, and to +permit persons to whom the Software is furnished to do so, subject to +the following conditions: + +The above copyright notice and this permission notice shall be included +in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. diff --git a/src/Codec/Arithmetic/Combinatorics.hs b/src/Codec/Arithmetic/Combinatorics.hs new file mode 100644 --- /dev/null +++ b/src/Codec/Arithmetic/Combinatorics.hs @@ -0,0 +1,282 @@ +{-# LANGUAGE BangPatterns #-} +-- | Optimal codes for combinatorial objects. +-- +-- The integer on which a combinatorial objects is mapped is typically +-- called its rank. Below are implementations of ranking and unranking +-- algorithms for the indexes of common combinatorial objects in the +-- lexicographic enumeration of objects of the same parameters. +module Codec.Arithmetic.Combinatorics + ( -- * Multiset Permutations + + -- | [Multiset permutations] + -- (https://en.wikipedia.org/wiki/Permutation#Permutations_of_multisets) + -- are ways to order the elements of a set where elements may appear + -- more than once. The number of such permutations is equal to the + -- multinomial coefficient with the same parameters: \[ {n \choose + -- k_{1}, k_{2}, \ldots, k_{m}} = \frac{n!}{k_{1}! k_{2}! \cdots + -- k_{m}!} \] This is the most general definition in this module, + -- of which all following objects are special cases. + + rankMultisetPermutation + , unrankMultisetPermutation + , multinomial + + -- * Permutations + + -- | A [permutation](https://en.wikipedia.org/wiki/Permutation) is an + -- ordering of all the objects of a set. The number of permutations of + -- a set of \(n\) elements is \(n!\). + + , rankPermutation + , unrankPermutation + + -- * Combinations + + -- | A [combination](https://en.wikipedia.org/wiki/Combination) is a + -- selection of \(k\) elements from a set of size \(n\). The number of + -- combinations for parameters \(n\) and \(k\) is given by the + -- binomial coefficient: \[ {n \choose k} = \frac{n!}{k! (n-k)!} \] + + , rankCombination + , unrankCombination + , choose + + -- * Distributions + + -- | A distribution (usually discussed under the name [stars and + -- bars](https://en.wikipedia.org/wiki/Stars_and_bars_(combinatorics\)) + -- is a way to distribute \(n\) equal elements (stars) among \(k\) + -- bins (i.e. \(k-1\) bars ). + + , rankDistribution + , unrankDistribution + + -- * Non-Empty Distributions + + -- | The class of distributions that have at least one element per + -- bin. + + , rankDistribution1 + , unrankDistribution1 + ) where + +import Control.Exception (assert) +import Data.Maybe (fromJust) +import qualified Data.Set as S +import qualified Data.List as L +import Data.Map.Strict (Map) +import qualified Data.Map.Strict as M +import Math.Combinatorics.Exact.Factorial (factorial) + +import qualified Codec.Arithmetic.Variety as V + +err :: String -> a +err = error . ("Combinatorics." ++) + +-- | Rank a multiset permutation. Returns the count of each element in +-- the set, the rank and the total number of permutations with those +-- counts (the multinomial coefficient). +rankMultisetPermutation :: Ord a => [a] -> ([(a,Int)], (Integer, Integer)) +rankMultisetPermutation msp = ( M.toList counts + , (index, coef0) ) + where + counts = L.foldl' (\m k -> M.insertWith (+) k 1 m) M.empty msp + total0 = sum counts + coef0 = factorial total0 + `div` product (factorial <$> counts) + index = sum $ go (fromIntegral total0) coef0 counts msp + + go :: Ord a => Integer -> Integer -> Map a Int -> [a] -> [Integer] + go _ _ _ [] = [] + go total coef m (a:as) = sum lowerSubCoefs : + go total' coef' m' as + where + (lt,eq,_) = M.splitLookup a m + total' = total - 1 -- decrement `total` by 1 + lowerSubCoefs = (`div` total) . (coef *) . fromIntegral <$> lt + n = fromJust eq + n' = n - 1 -- decrement `a`'s count by 1 + coef' = (coef * fromIntegral n) `div` total -- rm `n` factor from denom + m' = M.update (\_ -> if n' == 0 then Nothing else Just n') + a m + +-- | Reconstruct a multiset permutation, given the count of each element +-- in the set and a rank. +unrankMultisetPermutation :: Ord a => [(a,Int)] -> Integer -> [a] +unrankMultisetPermutation l i0 + | any ((< 0) . snd) l = err' "negative count" + | i0 < 0 || i0 >= coef0 = err' $ "out of bounds: " ++ show (i0,coef0) + | otherwise = go (fromIntegral total0) coef0 counts i0 + where + err' = err . ("unrankMultisetPermutation: " ++) + counts = M.fromList $ filter ((> 0) . snd) l + total0 = sum counts + coef0 = factorial total0 + `div` product (factorial <$> counts) + + go total coef m i | M.null m = [] + | otherwise = a : go total' coef' m' i' + where + total' = total - 1 -- decrement `total` by 1 + subCoefs = (`div` total) . (coef *) . fromIntegral <$> m + (a, lowerSubCoefsSum, coef') = findBin 0 $ M.toList subCoefs + i' = i - lowerSubCoefsSum -- update index to local bin + m' = M.update (\n -> if n == 1 then Nothing else Just $ n - 1) + a m + + findBin _ [] = err "impossible" + findBin acc ((el,subCoef):ascs) + | null ascs || acc' > i = (el, acc, subCoef) + | otherwise = findBin acc' ascs + where acc' = acc + subCoef + +-- | Computes the multinomial coefficient. +multinomial :: [Int] -> Integer +multinomial ns | any (< 0) ns = 0 + | otherwise = factorial (sum ns) + `div` product (factorial <$> ns) + +-- | Rank a permutation. Returns the rank and the total number of +-- permutations of sets with that size ( \(n!\) ). +rankPermutation :: Ord a => [a] -> (Integer, Integer) +rankPermutation p | length p /= n0 = err' "not unique elements" + | otherwise = V.fromValue val + where + err' = err . ("rankPermutation: " ++) + s0 = S.fromList p + n0 = S.size s0 + ns = fromIntegral <$> [n0,n0-1..1] + is = fromIntegral <$> go s0 p + val = assert (length is == length ns) + mconcat $ + zipWith V.mkValue is ns + + -- | Lookup element index in the set of remaining elements + go s [] = assert (S.null s) [] + go s (a:rest) = i : go s' rest + where i = S.findIndex a s + s' = S.delete a s + +-- | Reconstruct a permutation given a set of elements and a rank. The +-- order in which the elements of the set is given does not matter. +unrankPermutation :: Ord a => [a] -> Integer -> [a] +unrankPermutation as index + | length as /= n = err' "not unique elements" + | index < 0 || index >= base = err' $ "out of bounds" ++ show (index,base) + | otherwise = go set is + where + err' = err . ("unrankPermutation: " ++) + set = S.fromList as + n = S.size set + ns = fromIntegral <$> [n,n-1..1] + base = factorial $ fromIntegral n + bv = V.toBitVec $ V.mkValue index base + is = fromIntegral <$> V.decode ns bv + + -- | Successively delete elements at given indexes from a set + go s [] = assert (S.null s) [] + go s (i:rest) = S.elemAt i s : go (S.deleteAt i s) rest + +-- | Rank a combination in the form of a list of booleans. Returns the +-- \((n,k)\) parameters (where \(k\) is the number of `True` values and +-- \(n\) is the total), the rank and the total number of combinations +-- with those parameters (the binomial coefficient). +rankCombination :: [Bool] -> ((Int, Int), (Integer, Integer)) +rankCombination c = ( (n0, k0) + , (res, n0Ck0) ) + where + n0 = length c + k0 = sum $ fromEnum <$> c + n0Ck0 = n0 `choose` k0 + res = sum $ go (fromIntegral n0) (fromIntegral k0) n0Ck0 c + + go :: Integer -> Integer -> Integer -> [Bool] -> [Integer] + go _ _ _ [] = [] + go n k nCk (b:bs) = if b then nCk0 : go (n-1) (k-1) nCk1 bs + else go (n-1) k nCk0 bs + where + nCk0 = nCk - nCk1 -- sub coef if 0/False + nCk1 = (nCk * k) `div` n -- sub coef if 1/True + +-- | Reconstruct a combination given parameters \((n,k)\) and a rank. +unrankCombination :: (Int, Int) -> Integer -> [Bool] +unrankCombination nk@(n0,k0) i0 + | k0 > n0 || k0 < 0 || n0 < 0 = err' $ "invalid parameters: " ++ show nk + | i0 < 0 || i0 > n0Ck0 = err' $ "out of range: " ++ show (i0,n0Ck0) + | otherwise = go (fromIntegral n0) (fromIntegral k0) n0Ck0 i0 + + where + err' = err . ("unrankPermutation: " ++) + n0Ck0 = n0 `choose` k0 + go n k nCk i | n == 0 = [] + | i < nCk0 = False : go (n-1) k nCk0 i + | otherwise = True : go (n-1) (k-1) nCk1 (i-nCk0) + where + nCk0 = nCk - nCk1 -- sub coef if 0/False + nCk1 = (nCk * k) `div` n -- sub coef if 1/True + +-- | Computes the binomial coefficent given parameters \(n\) and \(k\). +choose :: Int -> Int -> Integer +choose n k | denom == 0 = 0 + | otherwise = num `div` denom + where num = factorial n + denom = factorial k * factorial (n-k) + +-- | Rank a distribution in the form of a list bin counts. Returns the +-- \((n,k)\) parameters (where \(n\) is the total number of elements and +-- \(k\) is the number of bins), the rank and the total number of +-- distributions with those parameters. +rankDistribution :: [Int] -> ((Int, Int), (Integer, Integer)) +rankDistribution [] = ((0,0),(0,1)) +rankDistribution (n0:ns) + | n0 < 0 || any (< 0) ns = err' "negative count" + | otherwise = ((bins,balls),(i,base)) + where + err' = err . ("rankDistribution: " ++) + comb = replicate n0 False -- 0s are stars, 1s are bars + ++ concatMap ((True:) . flip replicate False) ns + ((n,k),(i,base)) = rankCombination comb + bins = k + 1 + balls = n - bins + 1 + +-- | Reconstruct a distribution given parameters \((n,k)\) and a rank. +unrankDistribution :: (Int, Int) -> Integer -> [Int] +unrankDistribution (balls,bins) i + | balls < 0 || bins < 0 = err' $ "invalid parameters: " ++ show (balls,bins) + | i < 0 || i >= base = err' $ "out of range: " ++ show (i,base) + | bins == 0 = [] + | otherwise = countGaps 0 bs + where + err' = err . ("unrankDistribution: " ++) + n = balls + bins - 1 -- stars and bars + k = bins - 1 -- number of bars + base = if bins == 0 then 1 else n `choose` k + bs = unrankCombination (n,k) i + + countGaps !acc [] = [acc] + countGaps !acc (False:rest) = countGaps (acc + 1) rest + countGaps !acc (True:rest) = acc : countGaps 0 rest + +-- | Rank a non-empty distribution in the form of a list bin +-- counts. Returns the \((n,k)\) parameters (where \(n\) is the total +-- number of elements and \(k\) is the number of bins), the rank and the +-- total number of distributions with those parameters. +rankDistribution1 :: [Int] -> ((Int, Int), (Integer, Integer)) +rankDistribution1 ns + | any (< 1) ns = if any (< 0) ns then err' "negative count" + else err' "empty count" + | otherwise = ((balls,bins),(i,base)) + where + err' = err . ("rankDistribution1: " ++) + ((balls',bins),(i,base)) = rankDistribution $ (+(-1)) <$> ns + balls = balls' + bins + +-- | Reconstruct a distribution given parameters \((n,k)\) and a rank. +unrankDistribution1 :: (Int, Int) -> Integer -> [Int] +unrankDistribution1 (balls,bins) i + | balls < bins || bins < 0 = + err' $ "invalid parameters: " ++ show (balls,bins) + | otherwise = (+1) <$> unrankDistribution (balls',bins) i + where + err' = err . ("unrankDistribution1: " ++) + balls' = balls - bins diff --git a/src/Codec/Arithmetic/Variety.hs b/src/Codec/Arithmetic/Variety.hs new file mode 100644 --- /dev/null +++ b/src/Codec/Arithmetic/Variety.hs @@ -0,0 +1,130 @@ +{-# LANGUAGE BangPatterns, InstanceSigs #-} +-- | The optimal (shortest) binary code of a value in a domain of +-- uniform probability is simply the binary expansion of the index of +-- the value in that space. The optimal code of two such values is the +-- index of the pair in the cartesian product of both domains, and so on +-- for any number of values. This package defines a type `Value` with a +-- `Monoid` instance that performs this sort of composition. The only +-- difference with typical [arithmetic +-- coding](https://en.wikipedia.org/wiki/Arithmetic_coding) on a +-- rational number code is that for each operation, we operate on the +-- whole code with infinite precision. For an codec with finite +-- precision, see the @Variety.Bounded@ module. +module Codec.Arithmetic.Variety + ( -- * Value-base Interface + + encode + , codeLen + , decode + , encode1 + , codeLen1 + , decode1 + + -- * Value Type + + , Value(..) + , mkValue + , toBitVec + , compose + , maxValue + ) where + +import Codec.Arithmetic.Variety.BitVec (BitVec, bitVec) +import qualified Codec.Arithmetic.Variety.BitVec as BV + +err :: String -> a +err = error . ("Variety." ++) + +-- | Encode a series of value-base pairs into a single bit vector. A +-- base must be at least equal to @1@ and the associated value must +-- exist in the range @[0..base-1]@. +encode :: [(Integer,Integer)] -> BitVec +encode = toBitVec . mconcat . fmap (uncurry mkValue) + +-- | Return the length of the code of a sequence of values in the given +-- list of bases in bits. +codeLen :: [Integer] -> Int +codeLen = codeLen1 . product + +-- | Decode a bit vector given the same series of bases that was used to +-- encode it. Throws an error if the given vector's size doesn't match +-- the given bases. +decode :: [Integer] -> BitVec -> [Integer] +decode bases bv = case init $ scanr (*) 1 bases of -- last is 1 + [] -> [] + (base:ns) -> case compare len expectedLen of -- base == product bases + EQ -> go (BV.toInteger bv) ns + LT -> err "decode: not enough bits" + GT -> err "decode: too many bits" + where + len = BV.length bv + expectedLen = codeLen1 base + where + go i [] = [i] + go i2 (n1:ns) = i0 : go i1 ns + where (i0,i1) = quotRem i2 n1 + +-- | Consider a positive integer as a bit vector, given its base. The +-- base is only required to determine the number of leading 0s. +encode1 :: Integer -> Integer -> BitVec +encode1 = toBitVec .: mkValue + +-- | Return the length of the code of a single value in the given base +-- in bits. +codeLen1 :: Integer -> Int +codeLen1 n | n < 1 = err "codeLen: base must be positive and non-zero" + | otherwise = BV.bitLen $ n - 1 + +-- | Recover the value from a bit vector. +decode1 :: BitVec -> Integer +decode1 = BV.toInteger + +---------------- +-- VALUE TYPE -- +---------------- + +-- | A value with its base, or the number of possible values that could +-- be (i.e. radix, or +-- [variety](https://en.wikipedia.org/wiki/Variety_(cybernetics\))). The +-- value is like an index and ranges from [0..base-1] while the base is +-- a cardinality is always positive and non-zero. +newtype Value = Value { + -- | Recover the value and the base as @Integer@s + fromValue :: (Integer, Integer) +} deriving (Eq,Show,Read) + +-- | Construct from a value and a base. Throws an error if either is +-- negative or if the value is not strictly less than the base. +mkValue :: Integer -> Integer -> Value +mkValue i n | 0 <= i && i < n = Value (i,n) + | otherwise = err $ "mkValue: out of bounds: " ++ show (i,n) + +instance Semigroup Value where + (<>) :: Value -> Value -> Value + (<>) = compose + +instance Monoid Value where + mempty :: Value + mempty = Value (0,1) + +-- | Compose two values into a value of a greater base. This is +-- associative, but not commutative. +compose :: Value -> Value -> Value +compose (Value (i0,n0)) (Value (i1,n1)) = Value (i2, n2) + where + !i2 = i0 * n1 + i1 + !n2 = n0 * n1 + +-- | Maximal possible value as an @Integer@ in the given base. +maxValue :: Value -> Integer +maxValue = (+(-1)) . snd . fromValue + +-- | Drop the base and consider the value as a bit vector. The base +-- conceptually rounds to the next power of 2. +toBitVec :: Value -> BitVec +toBitVec (Value (i,n)) = bitVec (codeLen1 n) i + +(.:) :: (c -> d) -> (a -> b -> c) -> a -> b -> d +(.:) = (.) . (.) +infixr 8 .: +{-# INLINE (.:) #-} diff --git a/src/Codec/Arithmetic/Variety/BitVec.hs b/src/Codec/Arithmetic/Variety/BitVec.hs new file mode 100644 --- /dev/null +++ b/src/Codec/Arithmetic/Variety/BitVec.hs @@ -0,0 +1,211 @@ +{-# LANGUAGE InstanceSigs #-} +module Codec.Arithmetic.Variety.BitVec + ( BitVec + + -- * Construction + , bitVec + + -- * Conversion + , fromBits + , toBits + , fromBytes + , toBytes + , fromInteger + , toInteger + , fromString + , toString + + -- * Methods + , empty + , null + , length + , singleton + , append + , take + , drop + , splitAt + , replicate + , countLeadingZeros + , (!!) + , (!?) + + -- * Extra + , bitLen + ) where + +import Prelude hiding + (null, length, take, drop, splitAt, replicate, (!!), fromInteger, toInteger) +import GHC.Num (integerLog2) +import Control.Exception (assert) + +import Data.Bits ((.&.),(.|.)) +import qualified Data.Bits as Bits +import Data.Word (Word8) +import qualified Data.List as L +import qualified Data.List.Extra as L +import Data.ByteString.Lazy (ByteString) +import qualified Data.ByteString.Lazy as BS + +err :: String -> a +err = error . ("Variety.BitVec: " ++) + +-- | A vector of bits +data BitVec = BitVec !Int !Integer + deriving (Show,Read,Ord,Eq) + +instance Semigroup BitVec where + (<>) :: BitVec -> BitVec -> BitVec + (<>) = append + +instance Monoid BitVec where + mempty :: BitVec + mempty = empty + +-- | Construct a BitVec from a length and Integer. +bitVec :: Int -> Integer -> BitVec +bitVec = BitVec + +----------------- +-- CONVERSIONS -- +----------------- + +-- | Construct from a list of bits. `True` is @1@ and `False` is @0@. +fromBits :: [Bool] -> BitVec +fromBits bs = BitVec len $ L.foldl' Bits.setBit 0 ones + where + len = L.length bs + idxs = [len-1,len-2..0] + ones = fmap snd $ L.filter fst $ zip bs idxs + +-- | Return as a list of bits. `True` is @1@ and `False` is @0@. +toBits :: BitVec -> [Bool] +toBits (BitVec len int) = Bits.testBit int <$> [len-1,len-2..0] + +-- | Construct from a lazy @ByteString@ +fromBytes :: ByteString -> BitVec +fromBytes = fromBits . concatMap unpack8bits . BS.unpack + +-- | Pack the bits into a lazy @ByteString@. Pads the left with @0@s if +-- the length is not a multiple of 8. +toBytes :: BitVec -> ByteString +toBytes v@(BitVec len _) = BS.pack $ fmap pack8bits $ + L.chunksOf 8 $ pad ++ toBits v + where + padLen = (-len) `mod` 8 + pad = assert ((len + padLen) `mod` 8 == 0) $ + L.replicate padLen False + +-- | Read bits from the binary representation of an @Integer@. This +-- excludes the possibility of any leading zeros. Use `bitVec` for more +-- flexible construction. +fromInteger :: Integer -> BitVec +fromInteger int = BitVec sz int + where sz = bitLen int + +-- | Return the @Integer@ representation of the @BitVec@. +toInteger :: BitVec -> Integer +toInteger (BitVec _ int) = int + +-- | Read the code from a list of @0@ and @1@ chars. +fromString :: String -> BitVec +fromString = fromBits . fmap f + where f '0' = False + f '1' = True + f c = err $ "Non-binary char encountered: " ++ [c] + +-- | Return the bits as a list of @0@ and @1@ chars. +toString :: BitVec -> String +toString = fmap f . toBits + where f b = if b then '1' else '0' + +------------- +-- METHODS -- +------------- + +-- | The empty bit vector. +empty :: BitVec +empty = BitVec 0 0 + +-- | Returns `True` iff the bit vector is empty. +null :: BitVec -> Bool +null (BitVec 0 0) = True +null _ = False + +-- | Returns the number of bits in the vector. +length :: BitVec -> Int +length (BitVec len _) = len + +-- | Concatenate two bit vectors. +append :: BitVec -> BitVec -> BitVec +append (BitVec len0 int0) (BitVec len1 int1) = + BitVec (len0 + len1) (Bits.shiftL int0 len1 + int1) + +-- | A vector of length 1 with the given bit. +singleton :: Bool -> BitVec +singleton False = BitVec 1 0 +singleton True = BitVec 1 1 + +-- | @`take` n bv@ returns the bit vector consisting of the first @n@ +-- bits of @bv@. +take :: Int -> BitVec -> BitVec +take n bv@(BitVec len int) + | n <= 0 = empty + | n >= len = bv + | otherwise = BitVec n $ int `Bits.shiftR` (len - n) + +-- | @`drop` n bv@ returns @bv@ with the first @n@ bits removed. +drop :: Int -> BitVec -> BitVec +drop n bv@(BitVec len int) + | n <= 0 = bv + | n >= len = empty + | otherwise = BitVec (len - n) $ + int .&. ((1 `Bits.shiftL` (len - n)) - 1) + +-- | @`splitAt` n bv@ is equivalent to @(`take` n bv, `drop` n bv)@ +splitAt :: Int -> BitVec -> (BitVec, BitVec) +splitAt n bv = (take n bv, drop n bv) + +-- | @`replicate` n b@ constructs a bit vector of length @n@ with @b@ +-- the value of every bit. +replicate :: Int -> Bool -> BitVec +replicate n False = BitVec n 0 +replicate n True = BitVec n (Bits.bit n - 1) + +-- | Count the number of @0@ bits preceeding the first @1@ bit. +countLeadingZeros :: BitVec -> Int +countLeadingZeros (BitVec len int) = len - intLen + where intLen = bitLen int + +-- | Returns the value of a bit at a given index, with @0@ being the +-- index of the most significant bit. +(!!) :: BitVec -> Int -> Bool +(BitVec len int) !! i = Bits.testBit int (len - i - 1) +infixl 9 !! + +-- | Returns the value of a bit at a given index if within bounds, with +-- @0@ being the index of the most significant bit. +(!?) :: BitVec -> Int -> Maybe Bool +(BitVec len int) !? i + | i < 0 || i >= len = Nothing + | otherwise = Just $ Bits.testBit int (len - i - 1) +infixl 9 !? + +------------- +-- HELPERS -- +------------- + +-- | Pack exactly 8 bits into a byte +pack8bits :: [Bool] -> Word8 +pack8bits = L.foldl' f 0 + where f acc b = (acc `Bits.shiftL` 1) .|. fromIntegral (fromEnum b) + +-- | Return the 8 bits that make a byte +unpack8bits :: Word8 -> [Bool] +unpack8bits w = Bits.testBit w <$> [7,6,5,4,3,2,1,0] + +-- | The number of bits in the binary expansion of a positive +-- integer. For consistency with inductive definitions, leading zeros +-- are not considered and so @`bitLen` 0 == 0@. +bitLen :: Integer -> Int +bitLen 0 = 0 +bitLen n = fromIntegral (integerLog2 n) + 1 diff --git a/src/Codec/Arithmetic/Variety/Bounded.hs b/src/Codec/Arithmetic/Variety/Bounded.hs new file mode 100644 --- /dev/null +++ b/src/Codec/Arithmetic/Variety/Bounded.hs @@ -0,0 +1,82 @@ +-- | Since the arithmetic operations of composition might get +-- computationally expensive on very large codes, a similar interface is +-- provided here which produces and consumes bits in chunks whenever +-- spaces are about to grow beyond a certain size given in bytes, at the +-- cost of at most one bit per chunk. +-- +-- While the Haskell language standard defines `Integer` as having no +-- upper bound, GHC most commonly uses the GNU Multiple Precision +-- Arithmetic Library (GMP) as a backend for it, which incurs a limit of +-- 16GiB (or a little over 17GB) on the size of `Integer` values. +module Codec.Arithmetic.Variety.Bounded + ( encode + , codeLen + , decode + ) where + +import Data.Bits (Bits(bit)) + +import qualified Codec.Arithmetic.Variety as V +import Codec.Arithmetic.Variety.BitVec (BitVec) +import qualified Codec.Arithmetic.Variety.BitVec as BV + +err :: String -> a +err = error . ("Variety.Bounded: " ++) + +groupWithinPrec :: (a -> Integer) -> Int -> [a] -> [(Integer,[a])] +groupWithinPrec getBase prec + | prec < 0 = err "negative precision" + | otherwise = ffmap reverse . go 1 [] + where + maxBase = bit (prec*8 + 1) - 1 -- max val with `prec` bytes + go base group [] = filter (not . null . snd) [(base,group)] + go 1 group (a:as) = go (getBase a) (a:group) as + go base group (a:as) + | base' > maxBase = (base,group) : go b [a] as + | otherwise = go base' (a:group) as + where + b = getBase a + base' = base * b +{-# INLINE groupWithinPrec #-} + +-- | Given a max precision in bytes, encode a series of value-base pairs +-- into a single bit vector. Bases must be at least equal to @1@ and the +-- associated values must exist in the range @[0..base-1]@. +encode :: Int -> [(Integer,Integer)] -> BitVec +encode = mconcat + . fmap (V.encode . snd) + .: groupWithinPrec snd + +-- | Return the length of the code of a sequence of values in the given +-- precision and list of bases in bits. +codeLen :: Int -> [Integer] -> Int +codeLen = fromIntegral + . sum + . fmap (V.codeLen1 . fst) + .: groupWithinPrec id + +-- | Try to decode a sequence of values at the head of a bit vector +-- given the same precision and list of bases that was used to encode +-- it. If successful, returns the decoded values and the remainder of +-- the `BitVec`, with the sequence's code removed. Throws an error if +-- the given vector's size doesn't match the given bases. +decode :: Int -> [Integer] -> BitVec -> [Integer] +decode = go .: groupWithinPrec id + where + go [] bv | not (BV.null bv) = err "decode: too many bits" + | otherwise = [] + go ((base,bases):rest) bv + | BV.length hd /= len = err "decode: not enough bits" + | otherwise = V.decode bases hd ++ go rest tl + where + len = BV.bitLen (base - 1) + (hd,tl) = BV.splitAt len bv + +(.:) :: (b -> c) -> (a1 -> a2 -> b) -> a1 -> a2 -> c +(.:) = (.).(.) +infixr 8 .: +{-# INLINE (.:) #-} + +ffmap :: (Functor f, Functor g) => (a -> b) -> f (g a) -> f (g b) +ffmap = fmap . fmap +{-# INLINE ffmap #-} diff --git a/src/Codec/Elias.hs b/src/Codec/Elias.hs new file mode 100644 --- /dev/null +++ b/src/Codec/Elias.hs @@ -0,0 +1,119 @@ +-- | Elias codes are prefix codes for positive, non-zero integers with +-- no prior assumption to their size. +module Codec.Elias + ( -- * Gamma coding + + -- | An Elias gamma code consists of the binary expansion of an + -- integer, preceded by the unary encoding of the length of that + -- expansion in zeros. + + encodeGamma + , decodeGamma + + -- * Delta coding + + -- | An Elias delta code is like an Elias gamma code except that the + -- length is itself coded like a gamma code instead of simply a + -- unary encoding. + + , encodeDelta + , decodeDelta + + -- * Omega coding + + -- | An Elias omega code is the result of recursively encoding the + -- length of binary expansions until a length of @1@ is + -- reached. Since binary expansions are written without any leading + -- zeros, a single @0@ bit marks the end of the code. + + , encodeOmega + , decodeOmega + ) where + +import qualified Data.Bits as Bits +import Data.Bifunctor (Bifunctor(first)) +import Codec.Arithmetic.Variety.BitVec (BitVec) +import qualified Codec.Arithmetic.Variety.BitVec as BV + +boundsError :: a +boundsError = error "Elias: Number must be positive and non-zero" + +-- | Encode a number in a Elias gamma code. Throws an error if the input +-- is not positive and non-zero. +encodeGamma :: Integer -> BitVec +encodeGamma x | x > 0 = BV.replicate n False <> xBits + | otherwise = boundsError + where + xBits = BV.fromInteger x + n = BV.length xBits - 1 + +-- | Try to decode an Elias gamma code at the head of the given bit +-- vector. If successful, returns the decoded value and the remainder of +-- the `BitVec`, with the value code removed. Returns @Nothing@ if the +-- bit vector doesn't contain enough bits to define a number. +decodeGamma :: BitVec -> Maybe (Integer, BitVec) +decodeGamma bv | BV.length xBits /= xLen = Nothing + | otherwise = Just (x, bv'') + where + n = BV.countLeadingZeros bv + xLen = n + 1 + bv' = BV.bitVec (BV.length bv - n) $ BV.toInteger bv -- truncate + (xBits, bv'') = BV.splitAt xLen bv' + x = BV.toInteger xBits + +-- | Encode a number in a Elias delta code. Throws an error if the input +-- is not positive and non-zero. +encodeDelta :: Integer -> BitVec +encodeDelta x | x > 0 = encodeGamma (fromIntegral xLen) <> tailBits + | otherwise = boundsError + where + xBits = BV.fromInteger x + xLen = BV.length xBits + n = xLen - 1 + tailBits = BV.bitVec n $ Bits.clearBit x n -- without leading bit + +-- | Try to decode an Elias delta code at the head of the given bit +-- vector. If successful, returns the decoded value and the remainder of +-- the `BitVec`, with the value code removed. Returns @Nothing@ if the +-- bit vector doesn't contain enough bits to define a number. +decodeDelta :: BitVec -> Maybe (Integer, BitVec) +decodeDelta bv = do + (xLen, bv') <- first fromIntegral <$> decodeGamma bv + let n = xLen - 1 + (xTail, bv'') = BV.splitAt n bv' + xBits = BV.singleton True <> xTail + if BV.length xBits /= xLen then Nothing + else Just (BV.toInteger xBits, bv'') + +-- | Encode a number in a Elias omega code. Throws an error if the input +-- is not positive and non-zero. +encodeOmega :: Integer -> BitVec +encodeOmega x0 | x0 > 0 = go x0 eom + | otherwise = boundsError + where + eom = BV.bitVec 1 0 -- "0" + + go 1 = id -- end + go n = go (len-1) . (bv <>) + where + bv = BV.fromInteger n + len = fromIntegral $ BV.length bv + +-- | Try to decode an Elias omega code at the head of the given bit +-- vector. If successful, returns the decoded value and the remainder of +-- the `BitVec`, with the value code removed. Returns @Nothing@ if the +-- bit vector doesn't contain enough bits to define a number. +decodeOmega :: BitVec -> Maybe (Integer, BitVec) +decodeOmega = go 1 + where + go n bv = do + b <- bv BV.!? 0 -- head + case b of + True | BV.length valBits /= len -> Nothing + | otherwise -> go n' bv' + where + len = fromIntegral n + 1 + (valBits, bv') = BV.splitAt len bv + n' = BV.toInteger valBits + + False -> Just (n, BV.drop 1 bv) -- eom diff --git a/src/Codec/Elias/Natural.hs b/src/Codec/Elias/Natural.hs new file mode 100644 --- /dev/null +++ b/src/Codec/Elias/Natural.hs @@ -0,0 +1,81 @@ +-- | Elias codes are prefix codes for positive, non-zero integers with +-- no prior assumption to their size. +-- +-- To allow for encoding zero, functions of this module add @1@ at +-- encoding time and subtract @1@ at decoding time to support any +-- natural number. +module Codec.Elias.Natural + ( -- * Gamma coding + + -- | An Elias gamma code consists of the binary expansion of an + -- integer, preceded by the unary encoding of the length of that + -- expansion in zeros. + + encodeGamma + , decodeGamma + + -- * Delta coding + + -- | An Elias delta code is like an Elias gamma code except that the + -- length is itself coded like a gamma code instead of simply a + -- unary encoding. + + , encodeDelta + , decodeDelta + + -- * Omega coding + + -- | An Elias omega code is the result of recursively encoding the + -- length of binary expansions until a length of @1@ is + -- reached. Since binary expansions are written without any leading + -- zeros, a single @0@ bit marks the end of the code. + + , encodeOmega + , decodeOmega + ) where + +import Data.Bifunctor (Bifunctor(first)) +import Codec.Arithmetic.Variety.BitVec (BitVec) +import qualified Codec.Elias as E + +boundsError :: a +boundsError = error "Elias.Natural: negative number" + +-- | Encode a number in a Elias gamma code. Throws an error if the input +-- is negative. +encodeGamma :: Integer -> BitVec +encodeGamma n | n >= 0 = E.encodeGamma (n+1) + | otherwise = boundsError + +-- | Try to decode an Elias gamma code at the head of the given bit +-- vector. If successful, returns the decoded value and the remainder of +-- the `BitVec`, with the value code removed. Returns @Nothing@ if the +-- bit vector doesn't contain enough bits to define a number. +decodeGamma :: BitVec -> Maybe (Integer, BitVec) +decodeGamma = fmap (first (+(-1))) . E.decodeGamma + +-- | Encode a number in a Elias delta code. Throws an error if the input +-- is negative. +encodeDelta :: Integer -> BitVec +encodeDelta n | n >= 0 = E.encodeDelta (n+1) + | otherwise = boundsError + +-- | Try to decode an Elias delta code at the head of the given bit +-- vector. If successful, returns the decoded value and the remainder of +-- the `BitVec`, with the value code removed. Returns @Nothing@ if the +-- bit vector doesn't contain enough bits to define a number. +decodeDelta :: BitVec -> Maybe (Integer, BitVec) +decodeDelta = fmap (first (+(-1))) . E.decodeDelta + +-- | Encode a number in a Elias omega code. Throws an error if the input +-- is negative. +encodeOmega :: Integer -> BitVec +encodeOmega n | n >= 0 = E.encodeOmega n + | otherwise = boundsError + +-- | Try to decode an Elias omega code at the head of the given bit +-- vector. If successful, returns the decoded value and the remainder of +-- the `BitVec`, with the value code removed. Returns @Nothing@ if the +-- bit vector doesn't contain enough bits to define a number. +decodeOmega :: BitVec -> Maybe (Integer, BitVec) +decodeOmega = fmap (first (+(-1))) . E.decodeOmega diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -0,0 +1,186 @@ +{-# LANGUAGE ScopedTypeVariables #-} +module Main where + +import Test.HUnit +import Test.QuickCheck +import Control.Monad + +import qualified Codec.Elias as E +import qualified Codec.Arithmetic.Variety as V +import qualified Codec.Arithmetic.Variety.BitVec as BV +import qualified Codec.Arithmetic.Variety.Bounded as VB +import qualified Codec.Arithmetic.Combinatorics as Comb + +main :: IO () +main = do + -- HUnit + void $ runTestTT $ TestLabel "Combinatorics:msp" multisetPermutation + void $ runTestTT $ TestLabel "Combinatorics:permutation" permutation + void $ runTestTT $ TestLabel "Combinatorics:combination" combination + void $ runTestTT $ TestLabel "Combinatorics:distribution" distribution + void $ runTestTT $ TestLabel "Combinatorics:distribution1" distribution1 + void $ runTestTT $ TestLabel "Elias:reference" eliasReference + void $ runTestTT $ TestLabel "Elias:comprehensive" eliasComprehensive + + -- QuickTest + quickCheck (forAll genValueBasePairs prop_encode_decode_roundtrip) + quickCheck (forAll genValueBasePair (uncurry prop_encode1_decode1_roundtrip)) + quickCheck (forAll ((,) <$> choose (1, 10) <*> genValueBasePairs) + (uncurry prop_bounded_encode_decode_roundtrip)) + +multisetPermutation :: Test +multisetPermutation = "Multiset Permutation " ~: TestList $ (<$> samples) $ \s -> + let (params,(_,base)) = Comb.rankMultisetPermutation s + f = fst . snd . Comb.rankMultisetPermutation + . Comb.unrankMultisetPermutation params + in TestList $ (<$> [0..base-1]) $ \i -> show (params, i) ~: f i @?= i + where samples = [ [ ] + , [1] + , [1,1,2] + , [1,1,2,2] + , [1,1,2,3,3,3] ] :: [[Int]] + +----------------- +-- HUNIT TESTS -- +----------------- + +permutation :: Test +permutation = "Permutation " ~: TestList $ (<$> samples) $ \s -> + let (_,base) = Comb.rankPermutation s + f = fst . Comb.rankPermutation . Comb.unrankPermutation s + in TestList $ (<$> [0..base-1]) $ \i -> show (s, i) ~: f i @?= i + where samples = [ [ ] + , [0] + , [0,1] + , [0,1,2] + , [0,1,2,3,4] ] :: [[Int]] + +combination :: Test +combination = "Combination " ~: TestList $ (<$> samples) $ \(n,k)-> + let base = n `Comb.choose` k + f = fst . snd . Comb.rankCombination . Comb.unrankCombination (n,k) + in TestList $ (<$> [0..base-1]) $ \i -> show ((n,k), i) ~: f i @?= i + where samples = [ (0, 0) + , (5, 0) + , (5, 3) + , (5, 5) ] + +distribution :: Test +distribution = "Distribution " ~: TestList $ (<$> samples) $ \params -> + let (_,(_,base)) = Comb.rankDistribution $ Comb.unrankDistribution params 0 + f = fst . snd . Comb.rankDistribution . Comb.unrankDistribution params + in TestList $ (<$> [0..base-1]) $ \i -> show (params, i) ~: f i @?= i + where samples = [ (0, 0) + , (5, 0) + , (5, 3) + , (5, 5) + , (3, 5) ] + +distribution1 :: Test +distribution1 = "Non-empty Distribution " ~: TestList $ (<$> samples) $ \params -> + let (_,(_,base)) = Comb.rankDistribution1 $ Comb.unrankDistribution1 params 0 + f = fst . snd . Comb.rankDistribution1 . Comb.unrankDistribution1 params + in TestList $ (<$> [0..base-1]) $ \i -> show (params, i) ~: f i @?= i + where samples = [ (0, 0) + , (5, 0) + , (5, 3) + , (5, 5) ] + +eliasReference :: Test +eliasReference = TestList [ TestList gamma + , TestList delta + , TestList omega ] + where + gamma = [ "g(1)" ~: E.encodeGamma 1 @?= BV.fromString "1" + , "g(2)" ~: E.encodeGamma 2 @?= BV.fromString "010" + , "g(3)" ~: E.encodeGamma 3 @?= BV.fromString "011" + , "g(4)" ~: E.encodeGamma 4 @?= BV.fromString "00100" + , "g(5)" ~: E.encodeGamma 5 @?= BV.fromString "00101" + , "g(6)" ~: E.encodeGamma 6 @?= BV.fromString "00110" + , "g(7)" ~: E.encodeGamma 7 @?= BV.fromString "00111" + , "g(8)" ~: E.encodeGamma 8 @?= BV.fromString "0001000" + , "g(9)" ~: E.encodeGamma 9 @?= BV.fromString "0001001" + , "g(10)" ~: E.encodeGamma 10 @?= BV.fromString "0001010" + , "g(15)" ~: E.encodeGamma 15 @?= BV.fromString "0001111" + , "g(16)" ~: E.encodeGamma 16 @?= BV.fromString "000010000" + , "g(17)" ~: E.encodeGamma 17 @?= BV.fromString "000010001" ] + + delta = [ "d(1)" ~: E.encodeDelta 1 @?= BV.fromString "1" + , "d(2)" ~: E.encodeDelta 2 @?= BV.fromString "0100" + , "d(3)" ~: E.encodeDelta 3 @?= BV.fromString "0101" + , "d(4)" ~: E.encodeDelta 4 @?= BV.fromString "01100" + , "d(5)" ~: E.encodeDelta 5 @?= BV.fromString "01101" + , "d(6)" ~: E.encodeDelta 6 @?= BV.fromString "01110" + , "d(7)" ~: E.encodeDelta 7 @?= BV.fromString "01111" + , "d(8)" ~: E.encodeDelta 8 @?= BV.fromString "00100000" + , "d(9)" ~: E.encodeDelta 9 @?= BV.fromString "00100001" + , "d(10)" ~: E.encodeDelta 10 @?= BV.fromString "00100010" + , "d(15)" ~: E.encodeDelta 15 @?= BV.fromString "00100111" + , "d(16)" ~: E.encodeDelta 16 @?= BV.fromString "001010000" + , "d(17)" ~: E.encodeDelta 17 @?= BV.fromString "001010001" ] + + omega = [ "o(1)" ~: E.encodeOmega 1 @?= BV.fromString "0" + , "o(2)" ~: E.encodeOmega 2 @?= BV.fromString "100" + , "o(3)" ~: E.encodeOmega 3 @?= BV.fromString "110" + , "o(4)" ~: E.encodeOmega 4 @?= BV.fromString "101000" + , "o(5)" ~: E.encodeOmega 5 @?= BV.fromString "101010" + , "o(6)" ~: E.encodeOmega 6 @?= BV.fromString "101100" + , "o(7)" ~: E.encodeOmega 7 @?= BV.fromString "101110" + , "o(8)" ~: E.encodeOmega 8 @?= BV.fromString "1110000" + , "o(9)" ~: E.encodeOmega 9 @?= BV.fromString "1110010" + , "o(10)" ~: E.encodeOmega 10 @?= BV.fromString "1110100" + , "o(15)" ~: E.encodeOmega 15 @?= BV.fromString "1111110" + , "o(16)" ~: E.encodeOmega 16 @?= BV.fromString "10100100000" + , "o(17)" ~: E.encodeOmega 17 @?= BV.fromString "10100100010" ] + +eliasComprehensive :: Test +eliasComprehensive = TestList [ TestList gamma + , TestList delta + , TestList omega ] + where + gamma = (<$> [10..1000]) $ \i -> + "o(" ++ show i ++ ")" ~: E.decodeGamma (E.encodeGamma i) @?= Just (i, BV.empty) + delta = (<$> [10..1000]) $ \i -> + "o(" ++ show i ++ ")" ~: E.decodeDelta (E.encodeDelta i) @?= Just (i, BV.empty) + omega = (<$> [10..1000]) $ \i -> + "o(" ++ show i ++ ")" ~: E.decodeOmega (E.encodeOmega i) @?= Just (i, BV.empty) + +---------------- +-- QUICKCHECK -- +---------------- + +genValueBasePair :: Gen (Integer, Integer) +genValueBasePair = do + base <- choose (1, 100) + value <- choose (0, base - 1) + return (value, base) + +genValueBasePairs :: Gen [(Integer, Integer)] +genValueBasePairs = do + n <- choose (0, 100) + vectorOf n genValueBasePair + +prop_encode_decode_roundtrip :: [(Integer, Integer)] -> Property +prop_encode_decode_roundtrip pairs = + all (\(v, b) -> v >= 0 && v < b && b >= 1) pairs ==> + let encoded = V.encode pairs + bases = map snd pairs + expected = map fst pairs + decoded = V.decode bases encoded + in decoded === expected + +prop_encode1_decode1_roundtrip :: Integer -> Integer -> Property +prop_encode1_decode1_roundtrip value base = + value >= 0 && value < base && base >= 1 ==> + let encoded = V.encode1 value base + decoded = V.decode1 encoded + in decoded === value + +prop_bounded_encode_decode_roundtrip :: Int -> [(Integer, Integer)] -> Property +prop_bounded_encode_decode_roundtrip prec pairs = + prec > 0 && all (\(v, b) -> v >= 0 && v < b && b >= 1) pairs ==> + let encoded = VB.encode prec pairs + bases = map snd pairs + expected = map fst pairs + decoded = VB.decode prec bases encoded + in decoded === expected diff --git a/variety.cabal b/variety.cabal new file mode 100644 --- /dev/null +++ b/variety.cabal @@ -0,0 +1,73 @@ +cabal-version: 3.0 +name: variety +version: 0.1.0.0 +synopsis: integer arithmetic codes + +description: The @Variety@ module provides functions to optimally encode and decode sequences + of value-base pairs assuming uniform probability. + + If codes get too large and slow to process, @Variety.Bounded@ provides similar + interface with a precision parameter at small cost to code length. + + The @Combinatorics@ module provides functions to optimally encode and decode + common combinatorial objects through ranking and unranking. + + The @Elias@ module provides entirely non-parametric encoding and decoding of + positive integers. The usual definition doesn't allow for an encoding of 0, so a + mapping is baked into the functions in @Elias.Natural@ that shifts the number line + by 1. + +license: MIT +license-file: LICENSE +author: nbos +maintainer: nbos@nbos.ca +homepage: https://github.com/nbos/variety +bug-reports: https://github.com/nbos/variety/issues + +-- copyright: +category: Codec +build-type: Simple + +extra-doc-files: CHANGELOG.md + +source-repository head + type: git + location: https://github.com/nbos/variety + +common warnings + ghc-options: -Wall + +library + import: warnings + exposed-modules: Codec.Arithmetic.Variety + Codec.Arithmetic.Variety.BitVec + Codec.Arithmetic.Variety.Bounded + Codec.Arithmetic.Combinatorics + + Codec.Elias + Codec.Elias.Natural + -- other-modules: + -- other-extensions: + build-depends: base ^>=4.17.2.1 + , bytestring >= 0.11.5 && < 0.12 + , containers >= 0.6.7 && < 0.7 + , exact-combinatorics >= 0.2.0 && < 0.3 + , extra >= 1.8 && < 1.9 + + hs-source-dirs: src + default-language: Haskell2010 + +test-suite variety-test + import: warnings + default-language: Haskell2010 + -- other-modules: + -- other-extensions: + type: exitcode-stdio-1.0 + hs-source-dirs: test + main-is: Main.hs + build-depends: + base ^>=4.17.2.1, + variety, + QuickCheck, + HUnit +