diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # bitvec [![Hackage](https://img.shields.io/hackage/v/bitvec.svg)](https://hackage.haskell.org/package/bitvec) [![Stackage LTS](https://www.stackage.org/package/bitvec/badge/lts)](https://www.stackage.org/lts/package/bitvec) [![Stackage Nightly](https://www.stackage.org/package/bitvec/badge/nightly)](https://www.stackage.org/nightly/package/bitvec) -A newtype over `Bool` with a better `Vector` instance: 8x less memory, up to 1000x faster. +A newtype over `Bool` with a better `Vector` instance: 8x less memory, up to 3500x faster. The [`vector`](https://hackage.haskell.org/package/vector) package represents unboxed arrays of `Bool`s @@ -12,16 +12,21 @@ the most significant degradation happens for random writes (up to 10% slower). On the other hand, for certain bulk bit operations -`Vector Bit` is up to 1000x faster than `Vector Bool`. +`Vector Bit` is up to 3500x faster than `Vector Bool`. ## Thread safety -* `Data.Bit` is faster, but writes and flips are thread-unsafe. +* `Data.Bit` is faster, but writes and flips are not thread-safe. This is because naive updates are not atomic: they read the whole word from memory, then modify a bit, then write the whole word back. + Concurrently modifying non-intersecting slices of the same underlying array + may also lead to unexpected results, since they can share a word in memory. * `Data.Bit.ThreadSafe` is slower (usually 10-20%), but writes and flips are thread-safe. + Additionally, concurrently modifying non-intersecting slices of the same underlying array + works as expected. However, operations that affect multiple elements are not + guaranteed to be atomic. ## Quick start @@ -104,7 +109,7 @@ One may notice that the order of the inner traversal by `i` does not matter and get tempted to run it in several parallel threads. In this case it is vital to switch from `Data.Bit` to `Data.Bit.ThreadSafe`, -because the former is thread-unsafe with regards to writes. +because the former is not thread-safe with regards to writes. There is a moderate performance penalty (usually 10-20%) for using the thread-safe interface. @@ -120,12 +125,52 @@ * As a 64k-long unboxed `Vector Bool`, implementing union as `zipWith (||)`. * As a 64k-long unboxed `Vector Bit`, implementing union as `zipBits (.|.)`. -When the `libgmp` flag is enabled, +When the `simd` flag is enabled, according to our benchmarks (see `bench` folder), -the union of `Vector Bit` evaluates 24x-36x faster +the union of `Vector Bit` evaluates magnitudes faster than the union of not-too-sparse `IntSet`s -and stunningly outperforms `Vector Bool` by 500x-1000x. +and stunningly outperforms `Vector Bool`. +Here are benchmarks on MacBook M2: +``` +union + 16384 + Vector Bit: + 61.2 ns ± 3.2 ns + Vector Bool: + 96.1 μs ± 4.5 μs, 1570.84x + IntSet: + 2.15 μs ± 211 ns, 35.06x + 32768 + Vector Bit: + 143 ns ± 7.4 ns + Vector Bool: + 225 μs ± 16 μs, 1578.60x + IntSet: + 4.34 μs ± 429 ns, 30.39x + 65536 + Vector Bit: + 249 ns ± 18 ns + Vector Bool: + 483 μs ± 28 μs, 1936.42x + IntSet: + 8.77 μs ± 835 ns, 35.18x + 131072 + Vector Bit: + 322 ns ± 30 ns + Vector Bool: + 988 μs ± 53 μs, 3071.83x + IntSet: + 17.6 μs ± 1.6 μs, 54.79x + 262144 + Vector Bit: + 563 ns ± 27 ns + Vector Bool: + 2.00 ms ± 112 μs, 3555.36x + IntSet: + 36.8 μs ± 3.3 μs, 65.40x +``` + ## Binary polynomials Binary polynomials are polynomials with coefficients modulo 2. @@ -147,15 +192,9 @@ ## Package flags -* Flag `libgmp`, disabled by default. +* Flag `simd`, enabled by default. - Link against the [GMP](https://gmplib.org/) library for the ultimate performance of - `zipBits`, `invertBits` and `countBits`. GMP is readily available on most machines - ([`apt-get install libgmp-dev`](https://packages.ubuntu.com/focal/libgmp-dev) on Ubuntu, - [`brew install gmp`](https://formulae.brew.sh/formula/gmp) - or [`port install gmp`](https://ports.macports.org/port/gmp/summary) on macOS, - [`pacman -S mingw-w64-x86_64-gmp`](https://packages.msys2.org/package/mingw-w64-x86_64-gmp) on MinGW), - so users are strongly encouraged to enable it whenever possible. + Use a C SIMD implementation for the ultimate performance of `zipBits`, `invertBits` and `countBits`. ## Similar packages diff --git a/bench/Bench.hs b/bench/Bench.hs --- a/bench/Bench.hs +++ b/bench/Bench.hs @@ -4,6 +4,7 @@ import Test.Tasty.Patterns.Printer import Bench.BitIndex +import Bench.Common import Bench.GCD import Bench.Intersection import Bench.Invert @@ -18,27 +19,24 @@ main :: IO () main = defaultMain $ map (mapLeafBenchmarks addCompare) - [ bgroup "bitIndex" $ map benchBitIndex [5..14] - , bgroup "invert" $ map benchInvert [5..14] + [ bgroup "add" $ map benchAdd [5..14] + , bgroup "bitIndex" $ map benchBitIndex [5..14] + , bgroup "flip" $ map benchRandomFlip [5..14] , bgroup "gcdExt" $ map benchGCD [5..14] , bgroup "intersection" $ map benchIntersection [5..14] + , bgroup "invert" $ map benchInvert [5..14] , bgroup "product" $ map benchProduct [5..14] , bgroup "productShort" $ map benchProductShort [5..14] - , bgroup "square" $ map benchSquare [5..14] - , bgroup "write" $ map benchRandomWrite [5..14] - , bgroup "flip" $ map benchRandomFlip [5..14] , bgroup "read" $ map benchRandomRead [5..14] - , bgroup "remainder" $ map benchRemainder [5..14] + , bgroup "remainder" $ map benchRemainder [5..11] , bgroup "reverse" $ map benchReverse [5..14] - , bgroup "add" $ map benchAdd [5..14] + , bgroup "square" $ map benchSquare [5..14] , bgroup "sum" $ map benchSum [5..14] - , bgroup "union" $ map benchUnion [5..14] + , bgroup "union" $ map benchUnion [5..20] + , bgroup "write" $ map benchRandomWrite [5..14] ] -bitBenchName :: String -bitBenchName = "Bit" - addCompare :: ([String] -> Benchmark -> Benchmark) addCompare (name : path) - | name /= bitBenchName = bcompare (printAwkExpr (locateBenchmark (bitBenchName : path))) + | name /= labelBit = bcompare (printAwkExpr (locateBenchmark (labelBit : path))) addCompare _ = id diff --git a/bench/Bench/BitIndex.hs b/bench/Bench/BitIndex.hs --- a/bench/Bench/BitIndex.hs +++ b/bench/Bench/BitIndex.hs @@ -9,23 +9,16 @@ import qualified Data.Vector.Unboxed.Mutable as MU import Test.Tasty.Bench -randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec f k = U.generate n (\i -> f (i == n - 1)) +import Bench.Common + +allFalseButLast :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a +allFalseButLast f k = U.generate n (\i -> f (i == n - 1)) where n = 1 `shiftL` k benchBitIndex :: Int -> Benchmark benchBitIndex k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf bitIndexBit (randomVec Bit k) - , bench "BitTS" $ nf bitIndexBitTS (randomVec TS.Bit k) - , bench "Vector" $ nf elemIndexVector (randomVec id k) + [ bench labelBit $ nf (bitIndex (Bit True)) (allFalseButLast Bit k) + , bench labelBitTS $ nf (TS.bitIndex (TS.Bit True)) (allFalseButLast TS.Bit k) + , bench labelVector $ nf (U.elemIndex True) (allFalseButLast id k) ] - -bitIndexBit :: U.Vector Bit -> Maybe Int -bitIndexBit = bitIndex (Bit True) - -bitIndexBitTS :: U.Vector TS.Bit -> Maybe Int -bitIndexBitTS = TS.bitIndex (TS.Bit True) - -elemIndexVector :: U.Vector Bool -> Maybe Int -elemIndexVector = U.elemIndex True diff --git a/bench/Bench/Common.hs b/bench/Bench/Common.hs new file mode 100644 --- /dev/null +++ b/bench/Bench/Common.hs @@ -0,0 +1,95 @@ +module Bench.Common + ( labelBit + , labelBitTS + , labelVector + , labelIntSet + , labelInteger + , randomBools + , randomBools2 + , randomVec + , randomVec2 + , randomSet + , randomSet2 + , randomInteger + , randomInteger2 + , randomIndices + , randomIndicesAndBools + ) where + +import Data.Bit +import Data.Bits +import qualified Data.IntSet as IS +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Unboxed.Mutable as MU +import System.Random + +labelBit :: String +labelBit = "Vector Bit" + +labelBitTS :: String +labelBitTS = "Vector TS.Bit" + +labelVector :: String +labelVector = "Vector Bool" + +labelIntSet :: String +labelIntSet = "IntSet" + +labelInteger :: String +labelInteger = "Integer" + +seed1 :: Int +seed1 = 42 + +seed2 :: Int +seed2 = 123 + +mkRandomBools :: Int -> [Bool] +mkRandomBools seed = map (> (0 :: Int)) $ randoms $ mkStdGen seed + +randomBools :: [Bool] +randomBools = mkRandomBools seed1 + +randomBools2 :: [Bool] +randomBools2 = mkRandomBools seed2 + +mkRandomVec :: MU.Unbox a => Int -> (Bool -> a) -> Int -> U.Vector a +mkRandomVec seed f k = U.fromList $ map f $ take n $ mkRandomBools seed + where + n = 1 `shiftL` k + +randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a +randomVec = mkRandomVec seed1 + +randomVec2 :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a +randomVec2 = mkRandomVec seed2 + +mkRandomSet :: Int -> Int -> IS.IntSet +mkRandomSet seed k = IS.fromAscList (map fst (filter snd (zip [0..] (take n (mkRandomBools seed))))) + where + n = 1 `shiftL` k + +randomSet :: Int -> IS.IntSet +randomSet = mkRandomSet seed1 + +randomSet2 :: Int -> IS.IntSet +randomSet2 = mkRandomSet seed2 + +mkRandomInteger :: Int -> Int -> Integer +mkRandomInteger seed k = toInteger $ toF2Poly $ mkRandomVec seed Bit k + +randomInteger :: Int -> Integer +randomInteger = mkRandomInteger seed1 + +randomInteger2 :: Int -> Integer +randomInteger2 = mkRandomInteger seed2 + +randomIndices :: [Int] +randomIndices = map fst randomIndicesAndBools + +randomIndicesAndBools :: [(Int, Bool)] +randomIndicesAndBools + = map (\x -> if x > 0 then (x, True) else (x .&. maxBound, False)) + . randoms + . mkStdGen + $ 42 diff --git a/bench/Bench/GCD.hs b/bench/Bench/GCD.hs --- a/bench/Bench/GCD.hs +++ b/bench/Bench/GCD.hs @@ -5,26 +5,12 @@ import Data.Bit import qualified Data.Bit.ThreadSafe as TS import Data.Bits -import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomBools :: [Bool] -randomBools = map (> (0 :: Int)) $ randoms $ mkStdGen 42 - -randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec f k = U.fromList $ map f $ take n randomBools - where - n = 1 `shiftL` k - -randomVec' :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec' f k = U.fromList $ map f $ take n $ drop n randomBools - where - n = 1 `shiftL` k +import Bench.Common benchGCD :: Int -> Benchmark benchGCD k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf (uncurry gcdExt) ( toF2Poly $ randomVec Bit k, toF2Poly $ randomVec' Bit k) - , bench "BitTS" $ nf (uncurry TS.gcdExt) (TS.toF2Poly $ randomVec TS.Bit k, TS.toF2Poly $ randomVec' TS.Bit k) + [ bench labelBit $ nf (uncurry gcdExt) ( toF2Poly $ randomVec Bit k, toF2Poly $ randomVec2 Bit k) + , bench labelBitTS $ nf (uncurry TS.gcdExt) (TS.toF2Poly $ randomVec TS.Bit k, TS.toF2Poly $ randomVec2 TS.Bit k) ] diff --git a/bench/Bench/Intersection.hs b/bench/Bench/Intersection.hs --- a/bench/Bench/Intersection.hs +++ b/bench/Bench/Intersection.hs @@ -1,3 +1,6 @@ +{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} +{-# HLINT ignore "Avoid lambda" #-} + module Bench.Intersection ( benchIntersection ) where @@ -7,53 +10,14 @@ import Data.Bits import qualified Data.IntSet as IS import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomBools :: [Bool] -randomBools - = map (> (0 :: Int)) - . randoms - . mkStdGen - $ 42 - -randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec f k = U.fromList (map f (take n randomBools)) - where - n = 1 `shiftL` k - -randomVec2 :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec2 f k = U.fromList (map f (take n $ drop n randomBools)) - where - n = 1 `shiftL` k - -randomSet :: Int -> IS.IntSet -randomSet k = IS.fromAscList (map fst (filter snd (zip [0..] (take n randomBools)))) - where - n = 1 `shiftL` k - -randomSet2 :: Int -> IS.IntSet -randomSet2 k = IS.fromAscList (map fst (filter snd (zip [0..] (take n $ drop n randomBools)))) - where - n = 1 `shiftL` k +import Bench.Common benchIntersection :: Int -> Benchmark benchIntersection k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf (\x -> intersectionBit (randomVec Bit k) x) (randomVec2 Bit k) - , bench "BitTS" $ nf (\x -> intersectionBitTS (randomVec TS.Bit k) x) (randomVec2 TS.Bit k) - , bench "Vector" $ nf (\x -> intersectionVector (randomVec id k) x) (randomVec2 id k) - , bench "IntSet" $ nf (intersectionIntSet (randomSet k)) (randomSet2 k) + [ bench labelBit $ nf (\x -> zipBits (.&.) (randomVec Bit k) x) (randomVec2 Bit k) + , bench labelBitTS $ nf (\x -> TS.zipBits (.&.) (randomVec TS.Bit k) x) (randomVec2 TS.Bit k) + , bench labelVector $ nf (\x -> U.zipWith (&&) (randomVec id k) x) (randomVec2 id k) + , bench labelIntSet $ nf (IS.union (randomSet k)) (randomSet2 k) ] - -intersectionBit :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit -intersectionBit = zipBits (.&.) - -intersectionBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit -> U.Vector TS.Bit -intersectionBitTS = TS.zipBits (.&.) - -intersectionVector :: U.Vector Bool -> U.Vector Bool -> U.Vector Bool -intersectionVector = U.zipWith (&&) - -intersectionIntSet :: IS.IntSet -> IS.IntSet -> IS.IntSet -intersectionIntSet = IS.union diff --git a/bench/Bench/Invert.hs b/bench/Bench/Invert.hs --- a/bench/Bench/Invert.hs +++ b/bench/Bench/Invert.hs @@ -6,34 +6,13 @@ import qualified Data.Bit.ThreadSafe as TS import Data.Bits import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomBools :: [Bool] -randomBools - = map (> (0 :: Int)) - . randoms - . mkStdGen - $ 42 - -randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec f k = U.fromList (map f (take n randomBools)) - where - n = 1 `shiftL` k +import Bench.Common benchInvert :: Int -> Benchmark benchInvert k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf invertBit (randomVec Bit k) - , bench "BitTS" $ nf invertBitTS (randomVec TS.Bit k) - , bench "Vector" $ nf invertVector (randomVec id k) + [ bench labelBit $ nf invertBits (randomVec Bit k) + , bench labelBitTS $ nf TS.invertBits (randomVec TS.Bit k) + , bench labelVector $ nf (U.map not) (randomVec id k) ] - -invertBit :: U.Vector Bit -> U.Vector Bit -invertBit = invertBits - -invertBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit -invertBitTS = TS.invertBits - -invertVector :: U.Vector Bool -> U.Vector Bool -invertVector = U.map not diff --git a/bench/Bench/Product.hs b/bench/Bench/Product.hs --- a/bench/Bench/Product.hs +++ b/bench/Bench/Product.hs @@ -1,3 +1,6 @@ +{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} +{-# HLINT ignore "Avoid lambda" #-} + module Bench.Product ( benchProduct , benchProductShort @@ -8,52 +11,29 @@ import qualified Data.Bit.ThreadSafe as TS import Data.Bits import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomBools :: [Bool] -randomBools - = map (> (0 :: Int)) - . randoms - . mkStdGen - $ 42 - -randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec f k = U.fromList (map f (take n randomBools)) - where - n = 1 `shiftL` k - -randomVec2 :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec2 f k = U.fromList (map f (take n $ drop n randomBools)) - where - n = 1 `shiftL` k - -randomInteger :: Int -> Integer -randomInteger k = toInteger $ toF2Poly $ randomVec Bit k - -randomInteger2 :: Int -> Integer -randomInteger2 k = toInteger $ toF2Poly $ randomVec2 Bit k +import Bench.Common benchProduct :: Int -> Benchmark benchProduct k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf (\x -> (*) (toF2Poly $ randomVec Bit k) x) (toF2Poly $ randomVec2 Bit k) - , bench "BitTS" $ nf (\x -> (*) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec2 TS.Bit k) - , bench "Integer" $ nf (\x -> binMul (randomInteger k) x) (randomInteger2 k) + [ bench labelBit $ nf (\x -> (*) (toF2Poly $ randomVec Bit k) x) (toF2Poly $ randomVec2 Bit k) + , bench labelBitTS $ nf (\x -> (*) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec2 TS.Bit k) + , bench labelInteger $ nf (\x -> binMul (randomInteger k) x) (randomInteger2 k) ] benchProductShort :: Int -> Benchmark benchProductShort k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf (\x -> (*) (toF2Poly $ randomVec Bit k) x) (toF2Poly $ U.take 32 $ randomVec2 Bit k) - , bench "BitTS" $ nf (\x -> (*) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ U.take 32 $ randomVec2 TS.Bit k) - , bench "Integer" $ nf (\x -> binMul (randomInteger k) x) ((1 `shiftL` 32 - 1) .&. randomInteger2 k) + [ bench labelBit $ nf (\x -> (*) (toF2Poly $ randomVec Bit k) x) (toF2Poly $ U.take 32 $ randomVec2 Bit k) + , bench labelBitTS $ nf (\x -> (*) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ U.take 32 $ randomVec2 TS.Bit k) + , bench labelInteger $ nf (\x -> binMul (randomInteger k) x) ((1 `shiftL` 32 - 1) .&. randomInteger2 k) ] benchSquare :: Int -> Benchmark benchSquare k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf (\x -> (*) (toF2Poly $ randomVec Bit k) x) (toF2Poly $ randomVec Bit k) - , bench "BitTS" $ nf (\x -> (*) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec TS.Bit k) - , bench "Integer" $ nf (\x -> binMul (randomInteger k) x) (randomInteger k) + [ bench labelBit $ nf (\x -> (*) (toF2Poly $ randomVec Bit k) x) (toF2Poly $ randomVec Bit k) + , bench labelBitTS $ nf (\x -> (*) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec TS.Bit k) + , bench labelInteger $ nf (\x -> binMul (randomInteger k) x) (randomInteger k) ] binMul :: Integer -> Integer -> Integer diff --git a/bench/Bench/RandomFlip.hs b/bench/Bench/RandomFlip.hs --- a/bench/Bench/RandomFlip.hs +++ b/bench/Bench/RandomFlip.hs @@ -10,29 +10,23 @@ import Data.Foldable import qualified Data.IntSet as IS import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomFlips :: [Int] -randomFlips - = map abs - . randoms - . mkStdGen - $ 42 +import Bench.Common benchRandomFlip :: Int -> Benchmark benchRandomFlip k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf randomFlipBit k - , bench "BitTS" $ nf randomFlipBitTS k - , bench "Vector" $ nf randomFlipVector k - , bench "IntSet" $ nf randomFlipIntSet k + [ bench labelBit $ nf randomFlipBit k + , bench labelBitTS $ nf randomFlipBitTS k + , bench labelVector $ nf randomFlipVector k + , bench labelIntSet $ nf randomFlipIntSet k ] randomFlipBit :: Int -> Int randomFlipBit k = runST $ do let n = 1 `shiftL` k vec <- MU.new n - forM_ (take (mult * n) randomFlips) $ + forM_ (take (mult * n) randomIndices) $ \i -> unsafeFlipBit vec (i .&. (1 `shiftL` k - 1)) Bit i <- MU.unsafeRead vec 0 pure $ if i then 1 else 0 @@ -41,7 +35,7 @@ randomFlipBitTS k = runST $ do let n = 1 `shiftL` k vec <- MU.new n - forM_ (take (mult * n) randomFlips) $ + forM_ (take (mult * n) randomIndices) $ \i -> TS.unsafeFlipBit vec (i .&. (1 `shiftL` k - 1)) TS.Bit i <- MU.unsafeRead vec 0 pure $ if i then 1 else 0 @@ -50,7 +44,7 @@ randomFlipVector k = runST $ do let n = 1 `shiftL` k vec <- MU.new n - forM_ (take (mult * n) randomFlips) $ + forM_ (take (mult * n) randomIndices) $ \i -> MU.unsafeModify vec complement (i .&. (1 `shiftL` k - 1)) i <- MU.unsafeRead vec 0 pure $ if i then 1 else 0 @@ -62,7 +56,7 @@ vec = foldl' (\acc i -> let j = i .&. (1 `shiftL` k - 1) in (if IS.member j acc then IS.delete else IS.insert) j acc) mempty - (take (mult * n) randomFlips) + (take (mult * n) randomIndices) mult :: Int mult = 100 diff --git a/bench/Bench/RandomRead.hs b/bench/Bench/RandomRead.hs --- a/bench/Bench/RandomRead.hs +++ b/bench/Bench/RandomRead.hs @@ -10,67 +10,54 @@ -- import Data.List import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomVec :: [Bool] -randomVec - = map (> (0 :: Int)) - . randoms - . mkStdGen - $ 42 - -randomReads :: [Int] -randomReads - = map abs - . randoms - . mkStdGen - $ 42 +import Bench.Common benchRandomRead :: Int -> Benchmark benchRandomRead k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf randomReadBit k - , bench "BitTS" $ nf randomReadBitTS k - , bench "Vector" $ nf randomReadVector k - -- , bench "IntSet" $ nf randomReadIntSet k + [ bench labelBit $ nf randomReadBit k + , bench labelBitTS $ nf randomReadBitTS k + , bench labelVector $ nf randomReadVector k + -- , bench labelIntSet $ nf randomReadIntSet k ] randomReadBit :: Int -> Int randomReadBit k = runST $ do let n = 1 `shiftL` k - vec <- U.unsafeThaw (U.fromList (map Bit $ take n randomVec)) + vec <- U.unsafeThaw (U.fromList (map Bit $ take n randomBools)) let go acc [] = pure acc go acc (i : is) = do Bit b <- MU.unsafeRead vec (i .&. (1 `shiftL` k - 1)) go (acc + if b then 1 else 0) is - go 0 (take (mult * n) randomReads) + go 0 (take (mult * n) randomIndices) randomReadBitTS :: Int -> Int randomReadBitTS k = runST $ do let n = 1 `shiftL` k - vec <- U.unsafeThaw (U.fromList (map TS.Bit $ take n randomVec)) + vec <- U.unsafeThaw (U.fromList (map TS.Bit $ take n randomBools)) let go acc [] = pure acc go acc (i : is) = do TS.Bit b <- MU.unsafeRead vec (i .&. (1 `shiftL` k - 1)) go (acc + if b then 1 else 0) is - go 0 (take (mult * n) randomReads) + go 0 (take (mult * n) randomIndices) randomReadVector :: Int -> Int randomReadVector k = runST $ do let n = 1 `shiftL` k - vec <- U.unsafeThaw (U.fromList (take n randomVec)) + vec <- U.unsafeThaw (U.fromList (take n randomBools)) let go acc [] = pure acc go acc (i : is) = do b <- MU.unsafeRead vec (i .&. (1 `shiftL` k - 1)) go (acc + if b then 1 else 0) is - go 0 (take (mult * n) randomReads) + go 0 (take (mult * n) randomIndices) -- randomReadIntSet :: Int -> Int --- randomReadIntSet k = foldl' (+) 0 [ doRead (c + i `shiftL` 1 - i - c) | c <- [0 .. mult - 1], i <- randomReads ] +-- randomReadIntSet k = foldl' (+) 0 [ doRead (c + i `shiftL` 1 - i - c) | c <- [0 .. mult - 1], i <- randomIndices ] -- where -- n = 1 `shiftL` k -- vec = IS.fromDistinctAscList $ map fst $ filter snd --- $ zip [0..] $ take n randomVec +-- $ zip [0..] $ take n randomBools -- doRead i = if IS.member (i .&. (1 `shiftL` k - 1)) vec then 1 else 0 mult :: Int diff --git a/bench/Bench/RandomWrite.hs b/bench/Bench/RandomWrite.hs --- a/bench/Bench/RandomWrite.hs +++ b/bench/Bench/RandomWrite.hs @@ -10,29 +10,23 @@ import Data.Foldable import qualified Data.IntSet as IS import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomWrites :: [(Int, Bool)] -randomWrites - = map (\x -> if x > 0 then (x, True) else (negate x, False)) - . randoms - . mkStdGen - $ 42 +import Bench.Common benchRandomWrite :: Int -> Benchmark benchRandomWrite k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf randomWriteBit k - , bench "BitTS" $ nf randomWriteBitTS k - , bench "Vector" $ nf randomWriteVector k - , bench "IntSet" $ nf randomWriteIntSet k + [ bench labelBit $ nf randomWriteBit k + , bench labelBitTS $ nf randomWriteBitTS k + , bench labelVector $ nf randomWriteVector k + , bench labelIntSet $ nf randomWriteIntSet k ] randomWriteBit :: Int -> Int randomWriteBit k = runST $ do let n = 1 `shiftL` k vec <- MU.new n - forM_ (take (mult * n) randomWrites) $ + forM_ (take (mult * n) randomIndicesAndBools) $ \(i, b) -> MU.unsafeWrite vec (i .&. (1 `shiftL` k - 1)) (Bit b) Bit i <- MU.unsafeRead vec 0 pure $ if i then 1 else 0 @@ -41,7 +35,7 @@ randomWriteBitTS k = runST $ do let n = 1 `shiftL` k vec <- MU.new n - forM_ (take (mult * n) randomWrites) $ + forM_ (take (mult * n) randomIndicesAndBools) $ \(i, b) -> MU.unsafeWrite vec (i .&. (1 `shiftL` k - 1)) (TS.Bit b) TS.Bit i <- MU.unsafeRead vec 0 pure $ if i then 1 else 0 @@ -50,7 +44,7 @@ randomWriteVector k = runST $ do let n = 1 `shiftL` k vec <- MU.new n - forM_ (take (mult * n) randomWrites) $ + forM_ (take (mult * n) randomIndicesAndBools) $ \(i, b) -> MU.unsafeWrite vec (i .&. (1 `shiftL` k - 1)) b i <- MU.unsafeRead vec 0 pure $ if i then 1 else 0 @@ -62,7 +56,7 @@ vec = foldl' (\acc (i, b) -> (if b then IS.insert else IS.delete) (i .&. (1 `shiftL` k - 1)) acc) mempty - (take (mult * n) randomWrites) + (take (mult * n) randomIndicesAndBools) mult :: Int mult = 100 diff --git a/bench/Bench/Remainder.hs b/bench/Bench/Remainder.hs --- a/bench/Bench/Remainder.hs +++ b/bench/Bench/Remainder.hs @@ -1,5 +1,7 @@ {-# LANGUAGE CPP #-} {-# LANGUAGE MagicHash #-} +{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} +{-# HLINT ignore "Avoid lambda" #-} module Bench.Remainder ( benchRemainder @@ -8,45 +10,21 @@ import Data.Bit import qualified Data.Bit.ThreadSafe as TS import Data.Bits -import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Unboxed.Mutable as MU import GHC.Exts #ifdef MIN_VERSION_ghc_bignum import GHC.Num.Integer #else import GHC.Integer.Logarithms #endif -import System.Random import Test.Tasty.Bench -randomBools :: [Bool] -randomBools - = map (> (0 :: Int)) - . randoms - . mkStdGen - $ 42 - -randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec f k = U.fromList (map f (take (2 * n) randomBools)) - where - n = 1 `shiftL` k - -randomVec2 :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec2 f k = U.fromList (map f (take n $ drop (2 * n) randomBools)) - where - n = 1 `shiftL` k - -randomInteger :: Int -> Integer -randomInteger k = toInteger $ toF2Poly $ randomVec Bit k - -randomInteger2 :: Int -> Integer -randomInteger2 k = toInteger $ toF2Poly $ randomVec2 Bit k +import Bench.Common benchRemainder :: Int -> Benchmark benchRemainder k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf (\x -> rem (toF2Poly $ randomVec Bit k) x) (toF2Poly $ randomVec2 Bit k) - , bench "BitTS" $ nf (\x -> rem (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec2 TS.Bit k) - , bench "Integer" $ nf (\x -> binRem (randomInteger k) x) (randomInteger2 k) + [ bench labelBit $ nf (\x -> rem (toF2Poly $ randomVec Bit (2 * k)) x) (toF2Poly $ randomVec2 Bit k) + , bench labelBitTS $ nf (\x -> rem (TS.toF2Poly $ randomVec TS.Bit (2 * k)) x) (TS.toF2Poly $ randomVec2 TS.Bit k) + , bench labelInteger $ nf (\x -> binRem (randomInteger (2 * k)) x) (randomInteger2 k) ] binRem :: Integer -> Integer -> Integer diff --git a/bench/Bench/Reverse.hs b/bench/Bench/Reverse.hs --- a/bench/Bench/Reverse.hs +++ b/bench/Bench/Reverse.hs @@ -6,34 +6,13 @@ import qualified Data.Bit.ThreadSafe as TS import Data.Bits import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomBools :: [Bool] -randomBools - = map (> (0 :: Int)) - . randoms - . mkStdGen - $ 42 - -randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec f k = U.fromList (map f (take n randomBools)) - where - n = 1 `shiftL` k +import Bench.Common benchReverse :: Int -> Benchmark benchReverse k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf reverseBit (randomVec Bit k) - , bench "BitTS" $ nf reverseBitTS (randomVec TS.Bit k) - , bench "Vector" $ nf reverseVector (randomVec id k) + [ bench labelBit $ nf reverseBits (randomVec Bit k) + , bench labelBitTS $ nf TS.reverseBits (randomVec TS.Bit k) + , bench labelVector $ nf U.reverse (randomVec id k) ] - -reverseBit :: U.Vector Bit -> U.Vector Bit -reverseBit = reverseBits - -reverseBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit -reverseBitTS = TS.reverseBits - -reverseVector :: U.Vector Bool -> U.Vector Bool -reverseVector = U.reverse diff --git a/bench/Bench/Sum.hs b/bench/Bench/Sum.hs --- a/bench/Bench/Sum.hs +++ b/bench/Bench/Sum.hs @@ -1,3 +1,6 @@ +{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} +{-# HLINT ignore "Avoid lambda" #-} + module Bench.Sum ( benchAdd , benchSum @@ -7,44 +10,20 @@ import qualified Data.Bit.ThreadSafe as TS import Data.Bits import Data.Foldable -import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomBools :: [Bool] -randomBools - = map (> (0 :: Int)) - . randoms - . mkStdGen - $ 42 - -randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec f k = U.fromList (map f (take n randomBools)) - where - n = 1 `shiftL` k - -randomVec2 :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec2 f k = U.fromList (map f (take n $ drop n randomBools)) - where - n = 1 `shiftL` k - -randomInteger :: Int -> Integer -randomInteger k = toInteger $ toF2Poly $ randomVec Bit k - -randomInteger2 :: Int -> Integer -randomInteger2 k = toInteger $ toF2Poly $ randomVec2 Bit k +import Bench.Common benchAdd :: Int -> Benchmark benchAdd k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf (\x -> (+) (toF2Poly $ randomVec Bit k) x) (toF2Poly $ randomVec2 Bit k) - , bench "BitTS" $ nf (\x -> (+) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec2 TS.Bit k) - , bench "Integer" $ nf (\x -> xor (randomInteger k) x) (randomInteger2 k) + [ bench labelBit $ nf (\x -> (+) (toF2Poly $ randomVec Bit k) x) (toF2Poly $ randomVec2 Bit k) + , bench labelBitTS $ nf (\x -> (+) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec2 TS.Bit k) + , bench labelInteger $ nf (\x -> xor (randomInteger k) x) (randomInteger2 k) ] benchSum :: Int -> Benchmark benchSum k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf (foldl' (+) 0) [(1 :: F2Poly) .. fromInteger (1 `shiftL` k)] - , bench "BitTS" $ nf (foldl' (+) 0) [(1 :: TS.F2Poly) .. fromInteger (1 `shiftL` k)] - , bench "Integer" $ nf (foldl' xor 0) [(1 :: Integer) .. fromInteger (1 `shiftL` k)] + [ bench labelBit $ nf (foldl' (+) 0) [(1 :: F2Poly) .. fromInteger (1 `shiftL` k)] + , bench labelBitTS $ nf (foldl' (+) 0) [(1 :: TS.F2Poly) .. fromInteger (1 `shiftL` k)] + , bench labelInteger $ nf (foldl' xor 0) [(1 :: Integer) .. fromInteger (1 `shiftL` k)] ] diff --git a/bench/Bench/Union.hs b/bench/Bench/Union.hs --- a/bench/Bench/Union.hs +++ b/bench/Bench/Union.hs @@ -1,3 +1,6 @@ +{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} +{-# HLINT ignore "Avoid lambda" #-} + module Bench.Union ( benchUnion ) where @@ -7,53 +10,14 @@ import Data.Bits import qualified Data.IntSet as IS import qualified Data.Vector.Unboxed as U -import qualified Data.Vector.Unboxed.Mutable as MU -import System.Random import Test.Tasty.Bench -randomBools :: [Bool] -randomBools - = map (> (0 :: Int)) - . randoms - . mkStdGen - $ 42 - -randomVec :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec f k = U.fromList (map f (take n randomBools)) - where - n = 1 `shiftL` k - -randomVec2 :: MU.Unbox a => (Bool -> a) -> Int -> U.Vector a -randomVec2 f k = U.fromList (map f (take n $ drop n randomBools)) - where - n = 1 `shiftL` k - -randomSet :: Int -> IS.IntSet -randomSet k = IS.fromAscList (map fst (filter snd (zip [0..] (take n randomBools)))) - where - n = 1 `shiftL` k - -randomSet2 :: Int -> IS.IntSet -randomSet2 k = IS.fromAscList (map fst (filter snd (zip [0..] (take n $ drop n randomBools)))) - where - n = 1 `shiftL` k +import Bench.Common benchUnion :: Int -> Benchmark benchUnion k = bgroup (show (1 `shiftL` k :: Int)) - [ bench "Bit" $ nf (\x -> unionBit (randomVec Bit k) x) (randomVec2 Bit k) - , bench "BitTS" $ nf (\x -> unionBitTS (randomVec TS.Bit k) x) (randomVec2 TS.Bit k) - , bench "Vector" $ nf (\x -> unionVector (randomVec id k) x) (randomVec2 id k) - , bench "IntSet" $ nf (unionIntSet (randomSet k)) (randomSet2 k) + [ bench labelBit $ nf (\x -> zipBits (.|.) (randomVec Bit k) x) (randomVec2 Bit k) + , bench labelBitTS $ nf (\x -> TS.zipBits (.|.) (randomVec TS.Bit k) x) (randomVec2 TS.Bit k) + , bench labelVector $ nf (\x -> U.zipWith (||) (randomVec id k) x) (randomVec2 id k) + , bench labelIntSet $ nf (IS.union (randomSet k)) (randomSet2 k) ] - -unionBit :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit -unionBit = zipBits (.|.) - -unionBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit -> U.Vector TS.Bit -unionBitTS = TS.zipBits (.|.) - -unionVector :: U.Vector Bool -> U.Vector Bool -> U.Vector Bool -unionVector = U.zipWith (||) - -unionIntSet :: IS.IntSet -> IS.IntSet -> IS.IntSet -unionIntSet = IS.union diff --git a/bitvec.cabal b/bitvec.cabal --- a/bitvec.cabal +++ b/bitvec.cabal @@ -1,5 +1,5 @@ name: bitvec -version: 1.1.4.0 +version: 1.1.5.0 cabal-version: 2.0 build-type: Simple license: BSD3 @@ -9,7 +9,7 @@ homepage: https://github.com/Bodigrim/bitvec synopsis: Space-efficient bit vectors description: - A newtype over 'Bool' with a better 'Vector' instance: 8x less memory, up to 1000x faster. + A newtype over 'Bool' with a better 'Vector' instance: 8x less memory, up to 3500x faster. . The package represents unboxed arrays of 'Bool's @@ -21,16 +21,21 @@ the most significant degradation happens for random writes (up to 10% slower). On the other hand, for certain bulk bit operations - 'Vector' 'Bit' is up to 1000x faster than 'Vector' 'Bool'. + 'Vector' 'Bit' is up to 3500x faster than 'Vector' 'Bool'. . === Thread safety . - * "Data.Bit" is faster, but writes and flips are thread-unsafe. + * "Data.Bit" is faster, but writes and flips are not thread-safe. This is because naive updates are not atomic: they read the whole word from memory, then modify a bit, then write the whole word back. + Concurrently modifying non-intersecting slices of the same underlying array + may also lead to unexpected results, since they can share a word in memory. * "Data.Bit.ThreadSafe" is slower (usually 10-20%), but writes and flips are thread-safe. + Additionally, concurrently modifying non-intersecting slices of the same underlying array + works as expected. However, operations that affect multiple elements are not + guaranteed to be atomic. . === Similar packages . @@ -55,12 +60,11 @@ type: git location: git://github.com/Bodigrim/bitvec.git -flag libgmp +flag simd description: - Link against the GMP library for the ultimate performance of - `zipBits`, `invertBits` and `countBits`. Users are strongly encouraged - to enable this flag whenever possible. - default: False + Use a C SIMD implementation for the ultimate performance of `zipBits`, `invertBits` and `countBits`. + Disable this flag if there are problems with the C FFI. + default: True manual: True library @@ -69,8 +73,8 @@ Data.Bit.ThreadSafe build-depends: base >=4.11 && <5, - bytestring >=0.10 && <0.12, - deepseq <1.5, + bytestring >=0.10 && <0.13, + deepseq <1.6, primitive >=0.5, vector >=0.11 && <0.14 default-language: Haskell2010 @@ -94,11 +98,12 @@ else build-depends: ghc-bignum - if flag(libgmp) - extra-libraries: gmp - cpp-options: -DUseLibGmp + if flag(simd) && !arch(javascript) && !arch(wasm32) + c-sources: cbits/bitvec_simd.c + cc-options: -fopenmp-simd + cpp-options: -DUseSIMD other-modules: - Data.Bit.Gmp + Data.Bit.SIMD test-suite bitvec-tests type: exitcode-stdio-1.0 @@ -146,6 +151,7 @@ hs-source-dirs: bench other-modules: Bench.BitIndex + Bench.Common Bench.GCD Bench.Invert Bench.Intersection diff --git a/cbits/bitvec_simd.c b/cbits/bitvec_simd.c new file mode 100644 --- /dev/null +++ b/cbits/bitvec_simd.c @@ -0,0 +1,450 @@ +#include +#include + +#ifdef __x86_64__ +#include +#endif + +#include "HsFFI.h" + +HsInt _hs_bitvec_popcount(const uint32_t *src, HsInt len) { + HsInt count = 0; + #pragma omp simd + for (size_t i = 0; i < len; i++) { + uint32_t x = src[i]; + // count += popcount(t); + // https://bits.stephan-brumme.com/countBits.html + x = x - ((x >> 1) & 0x55555555); + x = (x & 0x33333333) + ((x >> 2) & 0x33333333); + x = (x + (x >> 4)) & 0x0f0f0f0f; + count += (x * 0x01010101) >> 24; + } + return count; +} + +void _hs_bitvec_com(uint8_t *dest, uint8_t *src, HsInt len) { + #pragma omp simd + for (size_t i = 0; i < len; i++) { + dest[i] = ~src[i]; + } +} + +void _hs_bitvec_and(uint8_t *dest, const uint8_t *src1, const uint8_t *src2, HsInt len) { + #pragma omp simd + for (size_t i = 0; i < len; i++) { + dest[i] = src1[i] & src2[i]; + } +} + +void _hs_bitvec_ior(uint8_t *dest, const uint8_t *src1, const uint8_t *src2, HsInt len) { + #pragma omp simd + for (size_t i = 0; i < len; i++) { + dest[i] = src1[i] | src2[i]; + } +} + +void _hs_bitvec_xor(uint8_t *dest, const uint8_t *src1, const uint8_t *src2, HsInt len) { + #pragma omp simd + for (size_t i = 0; i < len; i++) { + dest[i] = src1[i] ^ src2[i]; + } +} + +void _hs_bitvec_andn(uint8_t *dest, const uint8_t *src1, const uint8_t *src2, HsInt len) { + #pragma omp simd + for (size_t i = 0; i < len; i++) { + dest[i] = src1[i] & (~src2[i]); + } +} + +void _hs_bitvec_iorn(uint8_t *dest, const uint8_t *src1, const uint8_t *src2, HsInt len) { + #pragma omp simd + for (size_t i = 0; i < len; i++) { + dest[i] = src1[i] | (~src2[i]); + } +} + +void _hs_bitvec_nand(uint8_t *dest, const uint8_t *src1, const uint8_t *src2, HsInt len) { + #pragma omp simd + for (size_t i = 0; i < len; i++) { + dest[i] = ~(src1[i] & src2[i]); + } +} + +void _hs_bitvec_nior(uint8_t *dest, const uint8_t *src1, const uint8_t *src2, HsInt len) { + #pragma omp simd + for (size_t i = 0; i < len; i++) { + dest[i] = ~(src1[i] | src2[i]); + } +} + +void _hs_bitvec_xnor(uint8_t *dest, const uint8_t *src1, const uint8_t *src2, HsInt len) { + #pragma omp simd + for (size_t i = 0; i < len; i++) { + dest[i] = ~(src1[i] ^ src2[i]); + } +} + + +#ifdef __x86_64__ +static void reverse_bits_sse(uint64_t *dest, const uint64_t *src, HsInt len) { + __m128i mask1l = _mm_set1_epi32(0x55555555); + __m128i mask1r = _mm_set1_epi32(0xaaaaaaaa); + __m128i mask2l = _mm_set1_epi32(0x33333333); + __m128i mask2r = _mm_set1_epi32(0xcccccccc); + __m128i mask4l = _mm_set1_epi32(0x0f0f0f0f); + __m128i mask4r = _mm_set1_epi32(0xf0f0f0f0); + __m128i mask8l = _mm_set1_epi32(0x00ff00ff); + __m128i mask8r = _mm_set1_epi32(0xff00ff00); + __m128i mask16l = _mm_set1_epi32(0x0000ffff); + __m128i mask16r = _mm_set1_epi32(0xffff0000); + + size_t i = 0; + for (; i < (len & (~0x1)); i += 2) { + __m128i x = _mm_loadu_si128((const __m128i *) (src + i)); + + // reverse each word + x = _mm_or_si128(_mm_slli_epi32(_mm_and_si128(x, mask1l), 1), _mm_srli_epi32(_mm_and_si128(x, mask1r), 1)); + x = _mm_or_si128(_mm_slli_epi32(_mm_and_si128(x, mask2l), 2), _mm_srli_epi32(_mm_and_si128(x, mask2r), 2)); + x = _mm_or_si128(_mm_slli_epi32(_mm_and_si128(x, mask4l), 4), _mm_srli_epi32(_mm_and_si128(x, mask4r), 4)); + x = _mm_or_si128(_mm_slli_epi32(_mm_and_si128(x, mask8l), 8), _mm_srli_epi32(_mm_and_si128(x, mask8r), 8)); + x = _mm_or_si128(_mm_slli_epi32(_mm_and_si128(x, mask16l), 16), _mm_srli_epi32(_mm_and_si128(x, mask16r), 16)); + + // reverse order of words + x = _mm_shuffle_epi32(x, 0x1b); + + _mm_storeu_si128((__m128i *) (dest + len - 2 - i), x); + } + for (; i < len; i++) { + uint64_t x = src[i]; + x = ((x & 0x5555555555555555) << 1) | ((x & 0xaaaaaaaaaaaaaaaa) >> 1); + x = ((x & 0x3333333333333333) << 2) | ((x & 0xcccccccccccccccc) >> 2); + x = ((x & 0x0f0f0f0f0f0f0f0f) << 4) | ((x & 0xf0f0f0f0f0f0f0f0) >> 4); + x = ((x & 0x00ff00ff00ff00ff) << 8) | ((x & 0xff00ff00ff00ff00) >> 8); + x = ((x & 0x0000ffff0000ffff) << 16) | ((x & 0xffff0000ffff0000) >> 16); + x = ((x & 0x00000000ffffffff) << 32) | ((x & 0xffffffff00000000) >> 32); + dest[len - 1 - i] = x; + } +} + +__attribute__((target("avx2"))) +static void reverse_bits_avx(uint64_t *dest, const uint64_t *src, HsInt len) { + __m256i mask1l = _mm256_set1_epi32(0x55555555); + __m256i mask1r = _mm256_set1_epi32(0xaaaaaaaa); + __m256i mask2l = _mm256_set1_epi32(0x33333333); + __m256i mask2r = _mm256_set1_epi32(0xcccccccc); + __m256i mask4l = _mm256_set1_epi32(0x0f0f0f0f); + __m256i mask4r = _mm256_set1_epi32(0xf0f0f0f0); + __m256i mask8l = _mm256_set1_epi32(0x00ff00ff); + __m256i mask8r = _mm256_set1_epi32(0xff00ff00); + __m256i mask16l = _mm256_set1_epi32(0x0000ffff); + __m256i mask16r = _mm256_set1_epi32(0xffff0000); + + size_t i = 0; + for (; i < (len & (~0x3)); i += 4) { + __m256i x = _mm256_loadu_si256((const __m256i *) (src + i)); + + // reverse each word + x = _mm256_or_si256(_mm256_slli_epi32(_mm256_and_si256(x, mask1l), 1), _mm256_srli_epi32(_mm256_and_si256(x, mask1r), 1)); + x = _mm256_or_si256(_mm256_slli_epi32(_mm256_and_si256(x, mask2l), 2), _mm256_srli_epi32(_mm256_and_si256(x, mask2r), 2)); + x = _mm256_or_si256(_mm256_slli_epi32(_mm256_and_si256(x, mask4l), 4), _mm256_srli_epi32(_mm256_and_si256(x, mask4r), 4)); + x = _mm256_or_si256(_mm256_slli_epi32(_mm256_and_si256(x, mask8l), 8), _mm256_srli_epi32(_mm256_and_si256(x, mask8r), 8)); + x = _mm256_or_si256(_mm256_slli_epi32(_mm256_and_si256(x, mask16l), 16), _mm256_srli_epi32(_mm256_and_si256(x, mask16r), 16)); + + // reverse order of words + x = _mm256_permutevar8x32_epi32(x, _mm256_setr_epi32(7, 6, 5, 4, 3, 2, 1, 0)); + + _mm256_storeu_si256((__m256i *) (dest + len - 4 - i), x); + } + for (; i < len; i++) { + uint64_t x = src[i]; + x = ((x & 0x5555555555555555) << 1) | ((x & 0xaaaaaaaaaaaaaaaa) >> 1); + x = ((x & 0x3333333333333333) << 2) | ((x & 0xcccccccccccccccc) >> 2); + x = ((x & 0x0f0f0f0f0f0f0f0f) << 4) | ((x & 0xf0f0f0f0f0f0f0f0) >> 4); + x = ((x & 0x00ff00ff00ff00ff) << 8) | ((x & 0xff00ff00ff00ff00) >> 8); + x = ((x & 0x0000ffff0000ffff) << 16) | ((x & 0xffff0000ffff0000) >> 16); + x = ((x & 0x00000000ffffffff) << 32) | ((x & 0xffffffff00000000) >> 32); + dest[len - 1 - i] = x; + } +} +#endif + +void _hs_bitvec_reverse_bits(HsWord *dest, const HsWord *src, HsInt len) { +#ifdef __x86_64__ + if (__builtin_cpu_supports("avx2")) { + reverse_bits_avx(dest, src, len); + } else { + reverse_bits_sse(dest, src, len); + } +#else + if (sizeof(HsWord) == 8) { + // 64 bit + for (size_t i = 0; i < len; i++) { + uint64_t x = src[i]; + x = ((x & 0x5555555555555555) << 1) | ((x & 0xaaaaaaaaaaaaaaaa) >> 1); + x = ((x & 0x3333333333333333) << 2) | ((x & 0xcccccccccccccccc) >> 2); + x = ((x & 0x0f0f0f0f0f0f0f0f) << 4) | ((x & 0xf0f0f0f0f0f0f0f0) >> 4); + x = ((x & 0x00ff00ff00ff00ff) << 8) | ((x & 0xff00ff00ff00ff00) >> 8); + x = ((x & 0x0000ffff0000ffff) << 16) | ((x & 0xffff0000ffff0000) >> 16); + x = ((x & 0x00000000ffffffff) << 32) | ((x & 0xffffffff00000000) >> 32); + dest[len - 1 - i] = x; + } + } else { + // 32 bit + for (size_t i = 0; i < len; i++) { + uint32_t x = src[i]; + x = ((x & 0x55555555) << 1) | ((x & 0xaaaaaaaa) >> 1); + x = ((x & 0x33333333) << 2) | ((x & 0xcccccccc) >> 2); + x = ((x & 0x0f0f0f0f) << 4) | ((x & 0xf0f0f0f0) >> 4); + x = ((x & 0x00ff00ff) << 8) | ((x & 0xff00ff00) >> 8); + x = ((x & 0x0000ffff) << 16) | ((x & 0xffff0000) >> 16); + dest[len - 1 - i] = x; + } + } +#endif +} + + +#ifdef __x86_64__ +static HsInt bit_index_sse(const uint64_t *src, HsInt len, HsBool bit) { + __m128i zero = _mm_setzero_si128(); + __m128i bit_mask_128; + uint64_t bit_mask_64; + if (bit) { + bit_mask_128 = zero; + bit_mask_64 = 0; + } else { + bit_mask_128 = _mm_set1_epi64x(0xffffffffffffffff); + bit_mask_64 = 0xffffffffffffffff; + } + size_t i = 0; + for (; i < (len & (~0x1)); i += 2) { + __m128i x = _mm_xor_si128(_mm_loadu_si128((const __m128i *) (src + i)), bit_mask_128); + uint16_t mask = ~_mm_movemask_epi8(_mm_cmpeq_epi32(x, zero)); + if (mask != 0) { + size_t idx = __builtin_ctz(mask) >> 3; + uint64_t x = src[i + idx] ^ bit_mask_64; + return ((i + idx) << 6) + __builtin_ctzll(x); + } + } + for (; i < len; i++) { + uint64_t x = src[i] ^ bit_mask_64; + if (x != 0) { + return (i << 6) + __builtin_ctzll(x); + } + } + return -1; +} + +__attribute__((target("avx2"))) +static HsInt bit_index_avx(const uint64_t *src, HsInt len, HsBool bit) { + __m256i zero = _mm256_setzero_si256(); + __m256i bit_mask_256; + uint64_t bit_mask_64; + if (bit) { + bit_mask_256 = zero; + bit_mask_64 = 0; + } else { + bit_mask_256 = _mm256_set1_epi64x(0xffffffffffffffff); + bit_mask_64 = 0xffffffffffffffff; + } + size_t i = 0; + for (; i < (len & (~0x3)); i += 4) { + __m256i x = _mm256_xor_si256(_mm256_loadu_si256((const __m256i *) (src + i)), bit_mask_256); + uint32_t mask = ~_mm256_movemask_epi8(_mm256_cmpeq_epi32(x, zero)); + if (mask != 0) { + size_t idx = __builtin_ctzl(mask) >> 3; + uint64_t x = src[i + idx] ^ bit_mask_64; + return ((i + idx) << 6) + __builtin_ctzll(x); + } + } + for (; i < len; i++) { + uint64_t x = src[i] ^ bit_mask_64; + if (x != 0) { + return (i << 6) + __builtin_ctzll(x); + } + } + return -1; +} +#endif + +HsInt _hs_bitvec_bit_index(const HsWord *src, HsInt len, HsBool bit) { +#ifdef __x86_64__ + if (__builtin_cpu_supports("avx2")) { + return bit_index_avx(src, len, bit); + } else { + return bit_index_sse(src, len, bit); + } +#else + HsWord bit_mask; + if (bit) { + bit_mask = 0; + } else { + bit_mask = -1; + } + for (size_t i = 0; i < len; i++) { + HsWord x = src[i] ^ bit_mask; + if (x != 0) { + return (i << 3) * sizeof(HsWord) + __builtin_ctzll(x); + } + } + return -1; +#endif +} + + +#ifdef __x86_64__ +__attribute__((target("popcnt"))) +static HsInt nth_bit_index_popcnt(const uint64_t *src, HsInt len, HsBool bit, HsInt n) { + uint64_t bit_mask; + if (bit) { + bit_mask = 0; + } else { + bit_mask = -1; + } + for (size_t i = 0; i < len; i++) { + uint64_t x = src[i] ^ bit_mask; + + HsInt count = _mm_popcnt_u64(x); + if (n <= count) { + for (size_t i = 0; i < n - 1; i++) { + // clear lowest set bit + x &= x - 1; + } + return (i << 6) + __builtin_ctzll(x); + } else { + n -= count; + } + } + return -1; +} +#endif + +HsInt _hs_bitvec_nth_bit_index(const HsWord *src, HsInt len, HsBool bit, HsInt n) { +#ifdef __x86_64__ + if (__builtin_cpu_supports("popcnt")) { + return nth_bit_index_popcnt(src, len, bit, n); + } +#endif + HsWord bit_mask; + if (bit) { + bit_mask = 0; + } else { + bit_mask = -1; + } + for (size_t i = 0; i < len; i++) { + HsWord x = src[i] ^ bit_mask; + + // popcount + HsWord count = x - ((x >> 1) & 0x5555555555555555); + count = (count & 0x3333333333333333) + ((count >> 2) & 0x3333333333333333); + count = (count + (count >> 4)) & 0x0f0f0f0f0f0f0f0f; + count = (count * 0x101010101010101) >> 56; + + if (n <= count) { + for (size_t i = 0; i < n - 1; i++) { + // clear lowest set bit + x &= x - 1; + } + return (i << 3) * sizeof(HsWord) + __builtin_ctzll(x); + } else { + n -= count; + } + } + return -1; +} + + +#ifdef __x86_64__ +__attribute__((target("popcnt,bmi2"))) +static HsInt select_bits_pext(uint64_t *dest, const uint64_t *src, const uint64_t *mask, HsInt len, HsBool exclude) { + uint64_t bit_mask; + if (exclude) { + bit_mask = -1; + } else { + bit_mask = 0; + } + HsInt off = 0; // offset in bits into `dest` + for (size_t i = 0; i < len; i++) { + uint64_t x = src[i]; + uint64_t m = mask[i] ^ bit_mask; + HsInt count = _mm_popcnt_u64(m); + uint64_t y = _pext_u64(x, m); + HsInt off_words = off >> 6; + HsInt off_bits = off & 0x3f; + if (off_bits == 0) { + dest[off_words] = y; + } else { + dest[off_words] |= y << off_bits; + dest[off_words + 1] = y >> (64 - off_bits); + } + off += count; + } + return off; +} +#endif + +HsInt _hs_bitvec_select_bits(HsWord *dest, const HsWord *src, const HsWord *mask, HsInt len, HsBool exclude) { +#ifdef __x86_64__ + if (__builtin_cpu_supports("popcnt") && __builtin_cpu_supports("bmi2")) { + return select_bits_pext(dest, src, mask, len, exclude); + } +#endif + HsWord bit_mask; + if (exclude) { + bit_mask = -1; + } else { + bit_mask = 0; + } + HsInt off = 0; // offset in bits into `dest` + for (size_t i = 0; i < len; i++) { + HsWord x = src[i]; + HsWord m = mask[i] ^ bit_mask; + + // pext + HsWord y = 0; + HsInt count = 0; + if (m == -1) { + y = x; + count = sizeof(HsWord) * 8; + } else { + HsWord bb = 1; + for (; m != 0; bb <<= 1) { + if (x & m & -m) { + y |= bb; + } + m &= m - 1; + } + if (sizeof(HsWord) == 8) { + count = __builtin_ctzll(bb); + } else { + count = __builtin_ctzl(bb); + } + } + + if (sizeof(HsWord) == 8) { + // 64 bit + HsInt off_words = off >> 6; + HsInt off_bits = off & 0x3f; + if (off_bits == 0) { + dest[off_words] = y; + } else { + dest[off_words] |= y << off_bits; + dest[off_words + 1] = y >> (64 - off_bits); + } + off += count; + } else { + // 32 bit + HsInt off_words = off >> 5; + HsInt off_bits = off & 0x1f; + if (off_bits == 0) { + dest[off_words] = y; + } else { + dest[off_words] |= y << off_bits; + dest[off_words + 1] = y >> (32 - off_bits); + } + off += count; + } + } + return off; +} diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,15 @@ +# 1.1.5.0 + +* Make `zipBits` unconditionally strict in its second bit + vector argument (thanks to @treeowl). + +* Add `simd` flag (enabled by default) to use a C SIMD + implementation for `zipBits`, `invertBits`, `countBits`, + `bitIndex`, `nthBitIndex`, `selectBits`, `excludeBits`, + `reverseBits` (thanks to @konsumlamm). + +* Decomission `libgmp` flag. + # 1.1.4.0 * Include `Data.Bit.Gmp` only if `libgmp` flag is set. diff --git a/src/Data/Bit.hs b/src/Data/Bit.hs --- a/src/Data/Bit.hs +++ b/src/Data/Bit.hs @@ -8,7 +8,9 @@ -- Licence: BSD3 -- Maintainer: Andrew Lelechenko -- --- This module exposes an interface with thread-unsafe writes and flips. +-- This module exposes an interface with non-thread-safe writes and flips. +-- Additionally, concurrently modifying non-intersecting slices of the same underlying array +-- may lead to unexpected results. -- Consider using "Data.Bit.ThreadSafe", which is thread-safe, but slower -- (usually 10-20%, up to 50% for short vectors). -- @@ -22,8 +24,11 @@ -- Maintainer: Andrew Lelechenko -- -- This module exposes an interface with thread-safe writes and flips. +-- Additionally, concurrently modifying non-intersecting slices of the same underlying array +-- works as expected. However, operations that affect multiple elements are not +-- guaranteed to be atomic. -- Consider using "Data.Bit", which is faster --- (usually 10-20%, up to 50% for short vectors), but thread-unsafe. +-- (usually 10-20%, up to 50% for short vectors), but not thread-safe. -- -- @since 1.0 module Data.Bit.ThreadSafe diff --git a/src/Data/Bit/Gmp.hs b/src/Data/Bit/Gmp.hs deleted file mode 100644 --- a/src/Data/Bit/Gmp.hs +++ /dev/null @@ -1,101 +0,0 @@ -{-# LANGUAGE MagicHash #-} -{-# LANGUAGE UnliftedFFITypes #-} - -module Data.Bit.Gmp - ( mpnCom - , mpnPopcount - , mpnAndN - , mpnIorN - , mpnXorN - , mpnAndnN - , mpnIornN - , mpnNandN - , mpnNiorN - , mpnXnorN - ) where - -import Control.Monad.ST -import Control.Monad.ST.Unsafe -import Data.Primitive.ByteArray -import GHC.Exts -import System.IO.Unsafe - -foreign import ccall unsafe "__gmpn_com" - mpn_com :: MutableByteArray# s -> ByteArray# -> Int# -> IO () - -mpnCom :: MutableByteArray s -> ByteArray -> Int -> ST s () -mpnCom (MutableByteArray res#) (ByteArray arg#) (I# limbs#) = - unsafeIOToST (mpn_com res# arg# limbs#) -{-# INLINE mpnCom #-} - -foreign import ccall unsafe "__gmpn_popcount" - mpn_popcount :: ByteArray# -> Int# -> IO Word - -mpnPopcount :: ByteArray -> Int -> Word -mpnPopcount (ByteArray arg#) (I# limbs#) = - unsafeDupablePerformIO (mpn_popcount arg# limbs#) -{-# INLINE mpnPopcount #-} - -foreign import ccall unsafe "__gmpn_and_n" - mpn_and_n :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () - -mpnAndN :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () -mpnAndN (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# limbs#) = - unsafeIOToST (mpn_and_n res# arg1# arg2# limbs#) -{-# INLINE mpnAndN #-} - -foreign import ccall unsafe "__gmpn_ior_n" - mpn_ior_n :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () - -mpnIorN :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () -mpnIorN (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# limbs#) = - unsafeIOToST (mpn_ior_n res# arg1# arg2# limbs#) -{-# INLINE mpnIorN #-} - -foreign import ccall unsafe "__gmpn_xor_n" - mpn_xor_n :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () - -mpnXorN :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () -mpnXorN (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# limbs#) = - unsafeIOToST (mpn_xor_n res# arg1# arg2# limbs#) -{-# INLINE mpnXorN #-} - -foreign import ccall unsafe "__gmpn_andn_n" - mpn_andn_n :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () - -mpnAndnN :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () -mpnAndnN (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# limbs#) = - unsafeIOToST (mpn_andn_n res# arg1# arg2# limbs#) -{-# INLINE mpnAndnN #-} - -foreign import ccall unsafe "__gmpn_iorn_n" - mpn_iorn_n :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () - -mpnIornN :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () -mpnIornN (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# limbs#) = - unsafeIOToST (mpn_iorn_n res# arg1# arg2# limbs#) -{-# INLINE mpnIornN #-} - -foreign import ccall unsafe "__gmpn_nand_n" - mpn_nand_n :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () - -mpnNandN :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () -mpnNandN (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# limbs#) = - unsafeIOToST (mpn_nand_n res# arg1# arg2# limbs#) -{-# INLINE mpnNandN #-} - -foreign import ccall unsafe "__gmpn_nior_n" - mpn_nior_n :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () - -mpnNiorN :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () -mpnNiorN (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# limbs#) = - unsafeIOToST (mpn_nior_n res# arg1# arg2# limbs#) -{-# INLINE mpnNiorN #-} - -foreign import ccall unsafe "__gmpn_xnor_n" - mpn_xnor_n :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () - -mpnXnorN :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () -mpnXnorN (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# limbs#) = - unsafeIOToST (mpn_xnor_n res# arg1# arg2# limbs#) -{-# INLINE mpnXnorN #-} diff --git a/src/Data/Bit/Immutable.hs b/src/Data/Bit/Immutable.hs --- a/src/Data/Bit/Immutable.hs +++ b/src/Data/Bit/Immutable.hs @@ -1,6 +1,7 @@ {-# LANGUAGE CPP #-} {-# LANGUAGE BangPatterns #-} +{-# LANGUAGE BinaryLiterals #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE MagicHash #-} {-# LANGUAGE RankNTypes #-} @@ -44,8 +45,8 @@ import Control.Monad import Control.Monad.ST import Data.Bits -#if UseLibGmp -import Data.Bit.Gmp +#if UseSIMD +import Data.Bit.SIMD #endif #ifndef BITVEC_THREADSAFE import Data.Bit.Internal @@ -69,14 +70,9 @@ import GHC.Exts #endif -#if UseLibGmp -gmpLimbShift :: Int -gmpLimbShift = case wordSize of - 32 -> 2 - 64 -> 3 - _ -> error "gmpLimbShift: unknown architecture" -#endif - +-- | Note: For '(.&.)', '(.|.)' and 'xor', +-- if one input is larger than the other, the remaining bits will be ignored. +-- 'bitSize' is undefined (throws an exception). instance {-# OVERLAPPING #-} Bits (Vector Bit) where (.&.) = zipBits (.&.) (.|.) = zipBits (.|.) @@ -296,15 +292,20 @@ -- | Zip two vectors with the given function. -- Similar to 'Data.Vector.Unboxed.zipWith', --- but up to 1000x (!) faster. +-- but up to 3500x (!) faster. -- +-- Note: If one input is larger than the other, the remaining bits will be ignored. +-- -- For sufficiently dense sets, represented as bitmaps, --- 'zipBits' is up to 32x faster than +-- 'zipBits' is up to 64x faster than -- 'Data.IntSet.union', 'Data.IntSet.intersection', etc. -- --- Users are strongly encouraged to enable the --- @libgmp@ flag for the ultimate performance of 'zipBits'. +-- The function passed to zipBits may only use the following +-- 'Bits' methods: -- +-- '.&.', '.|.', 'xor', 'complement', 'zeroBits', and (likely uselessly) +-- 'bitSizeMaybe' and 'isSigned'. +-- -- >>> :set -XOverloadedLists -- >>> import Data.Bits -- >>> zipBits (.&.) [1,1,0] [0,1,1] -- intersection @@ -322,45 +323,98 @@ -> U.Vector Bit -> U.Vector Bit -> U.Vector Bit -zipBits _ (BitVec _ 0 _) _ = U.empty -zipBits _ _ (BitVec _ 0 _) = U.empty -#if UseLibGmp -zipBits f (BitVec 0 l1 arg1) (BitVec 0 l2 arg2) = runST $ do - let l = l1 `min` l2 +zipBits f = \xs ys -> case (xs, ys) of + (BitVec _ 0 _, !_) -> U.empty + (_, BitVec _ 0 _) -> U.empty +#if UseSIMD + (BitVec 0 l1 arg1, BitVec 0 l2 arg2) -> runST $ do + let + l = noinlineMin l1 l2 w = nWords l - b = w `shiftL` gmpLimbShift + b = wordsToBytes w brr <- newByteArray b - let ff = unBit $ f (Bit False) (Bit False) - ft = unBit $ f (Bit False) (Bit True) - tf = unBit $ f (Bit True) (Bit False) - tt = unBit $ f (Bit True) (Bit True) - case (ff, ft, tf, tt) of - (False, False, False, False) -> setByteArray brr 0 w (zeroBits :: Word) - (False, False, False, True) -> mpnAndN brr arg1 arg2 w - (False, False, True, False) -> mpnAndnN brr arg1 arg2 w - (False, False, True, True) -> copyByteArray brr 0 arg1 0 b - (False, True, False, False) -> mpnAndnN brr arg2 arg1 w - (False, True, False, True) -> copyByteArray brr 0 arg2 0 b - (False, True, True, False) -> mpnXorN brr arg1 arg2 w - (False, True, True, True) -> mpnIorN brr arg1 arg2 w - (True, False, False, False) -> mpnNiorN brr arg1 arg2 w - (True, False, False, True) -> mpnXnorN brr arg1 arg2 w - (True, False, True, False) -> mpnCom brr arg2 w - (True, False, True, True) -> mpnIornN brr arg1 arg2 w - (True, True, False, False) -> mpnCom brr arg1 w - (True, True, False, True) -> mpnIornN brr arg2 arg1 w - (True, True, True, False) -> mpnNandN brr arg1 arg2 w - (True, True, True, True) -> setByteArray brr 0 w (complement zeroBits :: Word) + -- We used to calculate (f False False, f False True, f True False, f True True). + -- Now we calculate all those in one go by passing all four possibilities within + -- a word. + case 0b1111 .&. (unBitsy $ f (Bitsy 0b0011) (Bitsy 0b0101)) of + 0b0000 -> setByteArray brr 0 w (zeroBits :: Word) + 0b0001 -> ompAnd brr arg1 arg2 b + 0b0010 -> ompAndn brr arg1 arg2 b + 0b0011 -> copyByteArray brr 0 arg1 0 b + 0b0100 -> ompAndn brr arg2 arg1 b + 0b0101 -> copyByteArray brr 0 arg2 0 b + 0b0110 -> ompXor brr arg1 arg2 b + 0b0111 -> ompIor brr arg1 arg2 b + 0b1000 -> ompNior brr arg1 arg2 b + 0b1001 -> ompXnor brr arg1 arg2 b + 0b1010 -> ompCom brr arg2 b + 0b1011 -> ompIorn brr arg1 arg2 b + 0b1100 -> ompCom brr arg1 b + 0b1101 -> ompIorn brr arg2 arg1 b + 0b1110 -> ompNand brr arg1 arg2 b + _0b1111 -> setByteArray brr 0 w (complement zeroBits :: Word) BitVec 0 l <$> unsafeFreezeByteArray brr #endif -zipBits f xs ys = runST $ do - let n = min (U.length xs) (U.length ys) - zs <- MU.new n - forM_ [0, wordSize .. n - 1] $ \i -> - writeWord zs i (f (indexWord xs i) (indexWord ys i)) - U.unsafeFreeze zs -{-# INLINABLE zipBits #-} + _ -> runST $ do + let n = noinlineMin (U.length xs) (U.length ys) + zs <- MU.new n + forM_ [0, wordSize .. n - 1] $ \i -> + writeWord zs i . unBitsy $ f (Bitsy $ indexWord xs i) (Bitsy $ indexWord ys i) + U.unsafeFreeze zs +{-# INLINE zipBits #-} +-- | This is hideous, but it keeps the code size down in applications of +-- 'zipBits'. Otherwise we end up taking different code paths depending +-- on how the comparison goes in the min calculation, and the Core gets +-- seriously ugly. Ugh! +noinlineMin :: Int -> Int -> Int +noinlineMin = min +{-# NOINLINE noinlineMin #-} + +-- | A version of 'Word' that only supports operations that make sense in +-- zipBits. This ensures that if someone does something overly silly in the function +-- they pass to zipBits, then they'll get a helpful (albeit run-time) error rather than just +-- weird garbage results. +newtype Bitsy = Bitsy {unBitsy :: Word} +instance Eq Bitsy where + _ == _ = notBitsy "==" +instance Bits Bitsy where + Bitsy x .&. Bitsy y = Bitsy (x .&. y) + Bitsy x .|. Bitsy y = Bitsy (x .|. y) + Bitsy x `xor` Bitsy y = Bitsy (x `xor` y) + complement (Bitsy x) = Bitsy (complement x) + zeroBits = Bitsy zeroBits + bitSizeMaybe _ = Nothing + isSigned _ = False -- Not useful, but not harmful + {-# INLINE (.&.) #-} + {-# INLINE (.|.) #-} + {-# INLINE xor #-} + {-# INLINE complement #-} + {-# INLINE zeroBits #-} + + shiftL _ _ = notBitsy "shiftL" + shiftR _ _ = notBitsy "shiftR" + shift _ _ = notBitsy "shift" + unsafeShiftL _ _ = notBitsy "unsafeShiftL" + unsafeShiftR _ _ = notBitsy "unsafeShiftR" + rotateL _ _ = notBitsy "rotateL" + rotateR _ _ = notBitsy "rotateR" + rotate _ _ = notBitsy "rotate" + bitSize _ = notBitsy "bitSize" + testBit _ _ = notBitsy "testBit" + bit _ = notBitsy "bit" + setBit _ _ = notBitsy "setBit" + clearBit _ _ = notBitsy "clearBit" + complementBit _ _ = notBitsy "complementBit" + popCount _ = notBitsy "popCount" + +{-# NOINLINE notBitsy #-} +notBitsy :: String -> a +notBitsy fun = error $ + "The function passed to zipBits may only use\n" ++ + ".&., .|., xor, complement, zeroBits, bitSizeMaybe, and isSigned.\n" ++ + "You used " ++ fun + -- | Map a vectors with the given function. -- Similar to 'Data.Vector.Unboxed.map', -- but faster. @@ -384,9 +438,6 @@ -- | Invert (flip) all bits. -- --- Users are strongly encouraged to enable the --- @libgmp@ flag for the ultimate performance of 'invertBits'. --- -- >>> :set -XOverloadedLists -- >>> invertBits [0,1,0,1,0] -- [1,0,1,0,1] @@ -396,11 +447,12 @@ :: U.Vector Bit -> U.Vector Bit invertBits (BitVec _ 0 _) = U.empty -#if UseLibGmp +#if UseSIMD invertBits (BitVec 0 l arg) = runST $ do let w = nWords l - brr <- newByteArray (w `shiftL` gmpLimbShift) - mpnCom brr arg w + b = wordsToBytes w + brr <- newByteArray b + ompCom brr arg b BitVec 0 l <$> unsafeFreezeByteArray brr #endif invertBits xs = runST $ do @@ -410,11 +462,13 @@ writeWord ys i (complement (indexWord xs i)) U.unsafeFreeze ys --- | For each set bit of the first argument, deposit +-- | For each set bit of the first argument, extract -- the corresponding bit of the second argument -- to the result. Similar to the --- [parallel bit deposit instruction (PDEP)](https://en.wikipedia.org/wiki/X86_Bit_manipulation_instruction_set#Parallel_bit_deposit_and_extract). +-- [parallel bit extract instruction (PEXT)](https://en.wikipedia.org/wiki/X86_Bit_manipulation_instruction_set#Parallel_bit_deposit_and_extract). -- +-- Note: If one input is larger than the other, the remaining bits will be ignored. +-- -- >>> :set -XOverloadedLists -- >>> selectBits [0,1,0,1,1] [1,1,0,0,1] -- [1,0,1] @@ -426,15 +480,25 @@ -- -- @since 0.1 selectBits :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit +#ifdef UseSIMD +selectBits (BitVec 0 iLen iArr) (BitVec 0 xLen xArr) | modWordSize len == 0 = runST $ do + marr <- newByteArray (len `shiftR` 3) + n <- selectBitsC marr xArr iArr (divWordSize len) False + BitVec 0 n <$> unsafeFreezeByteArray marr + where + len = min iLen xLen +#endif selectBits is xs = runST $ do xs1 <- U.thaw xs n <- selectBitsInPlace is xs1 U.unsafeFreeze (MU.take n xs1) --- | For each unset bit of the first argument, deposit +-- | For each unset bit of the first argument, extract -- the corresponding bit of the second argument -- to the result. -- +-- Note: If one input is larger than the other, the remaining bits will be ignored. +-- -- >>> :set -XOverloadedLists -- >>> excludeBits [0,1,0,1,1] [1,1,0,0,1] -- [1,0] @@ -446,6 +510,14 @@ -- -- @since 0.1 excludeBits :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit +#ifdef UseSIMD +excludeBits (BitVec 0 iLen iArr) (BitVec 0 xLen xArr) | modWordSize len == 0 = runST $ do + marr <- newByteArray (len `shiftR` 3) + n <- selectBitsC marr xArr iArr (divWordSize len) True + BitVec 0 n <$> unsafeFreezeByteArray marr + where + len = min iLen xLen +#endif excludeBits is xs = runST $ do xs1 <- U.thaw xs n <- excludeBitsInPlace is xs1 @@ -462,6 +534,12 @@ -- -- @since 1.0.1.0 reverseBits :: U.Vector Bit -> U.Vector Bit +#ifdef UseSIMD +reverseBits (BitVec 0 len arr) | modWordSize len == 0 = runST $ do + marr <- newByteArray (len `shiftR` 3) + reverseBitsC marr arr (divWordSize len) + BitVec 0 len <$> unsafeFreezeByteArray marr +#endif reverseBits xs = runST $ do let n = U.length xs ys <- MU.new n @@ -505,6 +583,11 @@ -- -- @since 1.0.0.0 bitIndex :: Bit -> U.Vector Bit -> Maybe Int +#if UseSIMD +bitIndex (Bit b) (BitVec 0 len arr) | modWordSize len == 0 = + let res = bitIndexC arr (divWordSize len) b + in if res < 0 then Nothing else Just res +#endif bitIndex b (BitVec off len arr) | len == 0 = Nothing | offBits == 0 = case modWordSize len of @@ -589,6 +672,11 @@ -- @since 1.0.0.0 nthBitIndex :: Bit -> Int -> U.Vector Bit -> Maybe Int nthBitIndex _ k _ | k <= 0 = error "nthBitIndex: n must be positive" +#if UseSIMD +nthBitIndex (Bit b) n (BitVec 0 len arr) | modWordSize len == 0 = + let res = nthBitIndexC arr (divWordSize len) b n + in if res < 0 then Nothing else Just res +#endif nthBitIndex b k (BitVec off len arr) | len == 0 = Nothing | offBits == 0 = either (const Nothing) Just $ case modWordSize len of @@ -666,9 +754,6 @@ -- | Return the number of set bits in a vector (population count, popcount). -- --- Users are strongly encouraged to enable the --- @libgmp@ flag for the ultimate performance of 'countBits'. --- -- >>> :set -XOverloadedLists -- >>> countBits [1,1,0,1,0,1] -- 4 @@ -680,9 +765,9 @@ -- @since 0.1 countBits :: U.Vector Bit -> Int countBits (BitVec _ 0 _) = 0 -#if UseLibGmp +#if UseSIMD countBits (BitVec 0 len arr) | modWordSize len == 0 = - fromIntegral (mpnPopcount arr (divWordSize len)) + ompPopcount arr (len `shiftR` 5) #endif countBits (BitVec off len arr) | offBits == 0 = case modWordSize len of 0 -> countBitsInWords (P.Vector offWords lWords arr) diff --git a/src/Data/Bit/Internal.hs b/src/Data/Bit/Internal.hs --- a/src/Data/Bit/Internal.hs +++ b/src/Data/Bit/Internal.hs @@ -216,7 +216,7 @@ -- | Write a word at the given bit offset in little-endian order (i.e., the LSB will correspond to the bit at the given address, the 2's bit will correspond to the address + 1, etc.). If the offset is such that the word extends past the end of the vector, the word is truncated and as many low-order bits as possible are written. writeWord :: PrimMonad m => U.MVector (PrimState m) Bit -> Int -> Word -> m () -writeWord (BitMVec _ 0 _) _ _ = pure () +writeWord (BitMVec _ 0 _) !_ !_ = pure () writeWord (BitMVec off len' arr) !i' !x | iMod == 0 = if len >= i + wordSize diff --git a/src/Data/Bit/Mutable.hs b/src/Data/Bit/Mutable.hs --- a/src/Data/Bit/Mutable.hs +++ b/src/Data/Bit/Mutable.hs @@ -116,6 +116,8 @@ -- rewriting the contents of the second argument. -- Cf. 'Data.Bit.zipBits'. -- +-- Note: If one input is larger than the other, the remaining bits will be ignored. +-- -- >>> :set -XOverloadedLists -- >>> import Data.Bits -- >>> Data.Vector.Unboxed.modify (zipInPlace (.&.) [1,1,0]) [0,1,1] @@ -240,10 +242,12 @@ writeWord xs i (complement x) {-# SPECIALIZE invertInPlace :: U.MVector s Bit -> ST s () #-} --- | Same as 'Data.Bit.selectBits', but deposit +-- | Same as 'Data.Bit.selectBits', but extract -- selected bits in-place. Returns the number of selected bits. -- It is the caller's responsibility to trim the result to this number. -- +-- Note: If one input is larger than the other, the remaining bits will be ignored. +-- -- >>> :set -XOverloadedLists -- >>> import Control.Monad.ST (runST) -- >>> import qualified Data.Vector.Unboxed as U @@ -263,11 +267,14 @@ let !(nSet, x') = selectWord (masked (n - i) (indexWord is i)) x writeWord xs ct x' loop (i + wordSize) (ct + nSet) +{-# SPECIALIZE selectBitsInPlace :: U.Vector Bit -> U.MVector s Bit -> ST s Int #-} --- | Same as 'Data.Bit.excludeBits', but deposit +-- | Same as 'Data.Bit.excludeBits', but extract -- excluded bits in-place. Returns the number of excluded bits. -- It is the caller's responsibility to trim the result to this number. -- +-- Note: If one input is larger than the other, the remaining bits will be ignored. +-- -- >>> :set -XOverloadedLists -- >>> import Control.Monad.ST (runST) -- >>> import qualified Data.Vector.Unboxed as U @@ -288,6 +295,7 @@ selectWord (masked (n - i) (complement (indexWord is i))) x writeWord xs ct x' loop (i + wordSize) (ct + nSet) +{-# SPECIALIZE excludeBitsInPlace :: U.Vector Bit -> U.MVector s Bit -> ST s Int #-} -- | Reverse the order of bits in-place. -- diff --git a/src/Data/Bit/SIMD.hs b/src/Data/Bit/SIMD.hs new file mode 100644 --- /dev/null +++ b/src/Data/Bit/SIMD.hs @@ -0,0 +1,156 @@ +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE UnliftedFFITypes #-} + +module Data.Bit.SIMD + ( ompPopcount + , ompCom + , ompAnd + , ompIor + , ompXor + , ompAndn + , ompIorn + , ompNand + , ompNior + , ompXnor + , reverseBitsC + , bitIndexC + , nthBitIndexC + , selectBitsC + ) where + +import Control.Monad.ST +import Control.Monad.ST.Unsafe (unsafeIOToST) +import Data.Primitive.ByteArray +import GHC.Exts + +foreign import ccall unsafe "_hs_bitvec_popcount" + omp_popcount :: ByteArray# -> Int# -> Int# + +-- | SIMD optimized popcount. The length is in 32 bit words. +ompPopcount :: ByteArray -> Int -> Int +ompPopcount (ByteArray arg#) (I# len#) = + I# (omp_popcount arg# len#) +{-# INLINE ompPopcount #-} + +foreign import ccall unsafe "_hs_bitvec_com" + omp_com :: MutableByteArray# s -> ByteArray# -> Int# -> IO () + +-- | SIMD optimized bitwise complement. The length is in bytes +-- and the result array should have at least that many bytes. +ompCom :: MutableByteArray s -> ByteArray -> Int -> ST s () +ompCom (MutableByteArray res#) (ByteArray arg#) (I# len#) = + unsafeIOToST (omp_com res# arg# len#) +{-# INLINE ompCom #-} + +foreign import ccall unsafe "_hs_bitvec_and" + omp_and :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () + +-- | SIMD optimized bitwise AND. The length is in bytes +-- and the result array should have at least that many bytes. +ompAnd :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () +ompAnd (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# len#) = + unsafeIOToST (omp_and res# arg1# arg2# len#) +{-# INLINE ompAnd #-} + +foreign import ccall unsafe "_hs_bitvec_ior" + omp_ior :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () + +-- | SIMD optimized bitwise OR. The length is in bytes +-- and the result array should have at least that many bytes. +ompIor :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () +ompIor (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# len#) = + unsafeIOToST (omp_ior res# arg1# arg2# len#) +{-# INLINE ompIor #-} + +foreign import ccall unsafe "_hs_bitvec_xor" + omp_xor :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () + +-- | SIMD optimized bitwise XOR. The length is in bytes +-- and the result array should have at least that many bytes. +ompXor :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () +ompXor (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# len#) = + unsafeIOToST (omp_xor res# arg1# arg2# len#) +{-# INLINE ompXor #-} + +foreign import ccall unsafe "_hs_bitvec_andn" + omp_andn :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () + +-- | SIMD optimized bitwise AND with the second argument inverted. The length is in bytes +-- and the result array should have at least that many bytes. +ompAndn :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () +ompAndn (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# len#) = + unsafeIOToST (omp_andn res# arg1# arg2# len#) +{-# INLINE ompAndn #-} + +foreign import ccall unsafe "_hs_bitvec_iorn" + omp_iorn :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () + +-- | SIMD optimized bitwise OR with the second argument inverted. The length is in bytes +-- and the result array should have at least that many bytes. +ompIorn :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () +ompIorn (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# len#) = + unsafeIOToST (omp_iorn res# arg1# arg2# len#) +{-# INLINE ompIorn #-} + +foreign import ccall unsafe "_hs_bitvec_nand" + omp_nand :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () + +-- | SIMD optimized bitwise NAND. The length is in bytes +-- and the result array should have at least that many bytes. +ompNand :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () +ompNand (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# len#) = + unsafeIOToST (omp_nand res# arg1# arg2# len#) +{-# INLINE ompNand #-} + +foreign import ccall unsafe "_hs_bitvec_nior" + omp_nior :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () + +-- | SIMD optimized bitwise NOR. The length is in bytes +-- and the result array should have at least that many bytes. +ompNior :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () +ompNior (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# len#) = + unsafeIOToST (omp_nior res# arg1# arg2# len#) +{-# INLINE ompNior #-} + +foreign import ccall unsafe "_hs_bitvec_xnor" + omp_xnor :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> IO () + +-- | SIMD optimized bitwise XNOR. The length is in bytes +-- and the result array should have at least that many bytes. +ompXnor :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> ST s () +ompXnor (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# len#) = + unsafeIOToST (omp_xnor res# arg1# arg2# len#) +{-# INLINE ompXnor #-} + +foreign import ccall unsafe "_hs_bitvec_reverse_bits" + reverse_bits :: MutableByteArray# s -> ByteArray# -> Int# -> IO () + +-- | The length is in words. +reverseBitsC :: MutableByteArray s -> ByteArray -> Int -> ST s () +reverseBitsC (MutableByteArray res#) (ByteArray arg#) (I# len#) = + unsafeIOToST (reverse_bits res# arg# len#) +{-# INLINE reverseBitsC #-} + +foreign import ccall unsafe "_hs_bitvec_bit_index" + bit_index :: ByteArray# -> Int# -> Bool -> Int# + +bitIndexC :: ByteArray -> Int -> Bool -> Int +bitIndexC (ByteArray arg#) (I# len#) bit = + I# (bit_index arg# len# bit) +{-# INLINE bitIndexC #-} + +foreign import ccall unsafe "_hs_bitvec_nth_bit_index" + nth_bit_index :: ByteArray# -> Int# -> Bool -> Int# -> Int# + +nthBitIndexC :: ByteArray -> Int -> Bool -> Int -> Int +nthBitIndexC (ByteArray arg#) (I# len#) bit (I# n#) = + I# (nth_bit_index arg# len# bit n#) +{-# INLINE nthBitIndexC #-} + +foreign import ccall unsafe "_hs_bitvec_select_bits" + select_bits_c :: MutableByteArray# s -> ByteArray# -> ByteArray# -> Int# -> Bool -> IO Int + +selectBitsC :: MutableByteArray s -> ByteArray -> ByteArray -> Int -> Bool -> ST s Int +selectBitsC (MutableByteArray res#) (ByteArray arg1#) (ByteArray arg2#) (I# len#) exclude = + unsafeIOToST (select_bits_c res# arg1# arg2# len# exclude) +{-# INLINE selectBitsC #-} diff --git a/test/Tests/SetOps.hs b/test/Tests/SetOps.hs --- a/test/Tests/SetOps.hs +++ b/test/Tests/SetOps.hs @@ -48,8 +48,8 @@ , testProperty "reverseInPlace middle" prop_reverseInPlace_middle , testProperty "reverseInPlaceLong middle" prop_reverseInPlaceLong_middle - , testProperty "selectBits" prop_selectBits_def - , testProperty "excludeBits" prop_excludeBits_def + , mkGroup2 "selectBits" prop_selectBits_def + , mkGroup2 "excludeBits" prop_excludeBits_def , mkGroup "countBits" prop_countBits_def ] @@ -68,6 +68,21 @@ prop (U.slice from len (U.generate (from + len + excess) (Bit . f))) propMiddleLong (NonNegative x) (NonNegative y) (NonNegative z) = propMiddle (NonNegative $ x * 31) (NonNegative $ y * 37) (NonNegative $ z * 29) + +mkGroup2 :: String -> (U.Vector Bit -> U.Vector Bit -> Property) -> TestTree +mkGroup2 name prop = testGroup name + [ testProperty "simple" prop + , testProperty "simple_long" (\(Large xs) (Large ys) -> prop xs ys) + , testProperty "middle" propMiddle + , testProperty "middle_long" propMiddleLong + ] + where + f m = let n = fromIntegral m :: Double in + odd (truncate (exp (abs (sin n) * 10)) :: Integer) + propMiddle (NonNegative from1) (NonNegative len1) (NonNegative excess1) (NonNegative from2) (NonNegative len2) (NonNegative excess2) = + prop (U.slice from1 len1 (U.generate (from1 + len1 + excess1) (Bit . f))) (U.slice from2 len2 (U.generate (from2 + len2 + excess2) (Bit . f))) + propMiddleLong (NonNegative x1) (NonNegative y1) (NonNegative z1) (NonNegative x2) (NonNegative y2) (NonNegative z2) = + propMiddle (NonNegative $ x1 * 31) (NonNegative $ y1 * 37) (NonNegative $ z1 * 29) (NonNegative $ x2 * 31) (NonNegative $ y2 * 37) (NonNegative $ z2 * 29) prop_generalize1 :: Fun Bit Bit -> Bit -> Property prop_generalize1 fun x =