packages feed

variety (empty) → 0.1.0.0

raw patch · 10 files changed

+1189/−0 lines, 10 filesdep +HUnitdep +QuickCheckdep +base

Dependencies added: HUnit, QuickCheck, base, bytestring, containers, exact-combinatorics, extra, variety

Files

+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# Revision history for variety++## 0.1.0.0 -- 2025-06-04++* First version.
+ LICENSE view
@@ -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.
+ src/Codec/Arithmetic/Combinatorics.hs view
@@ -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
+ src/Codec/Arithmetic/Variety.hs view
@@ -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 (.:) #-}
+ src/Codec/Arithmetic/Variety/BitVec.hs view
@@ -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
+ src/Codec/Arithmetic/Variety/Bounded.hs view
@@ -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 #-}
+ src/Codec/Elias.hs view
@@ -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
+ src/Codec/Elias/Natural.hs view
@@ -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
+ test/Main.hs view
@@ -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
+ variety.cabal view
@@ -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+