diff --git a/.travis.yml b/.travis.yml new file mode 100644 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,10 @@ +language: haskell + +notifications: + email: false + +install: + cabal install --only-dependencies --enable-tests --enable-benchmarks --force-reinstalls + +script: + cabal configure --enable-tests --enable-benchmark -f testing && cabal build && cabal test diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2013, Sam T. + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Sam T. nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,121 @@ +### Synopsis + +This package provides efficient integer interval sets. + +### Description + +> Persistent... is it trees? + +Yes, Radix trees. Trees are balanced by prefix bits, so we have fast +merge operations, such as union, intersection and difference. Chris +Okasaki and Andrew Gill [shows][int-maps] that Patricia tree based +integer maps might be order of magnitude faster than Red-Black tree +counterparts on this operations. The same apply to integer sets, we +just have keys, but don't have values. + +> That does mean the "dense"? + +That means we keep suffixes in bitmaps and we might pack, say 10, +integers which lies close together in one bitmap. This optimization +doesn't affect execution times for sparse sets, but makes dense sets +much more memory efficient — near 10-50 times less space usage +depending on machine word size and the actual density of the +set. Basically, this let us be 3-4 times less memory efficient +comparing with arrays of tightly packed bits, but see... + +> How suffix compaction is performed? + +There are exist a pretty simple algorithm used in memory allocators +called ["buddy memory allocator"][buddy-alloc]. In a nutshell, we have +a big block which is splitted by half when we remove from one of the +half, and merge then back when we insert. It's somewhat inverse to the +ordinary tree approach — in buddy tree we hold more information about +elements that it _doesn't_ contain, while in prefix tree we hold more +information about elements that it _does_ contain. It's easy to guess +that we should do with it — take the two structures then fuse them +into one to produce a new structure which perform _better_. + +Indeed, the key idea in the design is right here — we switch forth and +back between representations per subtree basis. We intersperse +different representations in different tree branches. It's like +chameleon: + +* If the some subset is _sparse_, we just keep a radix tree with + bitmaps at leafs. + +* If the some subset becomes _full_ we turn it into block. If some + buddy block appears, we join the buddy blocks into one. And so forth. + +That is, we just get a structure that dynamically choose the optimal +representation depending on _density_ of set. Moreover in best case +this lead to huge space savings: + + +``` haskell +> ppStats (fromList [0..123456]) +``` + +gives: + +``` +Bin count: 6 +Tip count: 1 +Fin count: 6 +Size in bytes: 408 +Saved space over dense set: 123072 +Saved space over bytestring: 11879 +Percent saved over dense set: 99.6695821185617% +Percent saved over bytestring: 96.67941727028567% +``` + +The `ppStats` is not an exposed function but you can play with it +using `cabal-dev ghci`. + +I don't know if it is an old idea, but this works just fine. + +> So when this data structure is good choice? + +In many situation. It might be used as persistent and compact +replacement for bool arrays or Data.IntSet with the following advantages: + +* Purity is extremely useful in multithreaded settings — + we could keep a set in a mutable transactional variable or an IORef + and atomically update/modify the set. So it could be used as + replacement for [TArray Int Bool][tarray] as well. +* By merging intervals together we achieve compactness. In best case + some of main operations will take O(1)time and space, so if you need + interval set it's here. +* Fast serizalization: if you are need conversion to/from bytestrings. + Because of bitmaps it's possible to do this conversion _extremely_ fast. + +> How this implementation relate to containers version? + +Heavely based. Essentially we just add the buddy interval compaction, +but it turns out that some operations becomes more complicated and +requires much more effort to implement — in order to maintain the all +tree invariants we need to take into account more cases. This is the +reason why some operations are not implemented yet (e.g. lack of +views), but I hope I'll fix it with the time. + +### Documentation + +For documentation see haddock generated documentation. + +### Build Status + +[![Build Status][status-img]][status-link] + +### Maintainer + +This library is written and maintained by Sam T. + +Feel free to report bugs and suggestions via +[github issue tracker][issues] or the mail. + +[status-img]: https://travis-ci.org/pxqr/intset.png +[status-link]: https://travis-ci.org/pxqr/intset +[issues]: https://github.com/pxqr/intset/issues/new + +[int-maps]: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.37.5452 +[buddy-alloc]: http://en.wikipedia.org/wiki/Buddy_memory_allocation +[tarray]: http://hackage.haskell.org/packages/archive/stm/2.2.0.1/doc/html/Control-Concurrent-STM-TArray.html diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/bench/Main.hs b/bench/Main.hs new file mode 100644 --- /dev/null +++ b/bench/Main.hs @@ -0,0 +1,220 @@ +-- | +-- Copyright : (c) Sam T. 2013 +-- License : BSD3 +-- Maintainer : pxqr.sta@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- Each function should be benchmarked at least in the following modes: +-- +-- * sparse — to see worst case performance. Taking into account +-- both implementations I think [0,64..N] is pretty sparse. +-- +-- * dense — to see expected performance. Again [0,2..N] is pretty +-- dense but not interval yet. +-- +-- * interval — to see best case performance. Set should be one +-- single interval like [0..N]. +-- +-- This should help us unify benchmarks and make it more infomative. +-- +{-# LANGUAGE BangPatterns #-} +module Main (main) where + +import Criterion.Main + +import Data.Bits +import Data.ByteString (ByteString) +import qualified Data.ByteString as B +import Data.IntSet as S +import Data.IntervalSet as SB +import Data.IntervalSet.ByteString as SB +import Data.List as L +import Data.Word + + + + +fromByteString :: ByteString -> S.IntSet +fromByteString = S.fromDistinctAscList . indices 0 . B.unpack + where + indices _ [] = [] + indices n (w : ws) = wordIxs n w ++ indices (n + 8) ws + + wordIxs n w = L.map ((n +) . fst) $ L.filter snd $ zip [0..] (bits w) + + bits i = L.map (testBit i) [0..bitSize (0 :: Word8) - 1] + +main :: IO () +main = defaultMain $ + [ bench "fromList/O-2.5K" $ nf S.fromList [0..2500] + , bench "fromList/O-5K" $ nf S.fromList [0..5000] + , bench "fromList/O-10K" $ nf S.fromList [0..10000] + , bench "fromList/O-20K" $ nf S.fromList [0..20000] + , bench "fromList/O-20K" $ nf S.fromList (L.map (* 10) [0..20000]) + , bench "fromList/S-2.5K" $ nf SB.fromList [0..2500] + , bench "fromList/S-5K" $ nf SB.fromList [0..5000] + , bench "fromList/S-10K" $ nf SB.fromList [0..10000] + , bench "fromList/S-20K" $ nf SB.fromList [0..20000] + , bench "fromList/S-20K-sparse" $ nf SB.fromList (L.map (* 10) [0..20000]) + + , bench "interval/O-200K" $ whnf S.fromDistinctAscList [0..200000] + , bench "interval/S-2M" $ whnf (SB.interval 0) 2000000 + + , let !s = S.fromList [1..50000] in + bench "toList/50K" $ nf S.toList s + + , let !s = SB.fromList [1..50000] in + bench "toList/50K" $ nf SB.toList s + + , let !bs = B.replicate 10000 255 in + bench "bitmap/from-10K-O" $ whnf Main.fromByteString bs + + , let !bs = B.replicate (2 ^ (17 :: Int)) 0 in + bench "bitmap/from-1M-S-empty" $ whnf SB.fromByteString bs + + , let !bs = B.replicate (2 ^ (17 :: Int)) 1 in + bench "bitmap/from-1M-S-dense" $ whnf SB.fromByteString bs + + , let !bs = B.replicate (2 ^ (17 :: Int)) 255 in + bench "bitmap/from-1M-S-buddy" $ whnf SB.fromByteString bs + + , let !s = SB.fromList [0, 2 .. 1000000 * 2] in + bench "bitmap/to-strict-1M-S-dense" $ whnf SB.toByteString s + + , let !s = interval 0 10000000 in + bench "bitmap/to-strict-10M-S-buddy" $ whnf SB.toByteString s + + , let !s = S.fromList [0..1000000] in + bench "member/1M" $ nf (L.all (`S.member` s)) [50000..100000] + + , let !s = SB.fromList [0..1000000] in + bench "member/1M" $ nf (L.all (`SB.member` s)) [50000..100000] + ] + ++ concat + [ splitBenchs + , splitGTBenchs + , splitLTBenchs + + , partBenchs + + , mergeTempl S.union SB.union "union" + , mergeTempl S.intersection SB.intersection "intersection" + , mergeTempl S.difference SB.difference "difference" + , mergeTempl symDiff' SB.symDiff "symmetric-difference" + ] + +symDiff' :: S.IntSet -> S.IntSet -> S.IntSet +symDiff' a b = (a `S.union` b) `S.difference` (a `S.intersection` b) + +partBenchs :: [Benchmark] +partBenchs = complexBench "partition" 10000 (snd . S.partition even) + (snd . SB.partition even) + +splitBenchs :: [Benchmark] +splitBenchs = complexBench "split" 1000000 (chunk S.split) (chunk SB.split) + where + chunk op s = L.foldr ((snd .) . op) s points + points = [100,200..1000000] + +splitGTBenchs :: [Benchmark] +splitGTBenchs = complexBench "split/GT" 1000000 (chunk S.split) (chunkGT SB.splitGT) + where + chunkGT op s = L.foldr op s points + chunk op s = L.foldr ((snd .) . op) s points + points = [100,200..1000000] + +splitLTBenchs :: [Benchmark] +splitLTBenchs = complexBench "split/LT" 1000000 (chunk S.split) (chunkGT SB.splitLT) + where + chunkGT op s = L.foldr op s points + chunk op s = L.foldr ((fst .) . op) s points + points = reverse [100,200..1000000] + +{-------------------------------------------------------------------- + Benchmark Templates +--------------------------------------------------------------------} + +type Name = String + +type Template = Name -- name of benchmark + -> Int -- size of int set + -> (S.IntSet -> S.IntSet) + -> (SB.IntSet -> SB.IntSet) + -> [Benchmark] + +type Gen a = Int -> a + +genericBench :: Name -> Gen S.IntSet -> Gen SB.IntSet -> Template +genericBench _type genA genB name n f g = + [ let !s = genA n in bench (name ++ "/O-" ++ show n ++ "-" ++ _type) $ whnf f s + , let !s = genB n in bench (name ++ "/S-" ++ show n ++ "-" ++ _type) $ whnf g s + ] + +sparseSB :: Gen SB.IntSet +sparseSB n = SB.fromList [0, 64..n * 64 ] + +sparseS :: Gen S.IntSet +sparseS n = S.fromDistinctAscList [0, 64..n * 64 ] + +denseSB :: Gen SB.IntSet +denseSB n = SB.fromList [0, 2 .. n * 2] + +denseS :: Gen S.IntSet +denseS n = S.fromDistinctAscList [0, 2 .. n * 2] + +intervalS :: Gen S.IntSet +intervalS n = S.fromDistinctAscList [0..n] + +intervalSB :: Gen SB.IntSet +intervalSB n = SB.fromList [0..n] + +sparseBench :: Template +sparseBench = genericBench "sparse" sparseS sparseSB + +denseBench :: Template +denseBench = genericBench "dense" denseS denseSB + +intervalBench :: Template +intervalBench = genericBench "interval" intervalS intervalSB + +complexBench :: Template +complexBench name n f g = L.concatMap (\t -> t name n f g) templs + where + templs = [sparseBench, denseBench, intervalBench] + + +mergeTempl :: (S.IntSet -> S.IntSet -> S.IntSet) + -> (SB.IntSet -> SB.IntSet -> SB.IntSet) + -> String -> [Benchmark] +mergeTempl sop bop n = + [ let (!a, !b) = (S.fromList [0,64..10000 * 64], S.fromList [1,65..10000 * 64]) in + bench (n ++"/O-10K-sparse-disjoint") $ whnf (uncurry sop) (a, b) + + , let (!a, !b) = (S.fromList [0,64..10000 * 64], S.fromList [0,64..10000 * 64]) in + bench (n ++"/O-10K-sparse-overlap") $ whnf (uncurry sop) (a, b) + + , let (!a, !b) = (SB.fromList [0,64..10000 * 64], SB.fromList [1,65..10000 * 64]) in + bench (n ++ "/S-10K-sparse-disjoint") $ whnf (uncurry bop) (a, b) + + , let (!a, !b) = (SB.fromList [0,64..10000 * 64], SB.fromList [0,64..10000 * 64]) in + bench (n ++ "/S-10K-sparse-overlap") $ whnf (uncurry bop) (a, b) + + , let (!a, !b) = (S.fromList [0,2..500000 * 2], S.fromList [1,3..500000 * 2]) in + bench (n ++ "/O-500K-dense-disjoint") $ whnf (uncurry sop) (a, b) + + , let (!a, !b) = (S.fromList [0,2..500000 * 2], S.fromList [0,2..500000 * 2]) in + bench (n ++ "/O-500K-dense-overlap") $ whnf (uncurry sop) (a, b) + + , let (!a, !b) = (SB.fromList [0,2..500000 * 2], SB.fromList [1,3..500000 * 2]) in + bench (n ++ "/S-500K-dense-disjoint") $ whnf (uncurry bop) (a, b) + + , let (!a, !b) = (SB.fromList [0,2..500000 * 2], SB.fromList [0,2..500000 * 2]) in + bench (n ++ "/S-500K-dense-overlap") $ whnf (uncurry bop) (a, b) + + , let (!a, !b) = (S.fromList [0..500000], S.fromList [0..500000]) in + bench (n ++ "/O-500K-buddy") $ whnf (uncurry sop) (a, b) + + , let (!a, !b) = (SB.fromList [0..500000], SB.fromList [0..500000]) in + bench (n ++ "/S-500K-buddy") $ whnf (uncurry bop) (a, b) + ] diff --git a/intset.cabal b/intset.cabal new file mode 100644 --- /dev/null +++ b/intset.cabal @@ -0,0 +1,95 @@ +name: intset +version: 0.1.0.0 +category: Data +license: BSD3 +license-file: LICENSE +author: Sam T. +maintainer: Sam T. +copyright: (c) 2013, Sam T. +homepage: https://github.com/pxqr/intset +bug-reports: https://github.com/pxqr/intset/issues +build-type: Simple +cabal-version: >=1.8 +synopsis: Pure, fast and memory efficient integer sets. +description: + + This library provides persistent, time and space efficient integer + sets implemented as dense big-endian patricia trees with buddy + suffixes compaction. In randomized settings this structure expected + to be as fast as Data.IntSet from containers, but if a sets is + likely to have long continuous intervals it should be much faster. + . + [/Release notes/] + . + * /0.1.0.0:/ Initial version. + + +extra-source-files: README.md + .travis.yml + + +source-repository head + type: git + location: git://github.com/pxqr/intset.git + + +flag testing + description: Enable testing stuff and expose internals. + default: False + + +library + exposed-modules: Data.IntervalSet + , Data.IntervalSet.ByteString + other-modules: Data.IntervalSet.Internal + build-depends: base == 4.* + , bits-extras >= 0.1.2 + , bytestring + , deepseq + + hs-source-dirs: src + ghc-options: -Wall + + if flag(testing) + cpp-options: -DTESTING + + +test-suite properties + type: exitcode-stdio-1.0 + main-is: Main.hs + hs-source-dirs: tests + build-depends: base == 4.* + , test-framework + , test-framework-quickcheck2 >= 0.3 + , QuickCheck + , intset + ghc-options: -Wall + if !flag(testing) + buildable: False + + +test-suite fusion + type: exitcode-stdio-1.0 + main-is: Fusion.hs + hs-source-dirs: tests + build-depends: base == 4.*, intset + ghc-options: -Wall -ddump-simpl-stats + -- to avoid -ddump-simpl-stats output + if !flag(testing) + buildable: False + + +benchmark benchmarks + type: exitcode-stdio-1.0 + main-is: Main.hs + hs-source-dirs: bench + build-depends: base == 4.* + , bytestring + , containers >= 0.5.2 + , criterion + , deepseq + , intset + ghc-options: -O3 -Wall -fno-warn-orphans +-- -fllvm + if !flag(testing) + buildable: False diff --git a/src/Data/IntervalSet.hs b/src/Data/IntervalSet.hs new file mode 100644 --- /dev/null +++ b/src/Data/IntervalSet.hs @@ -0,0 +1,113 @@ +-- | +-- Copyright : (c) Sam T. 2013 +-- License : BSD3 +-- Maintainer : pxqr.sta@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- An efficient implementation of dense integer sets based on +-- Big-Endian PATRICIA trees with buddy suffix compression. +-- +-- References: +-- +-- * Fast Mergeable Integer Maps (1998) by Chris Okasaki, Andrew Gill +-- +-- +-- This implementation performs espessially well then set contains +-- long integer invervals like @[0..2047]@ that are just merged into +-- one interval description. This allow to perform many operations +-- in constant time and space. However if set contain sparse +-- integers like @[1,12,7908,234,897]@ the same operations will take +-- /O(min(W, n))/ which is good enough in most cases. +-- +-- Conventions in complexity notation: +-- +-- * n — number of elements in a set; +-- +-- * W — number bits in a 'Key'. This is 32 or 64 at 32 and 64 bit +-- platforms respectively; +-- +-- * O(n) or O(k) — means this operation have complexity O(n) in +-- worst case (e.g. sparse set) or O(k) in best case (e.g. one +-- single interval). +-- +-- Note that some operations will take centuries to compute. For +-- exsample @map id universe@ will never end as well as 'filter' +-- applied to 'universe', 'naturals', 'positives' or 'negatives'. +-- +-- Also note that some operations like 'union', 'intersection' and +-- 'difference' have overriden from default fixity, so use these +-- operations with infix syntax carefully. +-- +{-# LANGUAGE CPP #-} + +#if __GLASGOW_HASKELL__ >= 720 +{-# LANGUAGE Safe #-} +#endif + +module Data.IntervalSet + ( + -- * Types + IntSet(..), Key + + -- * Query + -- ** Cardinality + , SB.null + , size + + -- ** Membership + , member, notMember + + -- ** Inclusion + , isSubsetOf, isSupersetOf +-- , isProperSubsetOf, isProperSupersetOf + + -- * Construction + , empty + , singleton + , interval +{- + , naturals + , negatives + , universe +-} + -- * Modification + , insert + , delete + + -- * Map Fold Filter + , SB.map + , SB.foldr + , SB.filter + + -- * Splits + , split, splitGT, splitLT + , partition + + -- * Min/Max + , findMin, findMax + + -- * Combine + , union, unions + , intersection, intersections + , difference, symDiff + + -- ** Monoids + , Union, Intersection, Difference + + -- * Conversion + -- *** Arbitary + , elems + , toList, fromList + + -- *** Ordered + , toAscList, toDescList + , fromAscList + +#if defined (TESTING) + -- * Debug + , isValid, splitFin +#endif + ) where + +import Data.IntervalSet.Internal as SB diff --git a/src/Data/IntervalSet/ByteString.hs b/src/Data/IntervalSet/ByteString.hs new file mode 100644 --- /dev/null +++ b/src/Data/IntervalSet/ByteString.hs @@ -0,0 +1,185 @@ +-- | +-- Copyright : (c) Sam T. 2013 +-- License : BSD3 +-- Maintainer : pxqr.sta@gmail.com +-- Stability : experimental +-- Portability : little endian arch +-- +-- Fast conversion from or to lazy and strict bytestrings. +-- Serialized IntSets are represented as single continious bitmap. +-- +-- This module is kept separated due safe considerations. +-- +{-# LANGUAGE CPP #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE OverloadedStrings #-} +module Data.IntervalSet.ByteString + ( fromByteString + , toByteString + ) where + +import Data.Bits +import Data.ByteString (ByteString) +import qualified Data.ByteString as BS +import qualified Data.ByteString.Internal as BS +import Foreign + +import Data.IntervalSet.Internal as S + + +#if defined(__GLASGOW_HASKELL__) +#include "MachDeps.h" +#endif + +{- + it seems like we have this conversion hella fast by desing + e.g. read by blocks(bitmaps), fast merge, fast 'bin' + + but we need to make memory access patterns linear and dense + e.g. read left subtree /before/ right subtree; + TODO carefully force this behaviour +-} + +-- | Unpack 'IntSet' from bitmap. +fromByteString :: ByteString -> IntSet +fromByteString bs = + let (fptr, off, len) = BS.toForeignPtr bs in + BS.inlinePerformIO $ withForeignPtr fptr $ \_ptr -> do + let ptr = _ptr `advancePtr` off + let !s = goFrom (castPtr ptr) len + return $! s + where + wordSize = sizeOf (0 :: Word) + + goFrom ptr len = go 0 empty -- goTree 0 len + where + go :: Int -> IntSet -> IntSet + go !x !acc + | x + wordSize <= len = do + let !bm = BS.inlinePerformIO (peekByteOff ptr x) -- TODO read little endian + let !s = unionBM (x * wordSize) bm acc + go (x + wordSize) s + | otherwise = goBytes x acc + + -- normally this loop should run only at the mostleft region of bitmap + -- note that the left index is not necessary multiple of a word size + goBytes :: Int -> IntSet -> IntSet + goBytes !i !s + | i < len = + let wbm = BS.inlinePerformIO (peekByteOff ptr i) + s' = foldrWord (i * 8) insert s wbm + in goBytes (i + 1) s' + | otherwise = s + +{- + goTree :: Int -> Int -> IntSet + goTree !l !r + | traceShow (l, r) False = undefined + | r - l > wordSize = + let !px = l `shiftL` 3 + !qx = r `shiftL` 3 + !msk = branchMask px qx + -- TODO fix mid + !mid = let br = branchMask l r in if br == r + then div (r + l) 2 + else br + in traceShow (l, mid, r, px, qx, msk) $ + bin px msk (goTree l mid) (goTree mid r) + + | r - l == wordSize = + let bm = BS.inlinePerformIO (peekByteOff ptr l) + in tip (l * wordSize) bm + + | otherwise = goBytes l r empty +-} + +foldrWord :: Int -> (Int -> a -> a) -> a -> Word8 -> a +foldrWord p f acc bm = go 0 + where + go i + | i == 8 = acc + | testBit bm i = f (p + i) (go (succ i)) + | otherwise = go (succ i) +{- +-- | Pack 'IntSet' as bitmap to the bytestring builder. +-- +-- NOTE: negative elements are ignored! +-- +toBuilder :: IntSet -> Builder +toBuilder _s = go (splitGT (-1) _s) (\_ -> BS.byteString "") 0 + where + indent n p = BS.byteString $ BS.replicate ((p - n) `shiftR` 3) 0 + {-# INLINE indent #-} + + {-# INLINE wordLE #-} +#if WORD_SIZE_IN_BITS == 64 + wordLE = BS.word64LE +#elif WORD_SIZE_IN_BITS == 32 + wordLE = BS.word32LE +#else +#error Unsupported platform +#endif + -- TODO preallocate buffer and write + -- TODO trim last zeroed bytes + go :: IntSet -> (Int -> Builder) -> (Int -> Builder) + go s c !n = case s of + Bin _ _ l r -> go l (go r c) n + Tip p bm -> indent n p <> + wordLE (fromIntegral bm) <> + c (p + WORD_SIZE_IN_BITS) + Fin p m -> indent n p <> + BS.byteString (BS.replicate (m `shiftR` 3) 255) <> + c (p + m) + Nil -> c n + +-- | Pack the 'IntSet' as bitmap to the lazy bytestring. +-- +-- NOTE: you should prefer 'toLazyByteString' over 'toByteString'. +-- +-- NOTE: negative elements are ignored! +-- +toLazyByteString :: IntSet -> BSL.ByteString +toLazyByteString = BS.toLazyByteString . toBuilder +{-# INLINE toLazyByteString #-} +-} +-- | Pack the 'IntSet' as bitmap to the strict bytestring. +-- +-- NOTE: negative elements are ignored! +-- +toByteString :: IntSet -> ByteString +toByteString snp = + let s = splitGT (-1) snp + maxEl = if S.null snp then 0 else findMax s + 1 + sizeWord = wordSize maxEl + sizeByte = byteSize maxEl + in BS.take sizeByte $ + BS.unsafeCreate (sizeWord * sizeOf (undefined :: BitMap)) + (`start` s) -- createAndTrim + where + wordSize x = (x `div` WORD_SIZE_IN_BITS) + + if (x `mod` WORD_SIZE_IN_BITS) == 0 then 0 else 1 + byteSize x = (x `div` 8) + if (x `mod` 8) == 0 then 0 else 1 + + indent :: Ptr Word8 -> Int -> Int -> IO () + indent ptr n p = void $ BS.memset (ptr `plusPtr` shiftR n 3) 0 + (fromIntegral (shiftR (p - n) 3)) + {-# INLINE indent #-} + + start :: Ptr Word8 -> IntSet -> IO () + start ptr s = void $ write s 0 + where + write :: IntSet -> Int -> IO Int + write s' !n = case s' of + Bin _ _ l r -> write l n >>= write r + Tip p bm -> do + indent ptr n p + pokeByteOff ptr (p `shiftR` 3) bm + return (p + WORD_SIZE_IN_BITS) + + Fin p m -> do + indent ptr n p + _ <- BS.memset (ptr `advancePtr` shiftR p 3) 255 + (fromIntegral (shiftR m 3)) + return (p + m) + + Nil -> return n diff --git a/src/Data/IntervalSet/Internal.hs b/src/Data/IntervalSet/Internal.hs new file mode 100644 --- /dev/null +++ b/src/Data/IntervalSet/Internal.hs @@ -0,0 +1,1543 @@ +-- | +-- Copyright : (c) Sam T. 2013 +-- License : BSD3 +-- Maintainer : pxqr.sta@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- See documentation for module header in Data.IntSet.Buddy. +-- +{-# LANGUAGE CPP #-} + +#if __GLASGOW_HASKELL__ +{-# LANGUAGE DeriveDataTypeable #-} +#endif + +-- TODO use 'seq' instead of bang patterns +{-# LANGUAGE BangPatterns #-} + +#if __GLASGOW_HASKELL__ >= 720 +-- TODO The only unsafe import is Data.Bits.Extras +{-# LANGUAGE Trustworthy #-} +#endif + +module Data.IntervalSet.Internal + ( + -- * Types + IntSet(..), Key + + -- * Monoids + , Union, Intersection, Difference + + -- * Query + -- ** Cardinality + , Data.IntervalSet.Internal.null + , size + + -- ** Membership + , member, notMember + + -- ** Inclusion + , isSubsetOf, isSupersetOf + , isProperSubsetOf, isProperSupersetOf + + -- * Construction + , empty + , singleton + , interval + + , naturals + , negatives + , universe + + -- * Modification + , insert + , delete + + -- * Map & Fold & Filter + , Data.IntervalSet.Internal.map + , Data.IntervalSet.Internal.foldr + , Data.IntervalSet.Internal.filter + + -- * Splits + , split, splitGT, splitLT + , partition + + -- * Min/Max + , findMin, findMax + + -- * Combine + , union, unions + , intersection, intersections + , difference, symDiff + + -- * Conversion + -- ** Lists + -- *** Arbitary + , elems + , toList, fromList + + -- *** Ordered + , toAscList, toDescList + , fromAscList + + -- TODO conversion to bitmap + -- ** Bitmap + + -- do not export this in Data.IntervalSet + -- * Internal + -- ** Types + , Prefix, Mask, BitMap + , finMask, nomatch, match, mask, insertFin, properSubsetOf + , intersectBM + + -- ** Smart constructors + , tip, tipI, tipD, bin + , insertBM + , unionBM + + -- ** Debug +-- , shorter, prefixOf, bitmapOf, branchMask, matchFin, +-- , finSubsetOf + , splitFin + + -- *** Stats + , binCount, tipCount, finCount + , wordCount + , savedSpace + , ppStats + + -- *** Invariants + , isValid + + -- *** Visualization + , showTree, showRaw + , putTree, putRaw + , symDiff' + ) where + + +import Control.DeepSeq +import Data.Bits as Bits +import Data.Bits.Extras +import Data.Data +import qualified Data.List as L +import Data.Monoid +import Data.Ord +import Data.Word +import Text.ParserCombinators.ReadP + + +-- machine specific properties of basic types +#if defined(__GLASGOW_HASKELL__) +#include "MachDeps.h" +#endif + + +-- | Prefix is used to distinguish subtrees by its prefix bits. When +-- new prefix is created its non prefix bits are zeroed. Prefix is +-- big endian. This means that we throw away only least significant +-- bits +type Prefix = Int + +-- | Mask is used to specify mask for branching bit. +-- For exsample if we have mask 0000100 that means: +-- +-- * We do not consider last three bits in prefix. +-- +-- * Branching bit is at position 3 starting from least +-- significant bit. +-- +-- * Prefix mask is 4 bits. (at left of the bitstring) +-- +type Mask = Int +-- TODO Mask = Word? + +-- | Bitmap is used to make intset dense. To achive this we throw away +-- last bits 6 or 7 bits from any(!) prefix and thus any(!) mask +-- should be more than word size in bits. Bitmap by itself contain +-- flags which indicates "is an element present in a set?" by marking +-- suffices indices. For exsample bitmap 01001010 contain elements 2, +-- 5 and 7. +-- +-- One bitmap might contain up to word size in bits (depending on arch) +-- elements. +-- +type BitMap = Word + +-- | Type of IntSet elements. +type Key = Int + + +-- | Integer set. +data IntSet + -- | Layout: prefix up to branching bit, mask for branching bit, + -- left subtree and right subtree. + -- + -- IntSet = Bin: contains elements of left and right subtrees thus + -- just merge to subtrees. All elements of left subtree is less + -- that from right subtree. Except non-negative numbers, they are + -- in left subtree of root bin, if any. + -- + = Bin {-# UNPACK #-} !Prefix {-# UNPACK #-} !Mask !IntSet !IntSet + + -- | Layout: Prefix up to /mask of bitmap size/, and bitmap + -- containing elements starting /from the prefix/. + -- + -- IntSet = Tip: contains elements + -- + | Tip {-# UNPACK #-} !Prefix {-# UNPACK #-} !BitMap + + -- | Layout: Prefix up to /mask of bitmap size/, and mask specifing + -- how large is set. There is no branching bit at all. + -- Tip is never full. + -- + -- IntSet = Fin: contains all elements from prefix to + -- (prefix - mask - 1) + -- + | Fin {-# UNPACK #-} !Prefix {-# UNPACK #-} !Mask + + -- | Empty set. Contains nothing. + | Nil + deriving + ( Eq +#if defined(__GLASGOW_HASKELL__) + , Typeable, Data +#endif + ) + +{-------------------------------------------------------------------- + Invariants +--------------------------------------------------------------------} + +-- | The following invariants should be hold in all subtrees of a set: +-- +-- 1. + Nil is never child of Bin; +-- +-- 2. + Bin is never contain two Fins with masks equal to mask of Bin; +-- +-- 3. - Mask becomes smaller at each child; +-- +-- 4. - Bin, Fin, Tip masks is always power of two; +-- +-- 5. + Fin mask is always greater or equal than size of bits in word; +-- +-- 6. - Bin left subtree contains element each of which is less than +-- each element of right subtree; +-- +-- 7. + Tip bitmap is never full; +-- +-- 8. + Tip mask is multiple of word bit count; +-- +-- See 'binI' to find out when two intsets should be merged into one. +-- +-- TODO check 3, 4, 6 invariants +-- +isValid :: IntSet -> Bool +isValid Nil = True +isValid (Tip p bm) = not (isFull bm) && (p `mod` WORD_SIZE_IN_BITS == 0) +isValid (Fin _ m ) = m >= WORD_SIZE_IN_BITS +isValid (Bin _ _ Nil _ ) = error "Bin _ _ Nil _" +isValid (Bin _ _ _ Nil) = error "Bin _ _ _ Nil" +isValid (Bin _ m (Fin _ m1) (Fin _ m2)) + = not (m == m1 && m == m2) +isValid (Bin _ _ l r) + = isValid l && isValid r + +--isPowerOf2 x = (x - 1 .&. x) == 0 + +{-------------------------------------------------------------------- + Instances +--------------------------------------------------------------------} + +instance Show IntSet where + showsPrec _ s = showString "{" . list (toList s) . showString "}" + where + list [] = showString "" + list [x] = shows x + list (x : xs) = shows x . showString "," . list xs + +instance Read IntSet where + readsPrec _ = readP_to_S $ do + "{" <- readS_to_P lex + xs <- readS_to_P reads `sepBy` (skipSpaces >> char ',') + "}" <- readS_to_P lex + return (fromList xs) + +instance Ord IntSet where + compare = comparing toList + -- TODO make it faster + +instance Monoid IntSet where + mempty = empty + mappend = union + mconcat = unions + +instance Num IntSet where + (+) = union + (*) = intersection + (-) = difference + + negate = Data.IntervalSet.Internal.complement + abs = error "IntervalSet.abs: not implemented" + signum = error "IntervalSet.singum: not implemented" + fromInteger = singleton . fromIntegral + +instance Bounded IntSet where + minBound = empty + maxBound = universe + +instance NFData IntSet where + +{-------------------------------------------------------------------- + Monoids +--------------------------------------------------------------------} + +-- | Monoid under 'union'. Used by default for 'IntSet'. +-- +-- You could use 'Sum' from 'Data.Monoid' as well. +-- +newtype Union = Union { getUnion :: IntSet } + deriving (Show, Read, Eq, Ord) + +instance Monoid Union where + mempty = Union empty + mappend a b = Union (getUnion a `union` getUnion b) + mconcat = Union . unions . L.map getUnion + +-- | Monoid under 'intersection'. +-- +-- You could use 'Product' from 'Data.Monoid' as well. +-- +newtype Intersection = Intersection { getIntersection :: IntSet } + deriving (Show, Read, Eq, Ord) + +instance Monoid Intersection where + mempty = Intersection universe + mappend a b = Intersection (getIntersection a `intersection` getIntersection b) + mconcat = Intersection . intersections . L.map getIntersection + +-- | Monoid under 'symDiff'. +-- +-- Don't mix up 'symDiff' with 'difference'. +-- +newtype Difference = Difference { getDifference :: IntSet } + deriving (Show, Read, Eq, Ord) + +instance Monoid Difference where + mempty = Difference empty + mappend a b = Difference (getDifference a `symDiff` getDifference b) + +{-------------------------------------------------------------------- + Query +--------------------------------------------------------------------} + +-- | /O(1)/. Is this the empty set? +null :: IntSet -> Bool +null Nil = True +null _ = False + +-- | /O(n)/ or /O(1)/. Cardinality of a set. +size :: IntSet -> Int +size (Bin _ _ l r) = size l + size r +size (Tip _ bm ) = popCount bm +size (Fin _ m ) | m > 0 = m + | otherwise = error "IntSet.size: int overflow" +size Nil = 0 + +-- | /O(min(W, n))/ or /O(1)/. +-- Test if the value is element of the set. +-- +member :: Key -> IntSet -> Bool +member !x = go + where + go (Bin p m l r) + | nomatch x p m = False + | zero x m = go l + | otherwise = go r + + go (Tip y bm) = prefixOf x == y && bitmapOf x .&. bm /= 0 + go (Fin p m) = p <= x && (x <= (p + m - 1)) + go Nil = False + +-- | /O(min(W, n))/ or /O(1)/. +-- Test if the value is not an element of the set. +-- +notMember :: Key -> IntSet -> Bool +notMember !x = not . member x +{-# INLINE notMember #-} + + +{-------------------------------------------------------------------- + Query/Inclusion +--------------------------------------------------------------------} + +-- | /O(n + m)/ or /O(1)/. Test if the second set contain each element +-- of the first. +isSubsetOf :: IntSet -> IntSet -> Bool +isSubsetOf t1@(Bin p1 m1 l1 r1) (Bin p2 m2 l2 r2) + | m1 `shorter` m2 = False + | m2 `shorter` m1 = match p1 p2 m2 && matchDown + | otherwise = p1 == p2 && isSubsetOf l1 l2 && isSubsetOf r1 r2 + where + matchDown + | zero p1 m2 = isSubsetOf t1 l2 + | otherwise = isSubsetOf t1 r2 + +isSubsetOf Bin {} Tip {} = False +isSubsetOf (Bin p1 m1 _ _) (Fin p2 m2) + = finMask m2 `shorterEq` m1 && match p1 p2 (finMask m2) + +isSubsetOf Bin {} Nil = False +isSubsetOf t1@(Tip p1 _ ) (Bin p2 m2 l r) + | nomatch p1 p2 m2 = False + | zero p1 m2 = isSubsetOf t1 l + | otherwise = isSubsetOf t1 r + +isSubsetOf (Tip p1 bm1) (Tip p2 bm2) + = p1 == p2 && isSubsetOfBM bm1 bm2 + +isSubsetOf (Tip p1 _ ) (Fin p2 m2 ) = match p1 p2 (finMask m2) +isSubsetOf Tip {} Nil = False +isSubsetOf t1@(Fin p1 m1 ) (Bin p2 m2 l r) + | finMask m1 `shorterEq` m2 = False + | nomatch p1 p2 m2 = False + | zero p1 m2 = isSubsetOf t1 l + | otherwise = isSubsetOf t1 r + +isSubsetOf Fin {} Tip {} = False +isSubsetOf (Fin p1 m1 ) (Fin p2 m2) + = m2 `shorterEq` m1 && match p1 p2 (finMask m2) + +isSubsetOf Fin {} Nil = False +isSubsetOf Nil _ = True + + +isSubsetOfBM :: BitMap -> BitMap -> Bool +isSubsetOfBM bm1 bm2 = bm1 .|. bm2 == bm2 +{-# INLINE isSubsetOfBM #-} + + +-- | /O(n + m)/ or /O(1)/. Test if the second set is subset of the +-- first. +isSupersetOf :: IntSet -> IntSet -> Bool +isSupersetOf = flip isSubsetOf +{-# INLINE isSupersetOf #-} + +-- | /O(n + m)/ or /O(1)/. Test if the first set is proper subset of +-- the other. +isProperSubsetOf :: IntSet -> IntSet -> Bool +isProperSubsetOf = error "isProper subset of" + +-- | /O(n + m)/ or /O(1)/. Test if the second set is proper subset of +-- the first. +isProperSupersetOf :: IntSet -> IntSet -> Bool +isProperSupersetOf = flip isProperSubsetOf +{-# INLINE isProperSupersetOf #-} + +{-------------------------------------------------------------------- + Construction +--------------------------------------------------------------------} + +-- | /O(1)/. The empty set. +empty :: IntSet +empty = Nil +{-# INLINE empty #-} + +-- | /O(1)/. A set containing one element. +singleton :: Key -> IntSet +singleton x = Tip (prefixOf x) (bitmapOf x) +{-# INLINE singleton #-} + +-- TODO make it faster +-- | /O(n)/. Set containing elements from the specified range. +-- +-- > interval a b = fromList [a..b] +-- +interval :: Key -> Key -> IntSet +interval l r + | l < 0 && r >= 0 = go l (-1) `union` go 0 r + | otherwise = go l r + where + go a b + | b < a = empty +-- | a == b = singleton a + | WORD_SIZE_IN_BITS `shorter` m = tip (prefixOf a) (intervalBM a b) + | otherwise = bin p m (interval a (mid - 1)) (interval mid b) + where + mid = p .|. m + p = mask a m + m = branchMask a b + +intervalBM :: Int -> Int -> BitMap +intervalBM a b = + let abm = bitmapOf a + bbm = bitmapOf b + in fromIntegral (Bits.complement (abm - 1) .&. ((bbm - 1) .|. bbm)) + + + +-- | /O(1)/. The set containing all natural numbers. +naturals :: IntSet +naturals = Fin 0 (bit (WORD_SIZE_IN_BITS - 1)) +{-# INLINE naturals #-} + +-- | /O(1)/. The set containing all negative numbers. +negatives :: IntSet +negatives = Fin (bit (WORD_SIZE_IN_BITS - 1)) (bit (WORD_SIZE_IN_BITS - 1)) +{-# INLINE negatives #-} + +-- | /O(1)/. The set containing the all integers it might contain. +universe :: IntSet +universe = Fin (bit (WORD_SIZE_IN_BITS - 1)) 0 +{-# INLINE universe #-} + + +{-------------------------------------------------------------------- + Modification +--------------------------------------------------------------------} + +-- | /O(min(W, n)/ or /O(1)/. Add a value to the set. +insert :: Key -> IntSet -> IntSet +insert !x = insertBM (prefixOf x) (bitmapOf x) + +insertBM :: Prefix -> BitMap -> IntSet -> IntSet +insertBM !kx !bm = go + where -- do not use worker; use wrapper in go + go t@(Bin p m l r) + | nomatch kx p m = join kx (Tip kx bm) p t + | zero kx m = binI p m (insertBM kx bm l) r + | otherwise = binI p m l (insertBM kx bm r) + + go t@(Tip kx' bm') + | kx' == kx = tipI kx (bm .|. bm') + | otherwise = join kx (Tip kx bm) kx' t + + go t@(Fin p m ) + | nomatch kx p (finMask m) = join kx (Tip kx bm) p t + | otherwise = t + + go Nil = Tip kx bm + + +-- | /O(min(n, W))/. Delete a value from the set. +delete :: Key -> IntSet -> IntSet +delete !x = deleteBM (prefixOf x) (bitmapOf x) + +deleteBM :: Prefix -> BitMap -> IntSet -> IntSet +deleteBM !kx !bm = go + where + go t@(Bin p m l r) + | nomatch kx p m = t + | zero kx m = binD p m (deleteBM kx bm l) r + | otherwise = binD p m l (deleteBM kx bm r) + + go t@(Tip kx' bm') + | kx == kx' = tipD kx (bm' .&. Bits.complement bm) + | otherwise = t + + go t@(Fin p m) -- TODO delete 0 universe doesn't work + | nomatch kx p (finMask m) = t + | otherwise = deleteBM kx bm (splitFin p m) + + go Nil = Nil + +{- Note that 'splitFin' always gives inconsistent intset. Resulting +tree doesn't hold Fin "buddy" or "Tip is never full" invariants. +Therefore this function should be used only when we sure we delete at +least one element from the tree later. +-} + +-- | Chunk fin to buddies or to single tip in the end. +splitFin :: Prefix -> Mask -> IntSet +splitFin p m + -- WARN here we have inconsistency - bitmap is full + -- but this will be fixed in next go for any kx bm + -- and this faster (I think) + | m == WORD_SIZE_IN_BITS = Tip p (Bits.complement 0) + | otherwise = Bin p m' (Fin p m') (Fin (p + m') m') + where + m' = intFromNat (natFromInt m `shiftR` 1) -- TODO endian independent +{-# INLINE splitFin #-} + +complement :: IntSet -> IntSet +complement Nil = universe +complement _ = error "complement: not implemented" + +{-------------------------------------------------------------------- + Combine +--------------------------------------------------------------------} + +{-------------------------------------------------------------------- + Union +--------------------------------------------------------------------} + +infixl 6 `union` + +-- | /O(n + m)/ or /O(1)/. Find set which contains elements of both +-- right and left sets. +union :: IntSet -> IntSet -> IntSet +union t1@(Bin p1 m1 l1 r1) t2@(Bin p2 m2 l2 r2) + | shorter m1 m2 = leftiest + | shorter m2 m1 = rightiest + | p1 == p2 = binI p1 m1 (l1 `union` l2) (r1 `union` r2) + | otherwise = join p1 t1 p2 t2 + where + leftiest + | nomatch p2 p1 m1 = join p1 t1 p2 t2 + | zero p2 m1 = binI p1 m1 (l1 `union` t2) r1 + | otherwise = binI p1 m1 l1 (r1 `union` t2) + + rightiest + | nomatch p1 p2 m2 = join p1 t1 p2 t2 + | zero p1 m2 = binI p2 m2 (t1 `union` l2) r2 + | otherwise = binI p2 m2 l2 (t1 `union` r2) + +union t@ Bin {} (Tip p bm) = insertBM p bm t +union t@ Bin {} (Fin p m ) = insertFin p m t +union t@ Bin {} Nil = t +union (Fin p m ) t = insertFin p m t +union (Tip p bm) t = insertBM p bm t +union Nil t = t + + +insertFin :: Prefix -> Mask -> IntSet -> IntSet +insertFin p1 m1 t2@(Bin p2 m2 l r) + | m2 `shorterEq` m1 && match p1 p2 m2 = + if zero p1 m2 + then binI p2 m2 (insertFin p1 m1 l) r + else binI p2 m2 l (insertFin p1 m1 r) + | match p2 p1 (finMask m1) = Fin p1 m1 + | otherwise = join p1 (Fin p1 m1) p2 t2 + +insertFin p1 m1 (Tip p bm) = insertBM p bm (Fin p1 m1) +insertFin p1 m1 (Fin p2 m2 ) -- TODO simplify + | isBuddy p1 m1 p2 m2 = Fin p1 (m1 * 2) + | isBuddy p2 m2 p1 m1 = Fin p2 (m1 * 2) + | properSubsetOf p1 m1 p2 m2 = Fin p2 m2 + | properSubsetOf p2 m2 p1 m1 = Fin p1 m1 + | m1 == m2 && p1 == p2 = Fin p1 m1 + | otherwise = join p1 (Fin p1 m1) p2 (Fin p2 m2) + +insertFin p m Nil = Fin p m + +-- | /O(max(n)^2 * spine)/ or /O(spine)/. +-- The union of list of sets. +unions :: [IntSet] -> IntSet +unions = L.foldl' union empty + + +-- test if the two Fins is good to merge +isBuddy :: Prefix -> Mask -> Prefix -> Mask -> Bool +isBuddy !p1 !m1 !p2 !m2 = m1 == m2 && xor p1 p2 == m1 && p1 .&. m1 == 0 +{-# INLINE isBuddy #-} + +-- used to if one Fin is subset of the another Fin +properSubsetOf :: Prefix -> Mask -> Prefix -> Mask -> Bool +properSubsetOf !p1 !m1 !p2 !m2 = (m2 `shorter` m1) && match p1 p2 (finMask m2) +{-# INLINE properSubsetOf #-} + +unionBM :: Prefix -> BitMap -> IntSet -> IntSet +unionBM !p !bm !t = case tip p bm of + Bin {} -> error "unionBM: impossible" + Fin p' m' -> insertFin p' m' t + Tip p' bm' -> insertBM p' bm' t + Nil -> t + +{-------------------------------------------------------------------- + Intersection +--------------------------------------------------------------------} + +infixl 7 `intersection` + +-- | /O(n + m)/ or /O(1)/. Find maximal common subset of the two given +-- sets. +intersection :: IntSet -> IntSet -> IntSet +intersection t1@(Bin p1 m1 l1 r1) t2@(Bin p2 m2 l2 r2) + | m1 `shorter` m2 = leftiest + | m2 `shorter` m1 = rightiest + | p1 == p2 = binD p1 m1 (intersection l1 l2) (intersection r1 r2) + | otherwise = Nil + where + leftiest + | nomatch p2 p1 m1 = Nil + | zero p2 m1 = intersection l1 t2 + | otherwise = intersection r1 t2 + + rightiest + | nomatch p1 p2 m2 = Nil + | zero p1 m2 = intersection t1 l2 + | otherwise = intersection t1 r2 + +intersection t@ Bin {} (Tip p bm) = intersectBM p bm t +intersection t@ Bin {} (Fin p m) = intersectFin p m t +intersection Bin {} Nil = Nil +intersection (Tip p bm) t = intersectBM p bm t +intersection (Fin p m) t = intersectFin p m t +intersection Nil _ = Nil + + +intersectFin :: Prefix -> Mask -> IntSet -> IntSet +intersectFin p1 m1 t@(Bin p2 m2 l r) + | m2 `shorterEq` m1 && match p1 p2 m2 + = if zero p1 m2 + then intersectFin p1 m1 l + else intersectFin p1 m1 r + | match p2 p1 (finMask m1) = t + | otherwise = Nil + +intersectFin p1 m1 (Tip p2 bm2) + | match p2 p1 (finMask m1) = Tip p2 bm2 + | otherwise = Nil + +-- Fins are never just intersects: +-- * one fin is either subset or superset of the other +-- * or they are disjoint +-- due power of two masks and prefixes +-- +intersectFin p1 m1 (Fin p2 m2) + | finSubsetOf p1 m1 p2 m2 = Fin p1 m1 + | finSubsetOf p2 m2 p1 m1 = Fin p2 m2 + | otherwise = Nil +intersectFin _ _ Nil = Nil + +-- not proper subset, just subset of +finSubsetOf :: Prefix -> Mask -> Prefix -> Mask -> Bool +finSubsetOf p1 m1 p2 m2 = (m2 `shorterEq` m1) && match p1 p2 (finMask m2) + + +intersectBM :: Prefix -> BitMap -> IntSet -> IntSet +intersectBM p1 bm1 (Bin p2 m2 l r) + | nomatch p1 p2 m2 = Nil + | zero p1 m2 = intersectBM p1 bm1 l + | otherwise = intersectBM p1 bm1 r + +intersectBM p1 bm1 (Tip p2 bm2 ) + | p1 == p2 = tipD p1 (bm1 .&. bm2) + | otherwise = Nil + +intersectBM p1 bm1 (Fin p2 m2) + | match p1 p2 (finMask m2) = Tip p1 bm1 + | otherwise = Nil + +intersectBM _ _ Nil = Nil + +-- | /O(max(n) * spine)/ or /O(spine)/. +-- Find out common subset of the list of sets. +intersections :: [IntSet] -> IntSet +intersections = L.foldl' intersection empty + +{-------------------------------------------------------------------- + Difference +--------------------------------------------------------------------} + +-- Since difference is not commutative it's simpler to match all patterns +-- See note for 'splitFin': we _should not_ split when this unnecessary. + +infixl 6 `difference` + +-- | /O(n + m)/ or /O(1)/. Find difference of the two sets. +difference :: IntSet -> IntSet -> IntSet +difference t1@(Bin p1 m1 l1 r1) t2@(Bin p2 m2 l2 r2) + | m1 `shorter` m2 = leftiest + | m2 `shorter` m1 = rightiest + | p1 == p2 = binD p1 m1 (difference l1 l2) (difference r1 r2) + | otherwise = t1 + where + leftiest + | nomatch p2 p1 m1 = t1 + | zero p2 m1 = binD p1 m1 (difference l1 t2) r1 + | otherwise = binD p1 m1 l1 (difference r1 t2) + + rightiest + | nomatch p1 p2 m2 = t1 + | zero p1 m2 = difference t1 l2 + | otherwise = difference t1 r2 + +difference t1@ Bin {} (Tip p bm) = deleteBM p bm t1 +difference t1@(Bin p1 m1 _ _) (Fin p2 m2) + | m1 `shorter` finMask m2 + = if match p2 p1 m1 + then difference t1 (splitFin p2 m2) + else t1 + + | finMask m2 `shorter` m1 + = if match p1 p2 (finMask m2) + then Nil + else t1 + + | p1 == p2 = Nil + | otherwise = t1 + +difference t1@ Bin {} Nil = t1 +difference t1@(Tip p _ ) (Bin p2 m2 l r) + | nomatch p p2 m2 = t1 + | zero p m2 = difference t1 l + | otherwise = difference t1 r + +difference t1@ Tip {} (Tip p bm) = deleteBM p bm t1 +difference t1@(Tip p1 _) (Fin p2 m2 ) -- + | nomatch p1 p2 (finMask m2) = t1 -- + | otherwise = Nil -- + +difference t1@(Tip _ _) Nil = t1 +difference t1@(Fin p1 m1) t2@(Bin p2 m2 l r) + | finMask m1 `shorter` m2 + = if match p2 p1 (finMask m1) + then difference (splitFin p1 m1) t2 + else t1 + + | m2 `shorter` finMask m1 = down + | p1 == p2 = difference (splitFin p1 m1) t2 + | otherwise = t1 + where + down + | nomatch p1 p2 m2 = t1 + | zero p1 m2 = difference t1 l + | otherwise = difference t1 r + +difference t1@(Fin _ _) (Tip p bm) = deleteBM p bm t1 +difference t1@(Fin p1 m1) t2@(Fin p2 m2) + | m1 `shorter` m2 + = if match p2 p1 (finMask m1) + then difference (splitFin p1 m1) t2 + else t1 + + | m2 `shorter` m1 = + if match p1 p2 (finMask m2) + then Nil + else t1 + + | p1 == p2 = Nil + | otherwise = t1 + +difference t1@(Fin _ _) Nil = t1 +difference Nil _ = Nil + +-- i can't see some useful use of difference applied to fold + +{-------------------------------------------------------------------- + Symmetric difference +--------------------------------------------------------------------} + + +symDiff' :: IntSet -> IntSet -> IntSet +symDiff' a b = (a `union` b) `difference` (a `intersection` b) +{-# INLINE symDiff #-} + +-- | /O(n + m)/ or /O(1)/. Find symmetric difference of the two sets: +-- resulting set containts elements that either in first or second +-- set, but not in both simultaneous. +-- +symDiff :: IntSet -> IntSet -> IntSet +symDiff t1@(Bin p1 m1 l1 r1) t2@(Bin p2 m2 l2 r2) + | m1 `shorter` m2 = leftiest + | m2 `shorter` m1 = rightiest + | p1 == p2 = bin p1 m1 (symDiff l1 l2) (symDiff r1 r2) + | otherwise = join p1 t1 p2 t2 + where + leftiest + | nomatch p2 p1 m1 = join p1 t1 p2 t2 + | zero p2 m1 = bin p1 m1 (symDiff l1 t2) r1 -- TODO tune (symDiff l1 t2) + | otherwise = bin p1 m1 l1 (symDiff r1 t2) + + rightiest + | nomatch p1 p2 m2 = join p2 t2 p1 t1 + | zero p1 m2 = bin p2 m2 (symDiff l2 t1) r2 -- TODO tune (symDiff l1 t2) + | otherwise = bin p2 m2 l2 (symDiff r2 t1) + +symDiff t1@ Bin {} (Tip p2 bm2 ) = symDiffTip p2 bm2 t1 +symDiff t1@ Bin {} (Fin p2 m2 ) = symDiffFin p2 m2 t1 +symDiff t1@ Bin {} Nil = t1 +symDiff (Tip p1 bm1 ) t2 = symDiffTip p1 bm1 t2 +symDiff (Fin p1 m1 ) t2 = symDiffFin p1 m1 t2 +symDiff Nil t2 = t2 + + +-- INVARIANT p1 and bm1 should form a Tip, not a Nil or Fin! +symDiffTip :: Prefix -> BitMap -> IntSet -> IntSet +-- {-# INLINE symDiffTip #-} +symDiffTip !p1 !bm1 = go + where + go t2@(Bin p2 m2 l r) + | nomatch p1 p2 m2 = join p1 (Tip p1 bm1) p2 t2 + | zero p1 m2 = bin p2 m2 (symDiffTip p1 bm1 l) r -- not use go + | otherwise = bin p2 m2 l (symDiffTip p1 bm1 r) -- not use go + + go t2@(Tip p2 bm2) + | p1 == p2 = tip p1 (bm1 `xor` bm2) + | otherwise = join p1 (Tip p1 bm1) p2 t2 + + go t2@(Fin p2 m2) + | nomatch p1 p2 (finMask m2) = join p1 (Tip p1 bm1) p2 t2 + | otherwise = symDiffTip p1 bm1 (splitFin p2 m2) -- not use go + + go Nil = Tip p1 bm1 + +symDiffFin :: Prefix -> Mask -> IntSet -> IntSet +symDiffFin !p1 !m1 = go + where + go t2@(Bin p2 m2 l r) + | finMask m1 `shorterEq` m2 + = if match p2 p1 (finMask m1) + then symDiff (splitFin p1 m1) t2 + else join p1 (Fin p1 m1) p2 t2 + + | otherwise = goDown -- TODO inline + where + goDown + | nomatch p1 p2 m2 = join p1 (Fin p1 m1) p2 t2 + | zero p1 m2 = bin p2 m2 (go l) r + | otherwise = bin p2 m2 l (go r) + + go (Fin p2 m2 ) -- TODO try use compare m1 m2 + | m1 `shorter` m2 = if match p2 p1 (finMask m1) + then symDiffFin p2 m2 (splitFin p1 m1) + else join p1 (Fin p1 m1) p2 (Fin p2 m2) + + | m2 `shorter` m1 = if match p1 p2 (finMask m2) + then symDiffFin p2 m2 (splitFin p1 m1) + else join p1 (Fin p1 m1) p2 (Fin p2 m2) + + | p1 == p2 = Nil + -- here we have (m1 == m2 && p1 /= p1) and should check if Fin's are buddies + | xor p1 p2 == m1 = if p1 < p2 + then Fin p1 (m1 * 2) + else Fin p2 (m1 * 2) + + | otherwise = join p1 (Fin p1 m1) p2 (Fin p2 m2) + + go (Tip p2 bm2) = symDiffTip p2 bm2 (Fin p1 m1) + go Nil = Fin p1 m1 + +{-------------------------------------------------------------------- + Strict Pair +--------------------------------------------------------------------} + +{- we use strict pair because almost all operations on intsets are + strict in almost all its arguments; so there is no reason to keep + result of split lazy; + + moreover this gives near 2 times performance improvements +-} + +data SPair a b = !a :*: !b + +unStrict :: SPair a b -> (a, b) +unStrict (a :*: b) = (a, b) + +{-------------------------------------------------------------------- + Splits +--------------------------------------------------------------------} + +-- | /O(min(W, n)/. Split the set such that the left projection of the +-- resulting pair contains elements less than the key and right +-- element contains greater than the key. The exact key is excluded +-- from result: +-- +-- > split 5 (fromList [0..10]) == (fromList [0..4], fromList [6..10]) +-- +-- Performance note: if need only lesser or greater keys, use +-- splitLT or splitGT respectively. +-- +split :: Key -> IntSet -> (IntSet, IntSet) +split !k = unStrict . splitBM (prefixOf k) (bitmapOf k) + +splitBM :: Prefix -> BitMap -> IntSet -> SPair IntSet IntSet +splitBM !px !tbm = root + where + root t@(Bin _ m l r) + | m >= 0 = go t + -- in last two clauses we have {l = positive} and {r = negative} + | px >= 0 = let posLT :*: posGT = go l in (r `union` posLT) :*: posGT + | otherwise = let negLT :*: negGT = go r in negLT :*: (negGT `union` l) + root t = go t + + go t@(Bin p m l r) + | nomatch px p m = if p < px then t :*: Nil else Nil :*: t + | zero px m = let ll :*: lr = go l in ll :*: (lr `union` r) + | otherwise = let rl :*: rr = go r in (l `union` rl) :*: rr + + go t@(Tip p bm) + | px < p = Nil :*: t + | p < px = t :*: Nil + | otherwise = tipD px (bm .&. lowBM) :*: tipD px (bm .&. hghBM) + where + lowBM = tbm - 1 + hghBM = Bits.complement (lowBM + tbm) + + go t@(Fin p m ) + | match px p (finMask m) = go (splitFin p m) + | p < px = t :*: Nil + | otherwise = Nil :*: t -- px < p + + go Nil = Nil :*: Nil + +-- | /O(min(W, n)/. Takes subset such that each element is greater +-- than the specified key. The exact key is excluded from result. +splitGT :: Key -> IntSet -> IntSet +splitGT !k = splitBMGT (prefixOf k) (bitmapOf k) + +splitBMGT :: Prefix -> BitMap -> IntSet -> IntSet +splitBMGT !px !tbm = root + where + root t@(Bin _ m l r) + | m >= 0 = go t + | px >= 0 = go l + | otherwise = let !r' = go r in union r' l + root t = go t + + go t@(Bin p m l r) + | nomatch px p m = if p < px then Nil else t + | zero px m = let !l' = go l in union l' r + | otherwise = go r + + go t@(Tip p bm) + | px < p = t + | p < px = Nil + | otherwise = tipD px (bm .&. hghBM) + where + lowBM = tbm - 1 + hghBM = Bits.complement (lowBM + tbm) + + go t@(Fin p m) + | match px p (finMask m) = go (splitFin p m) + | p < px = Nil + | otherwise = t + + go Nil = Nil + + +-- | /O(min(W, n)/. Takes subset such that each element is less +-- than the specified key. The exact key is excluded from result. +splitLT :: Key -> IntSet -> IntSet +splitLT !x = splitBMLT (prefixOf x) (bitmapOf x) + +splitBMLT :: Prefix -> BitMap -> IntSet -> IntSet +splitBMLT !px !tbm = root + where + root t@(Bin _ m l r) + | m >= 0 = go t + | px >= 0 = r `union` go l + | otherwise = go r + root t = go t + + go t@(Bin p m l r) + | nomatch px p m = if p < px then t else Nil + | zero px m = go l + | otherwise = l `union` go r + + go t@(Tip p bm) + | px < p = Nil + | p < px = t + | otherwise = tipD px (bm .&. lowBM) + where + lowBM = tbm - 1 + + go t@(Fin p m) + | match px p (finMask m) = go (splitFin p m) + | p < px = t + | otherwise = Nil + + go Nil = Nil + +{-------------------------------------------------------------------- + Partition +--------------------------------------------------------------------} + +-- | /O(n)/. Split a set using given predicate. +-- +-- > forall f. fst . partition f = filter f +-- > forall f. snd . partition f = filter (not . f) +-- +partition :: (Key -> Bool) -> IntSet -> (IntSet, IntSet) +partition f = unStrict . go + where + -- TODO use where clauses + go (Bin p m l r) = let ll :*: lr = go l + rl :*: rr = go r + -- in both cases we could have Nil and Fin + in bin p m ll rl :*: bin p m lr rr + go (Tip p bm) = let bm' = filterBitMap p f bm + -- in both cases we could have Nil and Fin + in tip p bm' :*: tip p (bm' `xor` bm) + go (Fin p m) = let (l, r) = L.partition f (listFin p m) + in fromList l :*: fromList r + go Nil = Nil :*: Nil + +{-------------------------------------------------------------------- + Min/max +--------------------------------------------------------------------} + +-- | /O(min(W, n))/ or /O(1)/. Find minimal element of the set. +-- If set is empty then min is undefined. +-- +findMin :: IntSet -> Key +findMin (Bin _ rootM l r) + | rootM < 0 = go r + | otherwise = go l + where + go (Bin _ _ lb _) = go lb + go (Tip p bm) = p + findMinBM bm + go (Fin p _) = p + go Nil = error "findMax.go: Bin Nil invariant failed" + +findMin (Tip p bm) = p + findMinBM bm +findMin (Fin p _) = p +findMin Nil = error "findMin: empty set" + +findMinBM :: BitMap -> Int +findMinBM = fromIntegral . trailingZeros +{-# INLINE findMinBM #-} + + +-- | /O(min(W, n))/ or /O(1)/. Find maximal element of the set. +-- Is set is empty then max is undefined. +-- +findMax :: IntSet -> Key +findMax (Bin _ rootM l r) + | rootM < 0 = go l + | otherwise = go r + where + go (Bin _ _ _ ri) = go ri + go (Tip p bm) = p + findMaxBM bm + go (Fin p m) = p + m - 1 + go Nil = error "findMax.go: Bin Nil invariant failed" + +findMax (Tip p bm) = p + findMaxBM bm +findMax (Fin p m ) = p + m - 1 +findMax Nil = error "findMax: empty set" + +findMaxBM :: BitMap -> Int +findMaxBM x = fromIntegral ((WORD_SIZE_IN_BITS - 1) - leadingZeros x) +{-# INLINE findMaxBM #-} + +{-------------------------------------------------------------------- + Conversion Fusion +--------------------------------------------------------------------} + +stream :: IntSet -> [Key] +stream = toList +{-# NOINLINE stream #-} + +unstream :: [Key] -> IntSet +unstream = fromList +{-# NOINLINE unstream #-} + +{-# RULES + "IntSet/stream/unstream" [~3] forall x. stream (unstream x) = x; + "IntSet/unstream/stream" [~3] forall x. unstream (stream x) = x; + "IntSet/stream/fromList" [ 3] forall x. stream (fromList x) = x; + "IntSet/unstream/toList" [ 3] forall x. toList (unstream x) = x + #-} + +{-------------------------------------------------------------------- + Map/fold/filter +--------------------------------------------------------------------} + +{-# RULES + "IntSet/map/id" Data.IntervalSet.Internal.map id = id + #-} + +-- | /O(n * min(W, n))/. +-- Apply the function to each element of the set. +-- +-- Do not use this operation with the 'universe', 'naturals' or +-- 'negatives' sets. +-- +map :: (Key -> Key) -> IntSet -> IntSet +map f = unstream . L.map f . stream +{-# INLINE map #-} + +-- | /O(n)/. Fold the element using the given right associative +-- binary operator. +-- +foldr :: (Key -> a -> a) -> a -> IntSet -> a +foldr f a = wrap + where + wrap (Bin _ m l r) + | m > 0 = go (go a r) l + | otherwise = go (go a l) r + wrap t = go a t + + go z (Bin _ _ l r) = go (go z r) l + go z (Tip p bm) = foldrBits p f z bm + go z (Fin p m) = L.foldr f z (listFin p m) + go z Nil = z + +-- | /O(n)/. Filter all elements that satisfy the predicate. +-- +-- Do not use this operation with the 'universe', 'naturals' or +-- 'negatives' sets. +-- +filter :: (Key -> Bool) -> IntSet -> IntSet +filter f = go + where + go (Bin p m l r) = binD p m (go l) (go r) + go (Tip p bm) = fromList $ L.filter f $ toList (Tip p bm) -- FIX use foldrBits + go (Fin p m) = fromList $ L.filter f $ listFin p m -- FIX fromDistinctAscList + go Nil = Nil + +listFin :: Prefix -> Mask -> [Key] +listFin p m = [p..(p + m) - 1] + +{-------------------------------------------------------------------- + List conversions +--------------------------------------------------------------------} + +{-# RULES + "IntSet/toList/fromList" forall x. fromList (toList x) = x; + "IntSet/toList/fromList/comp" fromList . toList = id; + "IntSet/fromList/toList" forall x. toList (fromList x) = x; + "IntSet/fromList/toList/comp" toList . fromList = id + #-} + +-- | /O(n * min(W, n))/ or /O(n)/. +-- Create a set from a list of its elements. +-- +fromList :: [Key] -> IntSet +fromList = L.foldl' (flip insert) empty +{-# NOINLINE [3] fromList #-} + +-- | /O(n)/. Convert the set to a list of its elements. +toList :: IntSet -> [Key] +toList = Data.IntervalSet.Internal.foldr (:) [] +{-# NOINLINE [3] toList #-} + + +-- | 'elems' is alias to 'toList' for compatibility. +elems :: IntSet -> [Key] +elems = toList +{-# INLINE elems #-} + +{-------------------------------------------------------------------- + List/Ordered +--------------------------------------------------------------------} + +-- | /O(n)/. +-- Convert the set to a list of its element in ascending order. +toAscList :: IntSet -> [Key] +toAscList = toList +{-# INLINE toAscList #-} + +-- TODO make it faster +-- | /O(n)/. +-- Convert the set to a list of its element in descending order. +toDescList :: IntSet -> [Key] +toDescList = reverse . toAscList +{-# INLINE toDescList #-} + +-- TODO make it faster +-- | Build a set from an ascending list of elements. +fromAscList :: [Key] -> IntSet +fromAscList = fromList +{-# INLINE fromAscList #-} + +{-------------------------------------------------------------------- + Smart constructors +--------------------------------------------------------------------} + +{- + we postfix smart constructors: + + * with 'I' - if we have inserted to the tree given as argument; + * with 'D' - if we have deleted from the tree given as argument; + * without - to denote that we either insert or delete from tree; + + Use more specific version of the constructor when possible to avoid + unneccesary branching (some branches that never executes due to + invariants) and less condition tests. + +-} + +-- used when we insert to the tip +tipI :: Prefix -> BitMap -> IntSet +tipI p bm + | isFull bm = Fin p WORD_SIZE_IN_BITS + | otherwise = Tip p bm +{-# INLINE tipI #-} + +-- used when we delete from the tip +tipD :: Prefix -> BitMap -> IntSet +tipD _ 0 = Nil +tipD p bm = Tip p bm +{-# INLINE tipD #-} + +-- used when we construct from unknown mask +tip :: Prefix -> BitMap -> IntSet +tip p bm + | bm == 0 = Nil + | isFull bm = Fin p WORD_SIZE_IN_BITS + | otherwise = Tip p bm +{-# INLINE tip #-} + +-- used when we insert in left or right subtree +binI :: Prefix -> Mask -> IntSet -> IntSet -> IntSet +-- DONE convert full Tip to Fin, then we can avoid this pattern matching +--binI _ _ (Tip p1 bm1) (Tip p2 bm2) +-- | isFull bm1 && isFull bm2 && xor p1 p2 == WORD_SIZE_IN_BITS +-- = Fin p1 128 + +binI p m (Fin _ m1) (Fin _ m2) + | m1 == m && m2 == m +-- TODO can we simplify this? | m1 == m2 + = Fin p (m * 2) + +binI p m l r = Bin p m l r + +-- used when we delete from left or right subtree +binD :: Prefix -> Mask -> IntSet -> IntSet -> IntSet +binD _ _ Nil r = r +binD _ _ l Nil = l +binD p m l r = Bin p m l r +{-# INLINE binD #-} + +-- used when we don't know kind of transformation +bin :: Prefix -> Mask -> IntSet -> IntSet -> IntSet +bin _ _ Nil r = r +bin _ _ l Nil = l +bin p m (Fin _ m1) (Fin _ m2) + | m1 == m && m2 == m +-- TODO can we simplify this? | m1 == m2 + = Fin p (m * 2) +bin p m l r = Bin p m l r +{-# INLINE bin #-} + + +-- note that join should not merge buddies +join :: Prefix -> IntSet -> Prefix -> IntSet -> IntSet +join p1 t1 p2 t2 + | zero p1 m = Bin p m t1 t2 + | otherwise = Bin p m t2 t1 + where + p = mask p1 m + m = branchMask p1 p2 +{-# INLINE join #-} + + +{-------------------------------------------------------------------- + Debug +--------------------------------------------------------------------} +binCount :: IntSet -> Int +binCount (Bin _ _ l r) = 1 + binCount l + binCount r +binCount _ = 0 + +tipCount :: IntSet -> Int +tipCount (Bin _ _ l r) = tipCount l + tipCount r +tipCount (Tip _ _) = 1 +tipCount _ = 0 + +finCount :: IntSet -> Int +finCount (Bin _ _ l r) = finCount l + finCount r +finCount (Fin _ _) = 1 +finCount _ = 0 + +wordCount :: IntSet -> Int +wordCount (Bin _ _ l r) = 5 + wordCount l + wordCount r +wordCount (Tip _ _) = 3 +wordCount (Fin _ _) = 3 +wordCount Nil = 1 + +origSize :: IntSet -> Int +origSize (Bin _ _ l r) = 5 + origSize l + origSize r +origSize (Tip _ _) = 3 +origSize (Fin _ m) = + let tips = m `div` WORD_SIZE_IN_BITS + bins = tips - 1 + in tips * 3 + bins * 5 +origSize Nil = 1 + +savedSpace :: IntSet -> Int +savedSpace s = origSize s - wordCount s + +bsSize :: IntSet -> Int +bsSize s = findMax s `div` 8 + +ppStats :: IntSet -> IO () +ppStats s = do + putStrLn $ "Bin count: " ++ show (binCount s) + putStrLn $ "Tip count: " ++ show (tipCount s) + putStrLn $ "Fin count: " ++ show (finCount s) + + let treeSize = wordCount s + putStrLn $ "Size in bytes: " ++ show (treeSize * 8) + + let savedSize = savedSpace s + let bssize = bsSize s + let savedSizeBS = bssize - treeSize * 8 + putStrLn $ "Saved space over dense set: " ++ show (savedSize * 8) + putStrLn $ "Saved space over bytestring: " ++ show savedSizeBS + + let orig = origSize s + let per = (fromIntegral savedSize / fromIntegral orig) * (100 :: Double) + let perBS = (fromIntegral savedSizeBS / fromIntegral bssize) * (100 :: Double) + putStrLn $ "Percent saved over dense set: " ++ show per ++ "%" + + putStrLn $ "Percent saved over bytestring: " ++ show perBS ++ "%" + +showTree :: IntSet -> String +showTree = go 0 + where + indent n = replicate (4 * n) ' ' + go n Nil = indent n ++ "{}" + go n (Fin p m) = indent n ++ show p ++ ".." ++ show (p + m - 1) + go n (Tip p bm) = indent n ++ show p ++ " " ++ show bm + go n (Bin p m l r) = concat + [ go (succ n) l, "\n" + , indent n, "+", show p, " ", show m, "\n" + , go (succ n) r + ] + +showRaw :: IntSet -> String +showRaw = go 0 + where + indent n = replicate (4 * n) ' ' + go n Nil = indent n ++ "Nil" + go n (Fin p m) = indent n ++ show p ++ " " ++ show m + go n (Tip p bm) = indent n ++ show p ++ " " ++ show bm + go n (Bin p m l r) = concat + [ go (succ n) l, "\n" + , indent n, "+", show p, " ", show m, "\n" + , go (succ n) r + ] + +putTree :: IntSet -> IO () +putTree = putStrLn . showTree + +putRaw :: IntSet -> IO () +putRaw = putStrLn . showRaw + +{-------------------------------------------------------------------- + Misc +--------------------------------------------------------------------} + +foldrBits :: Int -> (Int -> a -> a) -> a -> BitMap -> a +foldrBits p f acc bm = go 0 + where + go i + | i == WORD_SIZE_IN_BITS = acc + | testBit bm i = f (p + i) (go (succ i)) + | otherwise = go (succ i) + +filterBitMap :: Prefix -> (Key -> Bool) -> BitMap -> BitMap +filterBitMap px f bm = go 0 0 + where + go !i !acc + | i == WORD_SIZE_IN_BITS = acc + | testBit bm i && f (px + i) = go (succ i) (bitmapOfSuffix i .|. acc) + | otherwise = go (succ i) acc +{-# INLINE filterBitMap #-} + +isFull :: BitMap -> Bool +isFull x = x == Bits.complement 0 +{-# INLINE isFull #-} + + + +{-------------------------------------------------------------------- + Some of the later code is taken from Data.IntSet.Base +--------------------------------------------------------------------} +type Nat = Word + +natFromInt :: Int -> Nat +natFromInt = fromIntegral +{-# INLINE natFromInt #-} + +intFromNat :: Nat -> Int +intFromNat = fromIntegral +{-# INLINE intFromNat #-} + +{-------------------------------------------------------------------- + Endian independent bit twiddling +--------------------------------------------------------------------} +zero :: Int -> Mask -> Bool +zero i m = (natFromInt i .&. natFromInt m) == 0 +{-# INLINE zero #-} + +-- Suppose a is largest such that 2^a divides 2*m. +-- Then mask i m is i with the low a bits zeroed out. +mask :: Int -> Mask -> Prefix +mask i m + = maskW (natFromInt i) (natFromInt m) +{-# INLINE mask #-} + +match :: Int -> Prefix -> Mask -> Bool +match i p m = mask i m == p +{-# INLINE match #-} + +{- +matchFin :: Int -> Prefix -> Mask -> Bool +matchFin i p m = match i p m && (m .&. i == m .&. p) +-} + +nomatch :: Int -> Prefix -> Mask -> Bool +nomatch i p m = mask i m /= p +{-# INLINE nomatch #-} + +shorter :: Mask -> Mask -> Bool +shorter m1 m2 = natFromInt m1 > natFromInt m2 +{-# INLINE shorter #-} + +shorterEq :: Mask -> Mask -> Bool +shorterEq m1 m2 = natFromInt m1 >= natFromInt m2 +{-# INLINE shorterEq #-} + +branchMask :: Prefix -> Prefix -> Mask +branchMask p1 p2 + = intFromNat (highestBitMask (natFromInt p1 `xor` natFromInt p2)) +{-# INLINE branchMask #-} + +-- | Return a word where only the highest bit is set. +highestBitMask :: Word -> Word +highestBitMask x1 = + let x2 = x1 .|. x1 `shiftR` 1 + x3 = x2 .|. x2 `shiftR` 2 + x4 = x3 .|. x3 `shiftR` 4 + x5 = x4 .|. x4 `shiftR` 8 + x6 = x5 .|. x5 `shiftR` 16 +#if WORD_SIZE_IN_BITS == 64 + x7 = x6 .|. x6 `shiftR` 32 + in x7 `xor` (x7 `shiftR` 1) +#else + in x6 `xor` (x6 `shiftRL` 1) +#endif +{-# INLINE highestBitMask #-} + +{-------------------------------------------------------------------- + Big endian operations +--------------------------------------------------------------------} +maskW :: Nat -> Nat -> Prefix +maskW i m = intFromNat (i .&. (Bits.complement (m-1) `xor` m)) +{-# INLINE maskW #-} + +finMask :: Mask -> Mask +finMask m = m `shiftR` 1 +{-# INLINE finMask #-} + +{---------------------------------------------------------------------- + Functions that generate Prefix and BitMap of a Key or a Suffix. +----------------------------------------------------------------------} + +suffixBitMask :: Int +suffixBitMask = bitSize (undefined :: Word) - 1 +{-# INLINE suffixBitMask #-} + +prefixBitMask :: Int +prefixBitMask = Bits.complement suffixBitMask +{-# INLINE prefixBitMask #-} + +prefixOf :: Int -> Prefix +prefixOf x = x .&. prefixBitMask +{-# INLINE prefixOf #-} + +suffixOf :: Int -> Int +suffixOf x = x .&. suffixBitMask +{-# INLINE suffixOf #-} + +bitmapOfSuffix :: Int -> BitMap +bitmapOfSuffix s = 1 `shiftL` s +{-# INLINE bitmapOfSuffix #-} + +bitmapOf :: Int -> BitMap +bitmapOf x = bitmapOfSuffix (suffixOf x) +{-# INLINE bitmapOf #-} diff --git a/tests/Fusion.hs b/tests/Fusion.hs new file mode 100644 --- /dev/null +++ b/tests/Fusion.hs @@ -0,0 +1,13 @@ +module Main (main) where + +import Data.IntervalSet as S +import System.Exit + +-- should fuse to id +test :: [Int] -> [Int] +test x = toList (S.map id (fromList x)) + +main :: IO () +main = do + print $ test [] + exitSuccess diff --git a/tests/Main.hs b/tests/Main.hs new file mode 100644 --- /dev/null +++ b/tests/Main.hs @@ -0,0 +1,435 @@ +{-# OPTIONS -fno-warn-orphans #-} +module Main (main) where + +import Control.Applicative hiding (empty) +import Test.QuickCheck hiding ((.&.)) +import Test.Framework +import Test.Framework.Providers.QuickCheck2 + +import Data.List as L (sort, nub, map, filter, minimum, intersect) +import Data.IntervalSet as S +import Data.IntervalSet.ByteString as S +import Data.Monoid + + + +instance Arbitrary IntSet where + arbitrary = oneof [fromList <$> arbitrary, buddy <$> arbitrary] + where + buddy :: [Int] -> IntSet + buddy = fromList . concatMap mk + where + mk i = [i * 64 .. i * 64 + 64] + + shrink (Bin _ _ l r) = [l, r] + shrink (Fin p m) = [splitFin p m] + shrink _ = [] + +prop_empty :: [Int] -> Bool +prop_empty xs = (not . (`member` empty)) `all` xs + +prop_singleton :: Int -> [Int] -> Bool +prop_singleton e = all check + where + check x | x == e = member x (singleton e) + | otherwise = notMember x (singleton e) + +prop_insertLookup :: IntSet -> Int -> Bool +prop_insertLookup s i = member i (insert i s ) + +prop_unionLookup :: [Int] -> [Int] -> Bool +prop_unionLookup a b = all (`member` u) a && all (`member` u) b + where + u = fromList a `union` fromList b + +prop_sorted :: IntSet -> Bool +prop_sorted xs = toList xs == L.nub (L.sort (toList xs)) + +prop_valid :: IntSet -> Bool +prop_valid = isValid + +prop_unionSize :: IntSet -> IntSet -> Bool +prop_unionSize a b = size u >= size a + && size u >= size b + && size u <= size a + size b + where + u = a `union` b + +prop_unionComm :: IntSet -> IntSet -> Bool +prop_unionComm a b = a <> b == b <> a + +prop_unionAssoc :: IntSet -> IntSet -> IntSet -> Bool +prop_unionAssoc a b c = a <> (b <> c) == (a <> b) <> c + +prop_unionLeftId :: IntSet -> Bool +prop_unionLeftId a = mempty <> a == a + +prop_unionRightId :: IntSet -> Bool +prop_unionRightId a = mempty <> a == a + +prop_unionTop :: IntSet -> Bool +prop_unionTop a = a <> a == a + +prop_unionIdemp :: IntSet -> IntSet -> Bool +prop_unionIdemp a b = ((a <> b) <> b) == a <> b + +prop_intersectionSize :: IntSet -> IntSet -> Bool +prop_intersectionSize a b = size i <= size a && size i <= size b + where + i = intersection a b + +prop_intersection :: [Int] -> [Int] -> Bool +prop_intersection a b = + fromList a `intersection` fromList b + == fromList (a `L.intersect` b) + +prop_intersectComm :: IntSet -> IntSet -> Bool +prop_intersectComm a b = (a `intersection` b) == (b `intersection` a) + +prop_intersectAssoc :: IntSet -> IntSet -> IntSet -> Bool +prop_intersectAssoc a b c = + ((a `intersection` b) `intersection` c) + == (a `intersection` (b `intersection` c)) + +prop_intersectLeft :: IntSet -> Bool +prop_intersectLeft a = intersection empty a == empty + +prop_intersectRight :: IntSet -> Bool +prop_intersectRight a = intersection a empty == empty + +prop_intersectBot :: IntSet -> Bool +prop_intersectBot a = (a `intersection` a) == a + +prop_intersectIdemp :: IntSet -> IntSet -> Bool +prop_intersectIdemp a b = ((a `intersection` b) `intersection` b) + == intersection a b + +prop_showRead :: IntSet -> Bool +prop_showRead a = read (show a) == a + +prop_eq :: [Int] -> [Int] -> Bool +prop_eq a b + | a' == b' = fromList a' == fromList b' + | a' /= b' = fromList a' /= fromList b' + | otherwise = error "prop_eq: impossible" + where + a' = nub (sort a) + b' = nub (sort b) + +prop_eqRefl :: IntSet -> Bool +prop_eqRefl a = a == a + +prop_eqSym :: IntSet -> IntSet -> Bool +prop_eqSym a b = (a == b) == (b == a) + +prop_eqTrans :: IntSet -> IntSet -> IntSet -> Bool +prop_eqTrans a b c + | (a == b) && (b == c) = a == c + | otherwise = True + +prop_size :: [Int] -> Bool +prop_size xs = length (nub (sort xs)) == size (fromList xs) + +prop_insertDelete :: Int -> IntSet -> Bool +prop_insertDelete i = notMember i . delete i . insert i + +prop_mapPresSize :: IntSet -> Bool +prop_mapPresSize s = size (S.map (*2) s) == size s + +prop_mapLessSize :: IntSet -> Bool +prop_mapLessSize s = size (S.map (`div` 2) s) <= size s + +prop_mapping :: [Int] -> Bool +prop_mapping xs = toList (S.map (*2) (fromList xs)) == L.map (*2) (nub (sort xs)) + +prop_filterSize :: IntSet -> Bool +prop_filterSize s = size (S.filter even s) <= size s + +prop_filtering :: [Int] -> Bool +prop_filtering xs = S.filter even (fromList xs) == fromList (L.filter even xs) + +prop_min :: [Int] -> Bool +prop_min [] = True +prop_min xs = findMin (fromList xs) == L.minimum xs + +{- +prop_universeMember :: [Int] -> Bool +prop_universeMember = all (`member` universe) + +prop_universeDelete :: Int -> Bool +prop_universeDelete i = notMember i (delete i universe) + +prop_universeInsert :: Int -> Bool +prop_universeInsert i = insert i universe == universe + +prop_universeNatNeg :: Bool +prop_universeNatNeg = naturals <> negatives == universe + +prop_naturals :: [Int] -> Bool +prop_naturals = all check + where + check x | x >= 0 = member x naturals + | otherwise = True + +prop_negatives :: [Int] -> Bool +prop_negatives = all check + where + check x | x < 0 = member x negatives + | otherwise = True +-} + +prop_minInSet :: IntSet -> Bool +prop_minInSet s + | S.null s = True + | otherwise = member (findMin s) s + +prop_minIsTheLess :: IntSet -> Bool +prop_minIsTheLess s + | S.null s = True + | otherwise = all (findMin s <=) (toList s) + +prop_maxInSet :: IntSet -> Bool +prop_maxInSet s + | S.null s = True + | otherwise = member (findMax s) s + +prop_maxIsTheGreatest :: IntSet -> Bool +prop_maxIsTheGreatest s + | S.null s = True + | otherwise = all (<= findMax s) (toList s) + +prop_differenceMember :: IntSet -> IntSet -> Bool +prop_differenceMember a b = all (`notMember` difference a b) (toList b) + +prop_differenceIntersection :: IntSet -> IntSet -> Bool +prop_differenceIntersection a b = (difference a b `intersection` b) == empty + +prop_differenceSize :: IntSet -> IntSet -> Bool +prop_differenceSize a b = size (difference a b) <= size a + + +prop_differenceSubset :: IntSet -> IntSet -> Bool +prop_differenceSubset a b = (a `difference` b) `isSubsetOf` a + + +prop_differenceDeMorgan1 :: IntSet -> IntSet -> IntSet -> Bool +prop_differenceDeMorgan1 a b c = a - b * c == (a - b) + (a - c) + +prop_differenceDeMorgan2 :: IntSet -> IntSet -> IntSet -> Bool +prop_differenceDeMorgan2 a b c = a - (b + c) == (a - b) * (a - c) + +prop_differenceDistributive :: IntSet -> IntSet -> IntSet -> Bool +prop_differenceDistributive a b c = (a + b) - c == (a - c) + (b - c) + +prop_splitPivot :: IntSet -> Key -> Bool +prop_splitPivot s k = all (< k) (toList lt) && all (k <) (toList gt) + where + (lt, gt) = split k s + +prop_splitIntersect :: IntSet -> Key -> Bool +prop_splitIntersect s k = lt * gt == empty + where + (lt, gt) = split k s + +prop_splitGT :: IntSet -> Key -> Bool +prop_splitGT s k = all (k <) (toList (splitGT k s)) + +prop_splitLT :: IntSet -> Key -> Bool +prop_splitLT s k = all (< k) (toList (splitLT k s)) + +prop_interval :: Int -> Int -> Bool +prop_interval a s = interval l r == fromList [l..r] + where + l = a + r = l + min 10000 s + +prop_cmp :: IntSet -> IntSet -> Bool +prop_cmp a b = compare a b == compare (toList a) (toList b) + +prop_numInst :: Int -> Bool +prop_numInst i = fromIntegral i == singleton i + +prop_deleteEmpty :: Int -> Bool +prop_deleteEmpty k = delete k empty == empty + +prop_elems :: IntSet -> Bool +prop_elems s = toList s == elems s + +prop_combine :: [IntSet] -> Bool +prop_combine xs = all check xs + where + check x = us <> x == us && is * x == is + us = mconcat xs + is = intersections xs + +prop_sortIdemp :: [Int] -> Bool +prop_sortIdemp xs = let a = ssort xs in ssort a == a + where + ssort = toList . fromList + +prop_bitmapEncode :: IntSet -> Bool +prop_bitmapEncode xs = fromByteString (toByteString xs') == xs' + where -- we should restrict upper bound otherwise we might have out of memory + xs' = splitGT (-1) $ splitLT 1000000 xs + +prop_partition :: IntSet -> Bool +prop_partition s = fst (S.partition even s) == S.filter even s + && snd (S.partition even s) == S.filter odd s + +prop_subsetSize :: IntSet -> IntSet -> Bool +prop_subsetSize a b + | a `isSubsetOf` b = size a <= size b + | otherwise = True + +prop_supersetSize :: IntSet -> IntSet -> Bool +prop_supersetSize a b + | a `isSupersetOf` b = size a >= size b + | otherwise = True + +prop_subsetSuperset :: IntSet -> IntSet -> Bool +prop_subsetSuperset a b = (a `isSubsetOf` b) || (b `isSubsetOf` a) + || not (S.null (symDiff a b)) + +prop_subsetIntersection :: IntSet -> IntSet -> Bool +prop_subsetIntersection a b = (i `isSubsetOf` a) && (i `isSubsetOf` b) + where + i = a * b + +prop_subsetUnion :: IntSet -> IntSet -> Bool +prop_subsetUnion a b = (a `isSubsetOf` u) && (b `isSubsetOf` u) + where + u = a + b + + +prop_symDiffLeftNeutral :: IntSet -> Bool +prop_symDiffLeftNeutral a = symDiff empty a == a + +prop_symDiffRightNeutral :: IntSet -> Bool +prop_symDiffRightNeutral a = symDiff a empty == a + +prop_symDiffCommutative :: IntSet -> IntSet -> Bool +prop_symDiffCommutative a b = (a `symDiff` b) == b `symDiff` a + +prop_symDiffAssociative :: IntSet -> IntSet -> IntSet -> Bool +prop_symDiffAssociative a b c = ((a `symDiff` b) `symDiff` c) + == (a `symDiff` (b `symDiff` c)) + +prop_symDiffDistr :: IntSet -> IntSet -> IntSet -> Bool +prop_symDiffDistr a b c = (a `intersection` (b `symDiff` c)) + == ((a `intersection` b) `symDiff` (a `intersection` c)) + +prop_symDiffUnion :: IntSet -> IntSet -> Bool +prop_symDiffUnion a b = (a `difference` b) `union` (b `difference` a) + == symDiff a b + +prop_symDiffInter :: IntSet -> IntSet -> Bool +prop_symDiffInter a b = (a `union` b) `difference` (a `intersection` b) + == symDiff a b + +prop_symDiffSizeBound :: IntSet -> IntSet -> Bool +prop_symDiffSizeBound a b = size s <= size a + size b + where + s = symDiff a b + +main :: IO () +main = defaultMain + [ testProperty "empty" prop_empty + , testProperty "singleton" prop_singleton + , testProperty "insertLookup" prop_insertLookup + , testProperty "insert delete" prop_insertDelete + , testProperty "compare" prop_cmp + , testProperty "interval" prop_interval + + , testProperty "subset size" prop_subsetSize + , testProperty "superset size" prop_supersetSize + , testProperty "subset superset exclusion" prop_subsetSuperset + , testProperty "subset of intersection" prop_subsetIntersection + , testProperty "subset of union" prop_subsetUnion + + , testProperty "bitmap_encode" prop_bitmapEncode + + , testProperty "symmetric difference left neutral" prop_symDiffLeftNeutral + , testProperty "symmetric difference right neutral" prop_symDiffRightNeutral + , testProperty "symmetric difference commutative" prop_symDiffCommutative + , testProperty "symmetric difference associative" prop_symDiffAssociative + , testProperty "symmetric difference distributive" prop_symDiffDistr + , testProperty "symmetric difference union" prop_symDiffUnion + , testProperty "symmetric difference intersection" prop_symDiffInter + , testProperty "symmetric difference size upper bound" prop_symDiffSizeBound + + +{- + , testProperty "universe member" prop_universeMember + , testProperty "universe delete" prop_universeDelete + , testProperty "universe insert" prop_universeInsert + , testProperty "universe nat neg" prop_universeNatNeg + , testProperty "naturals" prop_naturals + , testProperty "negatives" prop_negatives +-} + + , testProperty "size" prop_size + , testProperty "sort" prop_sorted + , testProperty "sort idempotent" prop_sortIdemp + + + , testProperty "read . show == id" prop_showRead + + , testProperty "equality" prop_eq + , testProperty "equality reflexivity" prop_eqRefl + , testProperty "equality symmetry" prop_eqSym + , testProperty "equality transitivity" prop_eqTrans + + , testProperty "map preserve size" prop_mapPresSize + , testProperty "map not preserve siz" prop_mapLessSize + , testProperty "mapping" prop_mapping + + , testProperty "filter size" prop_filterSize + , testProperty "filtering" prop_filtering + + , testProperty "union size" prop_unionSize + , testProperty "unionLookup" prop_unionLookup + , testProperty "union commutative" prop_unionComm + , testProperty "union associative" prop_unionAssoc + , testProperty "union left identity" prop_unionLeftId + , testProperty "union right identity" prop_unionRightId + , testProperty "union top" prop_unionTop + , testProperty "union idemp" prop_unionIdemp + + , testProperty "minimal in set" prop_minInSet + , testProperty "maximal in set" prop_maxInSet + , testProperty "minimal is the less" prop_minIsTheLess + , testProperty "maximal is the greatest" prop_maxIsTheGreatest + + + , testProperty "intersection size" prop_intersectionSize + , testProperty "intersection lists" prop_intersection + , testProperty "intersection commutative" prop_intersectComm + , testProperty "intersection associative" prop_intersectAssoc + , testProperty "intersection idemp" prop_intersectIdemp + , testProperty "intersection left empty" prop_intersectLeft + , testProperty "intersection right empty" prop_intersectRight + , testProperty "intersection bot" prop_intersectBot + + , testProperty "difference member" prop_differenceMember + , testProperty "difference intersection" prop_differenceIntersection + , testProperty "difference size" prop_differenceSize + , testProperty "difference de morgan1" prop_differenceDeMorgan1 + , testProperty "difference de morgan2" prop_differenceDeMorgan2 + , testProperty "difference distributive" prop_differenceDistributive + , testProperty "difference subset" prop_differenceSubset + + , testProperty "split pivot" prop_splitPivot + , testProperty "split intersection" prop_splitIntersect + , testProperty "split greater than" prop_splitGT + , testProperty "split lesser than" prop_splitLT + , testProperty "partition filter" prop_partition + + , testProperty "min" prop_min + , testProperty "valid" prop_valid + + -- for coverage mostly + , testProperty "delete from empty" prop_deleteEmpty + , testProperty "combine" prop_combine + , testProperty "num instance" prop_numInst + , testProperty "elems" prop_elems + ]