diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2019 Andrew Lelechenko, 2012-2016 James Cook
+Copyright (c) 2019-2022 Andrew Lelechenko, 2012-2016 James Cook
 
 All rights reserved.
 
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -1,27 +1,32 @@
-# bitvec [![Build Status](https://travis-ci.org/Bodigrim/bitvec.svg)](https://travis-ci.org/Bodigrim/bitvec) [![Hackage](http://img.shields.io/hackage/v/bitvec.svg)](https://hackage.haskell.org/package/bitvec) [![Hackage CI](https://matrix.hackage.haskell.org/api/v2/packages/bitvec/badge)](https://matrix.hackage.haskell.org/package/bitvec) [![Stackage LTS](http://stackage.org/package/bitvec/badge/lts)](http://stackage.org/lts/package/bitvec) [![Stackage Nightly](http://stackage.org/package/bitvec/badge/nightly)](http://stackage.org/nightly/package/bitvec)
+# 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`
+package represents unboxed arrays of `Bool`s
 spending 1 byte (8 bits) per boolean.
 This library provides a newtype wrapper `Bit` and a custom instance
-of unboxed `Vector`, which packs bits densely,
-achieving __8x less memory footprint.__
+of an unboxed `Vector`, which packs bits densely,
+achieving an __8x smaller memory footprint.__
 The performance stays mostly the same;
 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:
-  read the whole word from memory,
+  they read the whole word from memory,
   then modify a bit, then write the whole word back.
-* `Data.Bit.ThreadSafe` is slower (up to 20%),
+  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
 
@@ -73,7 +78,7 @@
   U.unsafeFreeze sieve
 ```
 
-`Bit`-based implementation requires 8x less memory to store
+The `Bit`-based implementation requires 8x less memory to store
 the vector. For large sizes it allows to crunch more data in RAM
 without swapping. For smaller arrays it helps to fit into
 CPU caches.
@@ -84,8 +89,8 @@
 ```
 
 There are several high-level helpers, digesting bits in bulk,
-which makes them up to 64x faster than respective counterparts
-for `Vector Bool`. One can query population count (popcount)
+which makes them up to 64x faster than the respective counterparts
+for `Vector Bool`. One can query the population count (popcount)
 of a vector (giving us [the prime-counting function](https://en.wikipedia.org/wiki/Prime-counting_function)):
 
 ```haskell
@@ -104,13 +109,13 @@
 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.
-There is a moderate performance penalty (up to 20%)
+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.
 
 ## Sets
 
-Bit vectors can be used as a blazingly fast representation of sets
+Bit vectors can be used as a blazingly fast representation of sets,
 as long as their elements are `Enum`eratable and sufficiently dense,
 leaving `IntSet` far behind.
 
@@ -120,22 +125,64 @@
 * As a 64k-long unboxed `Vector Bool`, implementing union as `zipWith (||)`.
 * As a 64k-long unboxed `Vector Bit`, implementing union as `zipBits (.|.)`.
 
-In our benchmarks (see `bench` folder) for not-too-sparse sets
-the union of `Vector Bit` evaluates 24x-36x faster than the union of `IntSet`
-and stunningly outperforms `Vector Bool` 500x-1000x.
+When the `simd` flag is enabled,
+according to our benchmarks (see `bench` folder),
+the union of `Vector Bit` evaluates magnitudes faster
+than the union of not-too-sparse `IntSet`s
+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.
 Their applications include coding theory and cryptography.
-While one can successfully implement them with `poly` package,
+While one can successfully implement them with the [`poly`](https://hackage.haskell.org/package/poly) package,
 operating on `UPoly Bit`,
 this package provides even faster arithmetic routines
-exposed via `F2Poly` data type and its instances.
+exposed via the `F2Poly` data type and its instances.
 
 ```haskell
 > :set -XBinaryLiterals
-> -- (1 + x) (1 + x + x^2) = 1 + x^3 (mod 2)
+> -- (1 + x) * (1 + x + x^2) = 1 + x^3 (mod 2)
 > 0b11 * 0b111 :: F2Poly
 F2Poly {unF2Poly = [1,0,0,1]}
 ```
@@ -145,21 +192,9 @@
 
 ## Package flags
 
-This package supports the following flags to facilitate dependency management.
-Disabling them does not diminish `bitvec`'s capabilities, but makes certain operations slower.
-
-* Flag `integer-gmp`, enabled by default.
-
-  Depend on `integer-gmp` package and use it to speed up operations on binary polynomials.
-  Normally `integer-gmp` is shipped with core libraries anyways, so there is little to gain
-  from disabling it, unless you use a custom build of GHC.
-
-* Flag `libgmp`, enabled by default.
+* Flag `simd`, enabled by default.
 
-  Link against [GMP](https://gmplib.org/) library and use it to for ultimate performance of
-  `zipBits`, `invertBits` and `countBits`. GMP is readily available on most machines
-  (`brew install gmp` on macOS), but you may find useful to disable this flag working
-  with exotic setup.
+  Use a C SIMD implementation for the ultimate performance of `zipBits`, `invertBits` and `countBits`.
 
 ## Similar packages
 
@@ -170,3 +205,8 @@
 * [`array`](https://hackage.haskell.org/package/array)
   is memory-efficient for `Bool`, but lacks
   a handy `Vector` interface and is not thread-safe.
+
+## Additional resources
+
+* __Bit vectors without compromises__, Haskell Love, 31.07.2020:
+  [slides](https://github.com/Bodigrim/my-talks/raw/master/haskelllove2020/slides.pdf), [video](https://youtu.be/HhpH8DKFBls).
diff --git a/Setup.lhs b/Setup.lhs
deleted file mode 100644
--- a/Setup.lhs
+++ /dev/null
@@ -1,5 +0,0 @@
-#!/usr/bin/env runhaskell
-
-> import Distribution.Simple
-> main = defaultMain
-
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -1,8 +1,10 @@
-module Main where
+module Main (main) where
 
-import Gauge.Main
+import Test.Tasty.Bench
+import Test.Tasty.Patterns.Printer
 
 import Bench.BitIndex
+import Bench.Common
 import Bench.GCD
 import Bench.Intersection
 import Bench.Invert
@@ -16,17 +18,25 @@
 import Bench.Union
 
 main :: IO ()
-main = defaultMain
-  [ bgroup "bitIndex"     $ map benchBitIndex     [5..14]
-  , bgroup "invert"       $ map benchInvert       [5..14]
+main = defaultMain $ map (mapLeafBenchmarks addCompare)
+  [ 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 "randomWrite"  $ map benchRandomWrite  [5..14]
-  , bgroup "randomFlip"   $ map benchRandomFlip   [5..14]
-  , bgroup "randomRead"   $ map benchRandomRead   [5..14]
-  , bgroup "remainder"    $ map benchRemainder    [5..14]
+  , bgroup "productShort" $ map benchProductShort [5..14]
+  , bgroup "read"         $ map benchRandomRead   [5..14]
+  , bgroup "remainder"    $ map benchRemainder    [5..11]
   , bgroup "reverse"      $ map benchReverse      [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]
   ]
+
+addCompare :: ([String] -> Benchmark -> Benchmark)
+addCompare (name : 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
@@ -7,41 +7,18 @@
 import Data.Bits
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
-import Gauge.Main
+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/bitIndex"       $ nf bitIndexBit      (randomVec Bit k)
-  , bench "Bit/nthBitIndex"    $ nf nthBitIndexBit   (randomVec Bit k)
-  , bench "Bit/elemIndex"      $ nf elemIndexBit     (randomVec Bit k)
-  , bench "Bit.TS/bitIndex"    $ nf bitIndexBitTS    (randomVec TS.Bit k)
-  , bench "Bit.TS/nthBitIndex" $ nf nthBitIndexBitTS (randomVec TS.Bit k)
-  , bench "Bit.TS/elemIndex"   $ nf elemIndexBitTS   (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)
-
-nthBitIndexBit :: U.Vector Bit -> Maybe Int
-nthBitIndexBit = nthBitIndex (Bit True) 1
-
-elemIndexBit :: U.Vector Bit -> Maybe Int
-elemIndexBit = U.elemIndex (Bit True)
-
-bitIndexBitTS :: U.Vector TS.Bit -> Maybe Int
-bitIndexBitTS = TS.bitIndex (TS.Bit True)
-
-nthBitIndexBitTS :: U.Vector TS.Bit -> Maybe Int
-nthBitIndexBitTS = TS.nthBitIndex (TS.Bit True) 1
-
-elemIndexBitTS :: U.Vector TS.Bit -> Maybe Int
-elemIndexBitTS = U.elemIndex (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 Gauge.Main
-import System.Random
-
-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 Test.Tasty.Bench
 
-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/gcdExt"    $ nf (uncurry    gcdExt) (   toF2Poly $ randomVec    Bit k,    toF2Poly $ randomVec'    Bit k)
-  , bench "Bit.TS/gcdExt" $ 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,5 @@
+{- HLINT ignore "Avoid lambda" -}
+
 module Bench.Intersection
   ( benchIntersection
   ) where
@@ -7,61 +9,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 Gauge.Main
-import System.Random
-
-randomBools :: [Bool]
-randomBools
-  = map (\i -> if i > (0 :: Int) then True else False)
-  . 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
+import Test.Tasty.Bench
 
-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/zipBits"    $ nf (\x -> intersectionBit    (randomVec Bit k) x)    (randomVec2 Bit k)
-  , bench "Bit/zipWith"    $ nf (\x -> intersectionBit'   (randomVec Bit k) x)    (randomVec2 Bit k)
-  , bench "Bit.TS/zipBits" $ nf (\x -> intersectionBitTS  (randomVec TS.Bit k) x) (randomVec2 TS.Bit k)
-  , bench "Bit.TS/zipWith" $ 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 (.&.)
-
-intersectionBit' :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
-intersectionBit' = U.zipWith (\(Bit x) (Bit y) -> Bit (x && y))
-
-intersectionBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit -> U.Vector TS.Bit
-intersectionBitTS = TS.zipBits (.&.)
-
-intersectionBitTS' :: U.Vector TS.Bit -> U.Vector TS.Bit -> U.Vector TS.Bit
-intersectionBitTS' = U.zipWith (\(TS.Bit x) (TS.Bit y) -> TS.Bit (x && y))
-
-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,42 +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 Gauge.Main
-import System.Random
-
-randomBools :: [Bool]
-randomBools
-  = map (\i -> if i > (0 :: Int) then True else False)
-  . randoms
-  . mkStdGen
-  $ 42
+import Test.Tasty.Bench
 
-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/invertBits"        $ nf invertBit    (randomVec Bit k)
-  , bench "Bit/map-complement"    $ nf invertBit'   (randomVec Bit k)
-  , bench "Bit.TS/invertBits"     $ nf invertBitTS  (randomVec TS.Bit k)
-  , bench "Bit.TS/map-complement" $ 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
-
-invertBit' :: U.Vector Bit -> U.Vector Bit
-invertBit' = U.map complement
-
-invertBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit
-invertBitTS = TS.invertBits
-
-invertBitTS' :: U.Vector TS.Bit -> U.Vector TS.Bit
-invertBitTS' = U.map complement
-
-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,49 +1,38 @@
+{- HLINT ignore "Avoid lambda" -}
+
 module Bench.Product
   ( benchProduct
+  , benchProductShort
+  , benchSquare
   ) where
 
 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 Gauge.Main
-import System.Random
-
-randomBools :: [Bool]
-randomBools
-  = map (\i -> if i > (0 :: Int) then True else False)
-  . 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
+import Test.Tasty.Bench
 
-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/product"          $ nf (\x -> (*) (toF2Poly $ randomVec Bit k) x)    (toF2Poly $ randomVec2 Bit k)
-  , bench "Bit/productShort"     $ nf (\x -> (*) (toF2Poly $ randomVec Bit k) x)    (toF2Poly $ U.take 32 $ randomVec2 Bit k)
-  , bench "Bit/square"           $ nf (\x -> (*) (toF2Poly $ randomVec Bit k) x)    (toF2Poly $ randomVec Bit k)
-  , bench "Bit.TS/product"       $ nf (\x -> (*) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec2 TS.Bit k)
-  , bench "Bit.TS/productShort"  $ nf (\x -> (*) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ U.take 32 $ randomVec2 TS.Bit k)
-  , bench "Bit.TS/square"        $ nf (\x -> (*) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec TS.Bit k)
-  , bench "Integer/product"      $ nf (\x -> binMul (randomInteger k) x) (randomInteger2 k)
-  , bench "Integer/productShort" $ nf (\x -> binMul (randomInteger k) x) ((1 `shiftL` 32 - 1) .&. randomInteger2 k)
-  , bench "Integer/square"       $ nf (\x -> binMul (randomInteger k) x) (randomInteger 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 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 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
@@ -2,75 +2,50 @@
   ( benchRandomFlip
   ) where
 
+import Prelude hiding (Foldable(..))
 import Control.Monad
 import Control.Monad.ST
 import Data.Bit
 import qualified Data.Bit.ThreadSafe as TS
 import Data.Bits
+import Data.Foldable
 import qualified Data.IntSet as IS
-import Data.List
 import qualified Data.Vector.Unboxed.Mutable as MU
-import Gauge.Main
-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/flip"      $ nf randomFlipBit    k
-  , bench "Bit/modify"    $ nf randomFlipBit'   k
-  , bench "Bit.TS/flip"   $ nf randomFlipBitTS  k
-  , bench "Bit.TS/modify" $ 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
 
-randomFlipBit' :: Int -> Int
-randomFlipBit' k = runST $ do
-  let n = 1 `shiftL` k
-  vec <- MU.new n
-  forM_ (take (mult * n) randomFlips) $
-    \i -> MU.unsafeModify vec complement (i .&. (1 `shiftL` k - 1))
-  Bit i <- MU.unsafeRead vec 0
-  pure $ if i then 1 else 0
-
 randomFlipBitTS :: Int -> Int
 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
 
-randomFlipBitTS' :: Int -> Int
-randomFlipBitTS' k = runST $ do
-  let n = 1 `shiftL` k
-  vec <- MU.new n
-  forM_ (take (mult * n) randomFlips) $
-    \i -> MU.unsafeModify vec complement (i .&. (1 `shiftL` k - 1))
-  TS.Bit i <- MU.unsafeRead vec 0
-  pure $ if i then 1 else 0
-
 randomFlipVector :: Int -> Int
 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
@@ -82,7 +57,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 Gauge.Main
-import System.Random
-
-randomVec :: [Bool]
-randomVec
-  = map (\i -> if i > (0 :: Int) then True else False)
-  . randoms
-  . mkStdGen
-  $ 42
+import Test.Tasty.Bench
 
-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 "Bit.TS" $ 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
@@ -2,37 +2,32 @@
   ( benchRandomWrite
   ) where
 
+import Prelude hiding (Foldable(..))
 import Control.Monad
 import Control.Monad.ST
 import Data.Bit
 import qualified Data.Bit.ThreadSafe as TS
 import Data.Bits
+import Data.Foldable
 import qualified Data.IntSet as IS
-import Data.List
 import qualified Data.Vector.Unboxed.Mutable as MU
-import Gauge.Main
-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 "Bit.TS" $ 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 +36,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 +45,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 +57,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,4 +1,5 @@
 {-# LANGUAGE MagicHash #-}
+{- HLINT ignore "Avoid lambda" -}
 
 module Bench.Remainder
   ( benchRemainder
@@ -7,47 +8,23 @@
 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 Gauge.Main
 import GHC.Exts
-import GHC.Integer.Logarithms
-import System.Random
-
-randomBools :: [Bool]
-randomBools
-  = map (\i -> if i > (0 :: Int) then True else False)
-  . 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
+import GHC.Num.Integer
+import Test.Tasty.Bench
 
-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/remainder"     $ nf (\x -> rem (toF2Poly $ randomVec Bit k) x) (toF2Poly $ randomVec2 Bit k)
-  , bench "Bit.TS/remainder"  $ nf (\x -> rem (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec2 TS.Bit k)
-  , bench "Integer/remainder" $ 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
 binRem x y = go x
   where
-    binLog n = I# (integerLog2# n)
+    binLog n = I# (word2Int# (integerLog2# n))
     ly = binLog y
 
     go z = if lz < ly then z else go (z `xor` (y `shiftL` (lz - ly)))
diff --git a/bench/Bench/Reverse.hs b/bench/Bench/Reverse.hs
--- a/bench/Bench/Reverse.hs
+++ b/bench/Bench/Reverse.hs
@@ -6,42 +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 Gauge.Main
-import System.Random
-
-randomBools :: [Bool]
-randomBools
-  = map (\i -> if i > (0 :: Int) then True else False)
-  . randoms
-  . mkStdGen
-  $ 42
+import Test.Tasty.Bench
 
-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/reverseBits"    $ nf reverseBit    (randomVec Bit k)
-  , bench "Bit/reverse"        $ nf reverseBit'   (randomVec Bit k)
-  , bench "Bit.TS/reverseBits" $ nf reverseBitTS  (randomVec TS.Bit k)
-  , bench "Bit.TS/reverse"     $ 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
-
-reverseBit' :: U.Vector Bit -> U.Vector Bit
-reverseBit' = U.reverse
-
-reverseBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit
-reverseBitTS = TS.reverseBits
-
-reverseBitTS' :: U.Vector TS.Bit -> U.Vector TS.Bit
-reverseBitTS' = U.reverse
-
-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,45 +1,29 @@
+{- HLINT ignore "Avoid lambda" -}
+
 module Bench.Sum
-  ( benchSum
+  ( benchAdd
+  , benchSum
   ) where
 
+import Prelude hiding (Foldable(..))
 import Data.Bit
 import qualified Data.Bit.ThreadSafe as TS
 import Data.Bits
-import Data.List
-import qualified Data.Vector.Unboxed as U
-import qualified Data.Vector.Unboxed.Mutable as MU
-import Gauge.Main
-import System.Random
-
-randomBools :: [Bool]
-randomBools
-  = map (\i -> if i > (0 :: Int) then True else False)
-  . 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
+import Data.Foldable
+import Test.Tasty.Bench
 
-randomInteger :: Int -> Integer
-randomInteger k = toInteger $ toF2Poly $ randomVec Bit k
+import Bench.Common
 
-randomInteger2 :: Int -> Integer
-randomInteger2 k = toInteger $ toF2Poly $ randomVec2 Bit k
+benchAdd :: Int -> Benchmark
+benchAdd k = bgroup (show (1 `shiftL` k :: Int))
+  [ 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/add"     $ nf (\x -> (+) (toF2Poly $ randomVec Bit k) x)    (toF2Poly $ randomVec2 Bit k)
-  , bench "Bit/sum"     $ nf (foldl' (+) 0) [(1 :: F2Poly) .. fromInteger (1 `shiftL` k)]
-  , bench "Bit.TS/add"  $ nf (\x -> (+) (TS.toF2Poly $ randomVec TS.Bit k) x) (TS.toF2Poly $ randomVec2 TS.Bit k)
-  , bench "Bit.TS/sum"  $ nf (foldl' (+) 0) [(1 :: TS.F2Poly) .. fromInteger (1 `shiftL` k)]
-  , bench "Integer/add" $ nf (\x -> xor (randomInteger k) x) (randomInteger2 k)
-  , bench "Integer/sum" $ 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,5 @@
+{- HLINT ignore "Avoid lambda" -}
+
 module Bench.Union
   ( benchUnion
   ) where
@@ -7,61 +9,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 Gauge.Main
-import System.Random
-
-randomBools :: [Bool]
-randomBools
-  = map (\i -> if i > (0 :: Int) then True else False)
-  . 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
+import Test.Tasty.Bench
 
-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/zipBits"    $ nf (\x -> unionBit    (randomVec Bit k) x)    (randomVec2 Bit k)
-  , bench "Bit/zipWith"    $ nf (\x -> unionBit'   (randomVec Bit k) x)    (randomVec2 Bit k)
-  , bench "Bit.TS/zipBits" $ nf (\x -> unionBitTS  (randomVec TS.Bit k) x) (randomVec2 TS.Bit k)
-  , bench "Bit.TS/zipWith" $ 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 (.|.)
-
-unionBit' :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
-unionBit' = U.zipWith (\(Bit x) (Bit y) -> Bit (x || y))
-
-unionBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit -> U.Vector TS.Bit
-unionBitTS = TS.zipBits (.|.)
-
-unionBitTS' :: U.Vector TS.Bit -> U.Vector TS.Bit -> U.Vector TS.Bit
-unionBitTS' = U.zipWith (\(TS.Bit x) (TS.Bit y) -> TS.Bit (x || y))
-
-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,36 +1,41 @@
 name: bitvec
-version: 1.0.3.0
-cabal-version: >=1.10
+version: 1.1.6.0
+cabal-version: 2.0
 build-type: Simple
 license: BSD3
 license-file: LICENSE
-copyright: 2019 Andrew Lelechenko, 2012-2016 James Cook
+copyright: 2019-2022 Andrew Lelechenko, 2012-2016 James Cook
 maintainer: Andrew Lelechenko <andrew.lelechenko@gmail.com>
 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 <https://hackage.haskell.org/package/vector vector>
-  package represents unboxed arrays of 'Bool'
+  package represents unboxed arrays of 'Bool's
   spending 1 byte (8 bits) per boolean.
   This library provides a newtype wrapper 'Bit' and a custom instance
-  of unboxed 'Vector', which packs bits densely,
-  achieving __8x less memory footprint.__
+  of an unboxed 'Vector', which packs bits densely,
+  achieving an __8x smaller memory footprint.__
   The performance stays mostly the same;
   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:
-    read the whole word from memory,
+    they read the whole word from memory,
     then modify a bit, then write the whole word back.
-  * "Data.Bit.ThreadSafe" is slower (up to 20%),
+    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
   .
@@ -46,42 +51,38 @@
 author: Andrew Lelechenko <andrew.lelechenko@gmail.com>,
         James Cook <mokus@deepbondi.net>
 
-tested-with: GHC ==8.0.2 GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.1 GHC ==8.8.2 GHC ==8.8.3 GHC ==8.10.1
-extra-source-files:
+tested-with: GHC ==9.0.2 GHC ==9.2.8 GHC ==9.4.8 GHC ==9.6.7 GHC ==9.8.4 GHC ==9.10.3 GHC ==9.12.2 GHC ==9.14.1
+extra-doc-files:
   changelog.md
   README.md
 
 source-repository head
   type: git
-  location: git://github.com/Bodigrim/bitvec.git
-
-flag integer-gmp
-  description: Use integer-gmp package for binary polynomials
-  default: True
+  location: https://github.com/Bodigrim/bitvec.git
 
-flag libgmp
-  description: Link against GMP library
+flag simd
+  description:
+    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
   exposed-modules:
     Data.Bit
     Data.Bit.ThreadSafe
   build-depends:
-    base >=4.9 && <5,
-    deepseq,
-    ghc-prim,
-    primitive >=0.5,
-    vector >=0.11
-  if impl(ghc <8.0)
-    build-depends:
-      semigroups >=0.8
+    base >=4.15 && <5,
+    bytestring >=0.10 && <0.13,
+    deepseq <1.6,
+    ghc-bignum <1.5,
+    primitive >=0.5 && <0.10,
+    vector >=0.13 && <0.14
   default-language: Haskell2010
   hs-source-dirs: src
   other-modules:
     Data.Bit.F2Poly
     Data.Bit.F2PolyTS
-    Data.Bit.Gmp
     Data.Bit.Immutable
     Data.Bit.ImmutableTS
     Data.Bit.Internal
@@ -90,14 +91,14 @@
     Data.Bit.MutableTS
     Data.Bit.PdepPext
     Data.Bit.Utils
-  ghc-options: -O2 -Wall
+  ghc-options: -O2 -Wall -Wcompat
   include-dirs: src
-  if flag(integer-gmp) && impl(ghc >=8.0.1)
-    build-depends: integer-gmp
-    cpp-options: -DUseIntegerGmp
-  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.SIMD
 
 test-suite bitvec-tests
   type: exitcode-stdio-1.0
@@ -105,37 +106,38 @@
   build-depends:
     base,
     bitvec,
-    integer-gmp,
-    primitive >=0.5,
-    quickcheck-classes >=0.6.1,
+    ghc-bignum,
+    primitive >=0.5 && <0.10,
+    quickcheck-classes-base <0.7,
+    quickcheck-classes >=0.6.1 && <0.7,
     vector >=0.11,
-    tasty,
-    tasty-hunit,
-    tasty-quickcheck
-  if impl(ghc <8.0)
-    build-depends:
-      semigroups >=0.8
+    tasty <1.6,
+    tasty-quickcheck <0.12
   default-language: Haskell2010
   hs-source-dirs: test
   other-modules:
     Support
     Tests.Conc
+    Tests.F2Poly
     Tests.MVector
     Tests.MVectorTS
     Tests.SetOps
     Tests.SetOpsTS
     Tests.Vector
-  ghc-options: -Wall -threaded -rtsopts
+  ghc-options: -Wall -Wcompat
+  if !arch(wasm32)
+    ghc-options: -threaded -rtsopts
   include-dirs: test
 
-benchmark gauge
+benchmark bitvec-bench
   build-depends:
     base,
     bitvec,
-    containers,
-    gauge,
-    integer-gmp,
-    random,
+    containers <0.9,
+    ghc-bignum,
+    random <1.4,
+    tasty,
+    tasty-bench >=0.4 && <0.6,
     vector
   type: exitcode-stdio-1.0
   main-is: Bench.hs
@@ -143,6 +145,7 @@
   hs-source-dirs: bench
   other-modules:
     Bench.BitIndex
+    Bench.Common
     Bench.GCD
     Bench.Invert
     Bench.Intersection
@@ -154,4 +157,4 @@
     Bench.Reverse
     Bench.Sum
     Bench.Union
-  ghc-options: -O2 -Wall
+  ghc-options: -O2 -Wall -Wcompat
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 <inttypes.h>
+#include <stddef.h>
+
+#ifdef __x86_64__
+#include <immintrin.h>
+#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,47 @@
+# 1.1.6.0
+
+* Make `basicOverlaps` precise instead of rounding bounds up to the next `Word`.
+
+# 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.
+* Tweak inlining pragmas to inline less aggressively.
+
+# 1.1.3.0
+
+* Fix malformed `signum` for `F2Poly`.
+
+# 1.1.2.0
+
+* Fix `setBit`, `clearBit`, `complementBit` to preserve vector's length.
+* Fix various issues on big-endian architectures.
+* Fix Cabal 3.7+ incompatibility.
+
+# 1.1.1.0
+
+* Export `BitVec` and `BitMVec` constructors.
+
+# 1.1.0.0
+
+* Fix a grave bug in `bitIndex`.
+* Remove `integer-gmp` flag.
+* Make `libgmp` flag disabled by default.
+  Users are strongly encouraged to enable it whenever possible.
+* Add `mapBits` and `mapInPlace` functions.
+* Add `cloneToByteString` and `cloneFromByteString` functions.
+
 # 1.0.3.0
 
 * Add `Bits (Vector Bit)` instance.
diff --git a/src/Data/Bit.hs b/src/Data/Bit.hs
--- a/src/Data/Bit.hs
+++ b/src/Data/Bit.hs
@@ -1,27 +1,41 @@
 {-# LANGUAGE CPP #-}
+{-# OPTIONS_HADDOCK prune #-}
 
 #ifndef BITVEC_THREADSAFE
 -- |
 -- Module:      Data.Bit
--- Copyright:   (c) 2019 Andrew Lelechenko, 2012-2016 James Cook
+-- Copyright:   (c) 2019-2022 Andrew Lelechenko, 2012-2016 James Cook
 -- Licence:     BSD3
 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
 --
--- This module exposes an interface with thread-unsafe writes and flips.
--- Consider using "Data.Bit.ThreadSafe", which is thread-safe, but slower (up to 20%).
+-- 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).
+--
+-- @since 0.1
 module Data.Bit
 #else
 -- |
 -- Module:      Data.Bit.ThreadSafe
--- Copyright:   (c) 2019 Andrew Lelechenko, 2012-2016 James Cook
+-- Copyright:   (c) 2019-2022 Andrew Lelechenko, 2012-2016 James Cook
 -- Licence:     BSD3
 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
 --
 -- This module exposes an interface with thread-safe writes and flips.
--- Consider using "Data.Bit", which is faster (up to 20%), but thread-unsafe.
+-- 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 not thread-safe.
+--
+-- @since 1.0
 module Data.Bit.ThreadSafe
 #endif
   ( Bit(..)
+  , U.Vector(BitVec)
+  , U.MVector(BitMVec)
 
   , unsafeFlipBit
   , flipBit
@@ -35,8 +49,12 @@
   , castToWords8
   , cloneToWords8
 
+  , cloneFromByteString
+  , cloneToByteString
+
   -- * Immutable operations
   , zipBits
+  , mapBits
   , invertBits
   , reverseBits
   , bitIndex
@@ -53,6 +71,7 @@
 
   -- * Mutable operations
   , zipInPlace
+  , mapInPlace
   , invertInPlace
   , reverseInPlace
   , selectBitsInPlace
@@ -66,6 +85,7 @@
   ) where
 
 import Prelude hiding (and, or)
+import qualified Data.Vector.Unboxed as U
 
 #ifndef BITVEC_THREADSAFE
 import Data.Bit.F2Poly
diff --git a/src/Data/Bit/F2Poly.hs b/src/Data/Bit/F2Poly.hs
--- a/src/Data/Bit/F2Poly.hs
+++ b/src/Data/Bit/F2Poly.hs
@@ -1,10 +1,8 @@
 {-# LANGUAGE CPP                        #-}
 
 {-# LANGUAGE BangPatterns               #-}
-{-# LANGUAGE DeriveDataTypeable         #-}
 {-# LANGUAGE DeriveGeneric              #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE LambdaCase                 #-}
 {-# LANGUAGE MagicHash                  #-}
 {-# LANGUAGE RankNTypes                 #-}
 
@@ -37,72 +35,88 @@
 import Data.Char
 import Data.Coerce
 import Data.Primitive.ByteArray
-import Data.Typeable
+import qualified Data.Vector.Primitive as P
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
+import GHC.Exts
 import GHC.Generics
+import GHC.Num.BigNat
+import GHC.Num.Integer
 import Numeric
 
-#if UseIntegerGmp
-import qualified Data.Vector.Primitive as P
-import GHC.Exts
-import GHC.Integer.GMP.Internals
-import GHC.Integer.Logarithms
-import Unsafe.Coerce
-#endif
-
 -- | Binary polynomials of one variable, backed
--- by an unboxed 'Data.Vector.Unboxed.Vector' 'Bit'.
+-- by an unboxed 'Data.Vector.Unboxed.Vector' t'Bit'.
 --
 -- Polynomials are stored normalized, without leading zero coefficients.
 --
--- 'Ord' instance does not make much sense mathematically,
+-- The 'Ord' instance does not make much sense mathematically,
 -- it is defined only for the sake of 'Data.Set.Set', 'Data.Map.Map', etc.
 --
 -- >>> :set -XBinaryLiterals
--- >>> -- (1 + x) (1 + x + x^2) = 1 + x^3 (mod 2)
+-- >>> -- (1 + x) * (1 + x + x^2) = 1 + x^3 (mod 2)
 -- >>> 0b11 * 0b111 :: F2Poly
 -- 0b1001
+--
+-- @since 1.0.1.0
 newtype F2Poly = F2Poly {
   unF2Poly :: U.Vector Bit
-  -- ^ Convert 'F2Poly' to a vector of coefficients
+  -- ^ Convert an t'F2Poly' to a vector of coefficients
   -- (first element corresponds to a constant term).
+  --
+  -- >>> :set -XBinaryLiterals
+  -- >>> unF2Poly 0b1101
+  -- [1,0,1,1]
+  --
+  -- @since 1.0.1.0
   }
-  deriving (Eq, Ord, Typeable, Generic, NFData)
+  deriving (Eq, Ord, Generic, NFData)
 
--- | Make 'F2Poly' from a list of coefficients
+-- | Make an t'F2Poly' from a list of coefficients
 -- (first element corresponds to a constant term).
+--
+-- >>> :set -XOverloadedLists
+-- >>> toF2Poly [1,0,1,1,0,0]
+-- 0b1101
+--
+-- @since 1.0.1.0
 toF2Poly :: U.Vector Bit -> F2Poly
 toF2Poly xs = F2Poly $ dropWhileEnd $ castFromWords $ cloneToWords xs
 
--- | Valid 'F2Poly' has offset 0 and no trailing garbage.
-_isValid :: F2Poly -> Bool
-_isValid (F2Poly (BitVec o l arr)) = o == 0 && l == l'
-  where
-    l' = U.length $ dropWhileEnd $ BitVec 0 (sizeofByteArray arr `shiftL` 3) arr
+zero :: F2Poly
+zero = F2Poly $ BitVec 0 0 $
+  ByteArray (unBigNat bigNatZero)
 
+one :: F2Poly
+one = F2Poly $ BitVec 0 1 $
+  ByteArray (unBigNat bigNatOne)
+
+-- -- | A valid t'F2Poly' has offset 0 and no trailing garbage.
+-- _isValid :: F2Poly -> Bool
+-- _isValid (F2Poly (BitVec o l arr)) = o == 0 && l == l'
+--   where
+--     l' = U.length $ dropWhileEnd $ BitVec 0 (sizeofByteArray arr `shiftL` 3) arr
+
 -- | Addition and multiplication are evaluated modulo 2.
 --
 -- 'abs' = 'id' and 'signum' = 'const' 1.
 --
 -- 'fromInteger' converts a binary polynomial, encoded as 'Integer',
--- to 'F2Poly' encoding.
+-- to t'F2Poly' encoding.
 instance Num F2Poly where
   (+) = coerce xorBits
   (-) = coerce xorBits
   negate = id
   abs    = id
-  signum = const (F2Poly (U.singleton (Bit True)))
+  signum = const one
   (*) = coerce ((dropWhileEnd .) . karatsuba)
-#if UseIntegerGmp
   fromInteger !n = case n of
-    S# i#   -> F2Poly $ BitVec 0 (wordSize - I# (word2Int# (clz# (int2Word# i#))))
-                      $ fromBigNat $ wordToBigNat (int2Word# i#)
-    Jp# bn# -> F2Poly $ BitVec 0 (I# (integerLog2# n) + 1) $ fromBigNat bn#
-    Jn#{}   -> error "F2Poly.fromInteger: argument must be non-negative"
-#else
-  fromInteger = F2Poly . dropWhileEnd . integerToBits
-#endif
+    IS i#
+      | n < 0     -> throw Underflow
+      | otherwise -> F2Poly $ BitVec 0 (wordSize - I# (word2Int# (clz# (int2Word# i#))))
+                     $ ByteArray (bigNatFromWord# (int2Word# i#))
+    IP bn# -> F2Poly $ BitVec 0 (I# (word2Int# (integerLog2# n)) + 1) $ ByteArray bn#
+    IN{}   -> throw Underflow
+  {-# INLINE fromInteger #-}
 
   {-# INLINE (+)         #-}
   {-# INLINE (-)         #-}
@@ -110,24 +124,19 @@
   {-# INLINE abs         #-}
   {-# INLINE signum      #-}
   {-# INLINE (*)         #-}
-  {-# INLINE fromInteger #-}
 
 instance Enum F2Poly where
   fromEnum = fromIntegral
-#if UseIntegerGmp
-  toEnum !(I# i#) = F2Poly $ BitVec 0 (wordSize - I# (word2Int# (clz# (int2Word# i#))))
-                           $ fromBigNat $ wordToBigNat (int2Word# i#)
-#else
-  toEnum = fromIntegral
-#endif
+  toEnum (I# i#) = F2Poly $ BitVec 0 (wordSize - I# (word2Int# (clz# (int2Word# i#))))
+                           $ ByteArray (bigNatFromWord# (int2Word# i#))
 
 instance Real F2Poly where
   toRational = fromIntegral
 
--- | 'toInteger' converts a binary polynomial, encoded as 'F2Poly',
--- to 'Integer' encoding.
+-- | 'toInteger' converts a binary polynomial, encoded as t'F2Poly',
+-- to an 'Integer' encoding.
 instance Integral F2Poly where
-  toInteger = bitsToInteger . unF2Poly
+  toInteger xs = integerFromBigNat# (bitsToByteArray (unF2Poly xs))
   quotRem (F2Poly xs) (F2Poly ys) = (F2Poly (dropWhileEnd qs), F2Poly (dropWhileEnd rs))
     where
       (qs, rs) = quotRemBits xs ys
@@ -144,16 +153,14 @@
   -> U.Vector Bit
 xorBits (BitVec _ 0 _) ys = ys
 xorBits xs (BitVec _ 0 _) = xs
-#if UseIntegerGmp
 -- GMP has platform-dependent ASM implementations for mpn_xor_n,
 -- which are impossible to beat by native Haskell.
-xorBits (BitVec 0 lx xarr) (BitVec 0 ly yarr) = case lx `compare` ly of
+xorBits (BitVec 0 lx (ByteArray xarr)) (BitVec 0 ly (ByteArray yarr)) = case lx `compare` ly of
   LT -> BitVec 0 ly zs
   EQ -> dropWhileEnd $ BitVec 0 (lx `min` (sizeofByteArray zs `shiftL` 3)) zs
   GT -> BitVec 0 lx zs
   where
-    zs = fromBigNat (toBigNat xarr `xorBigNat` toBigNat yarr)
-#endif
+    zs = ByteArray (xarr `bigNatXor` yarr)
 xorBits xs ys = dropWhileEnd $ runST $ do
   let lx = U.length xs
       ly = U.length ys
@@ -170,33 +177,31 @@
 
 karatsuba :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
 karatsuba xs ys
-  | karatsubaThreshold < 2 * wordSize
-  = error $ "karatsubaThreshold must be >= " ++ show (2 * wordSize)
   | xs == ys = sqrBits xs
   | lenXs <= karatsubaThreshold || lenYs <= karatsubaThreshold
   = mulBits xs ys
   | otherwise = runST $ do
     zs <- MU.unsafeNew lenZs
-    forM_ [0, wordSize .. lenZs - 1] $ \k -> do
+    forM_ [0 .. divWordSize (lenZs - 1)] $ \k -> do
       let z0  = indexWord0 zs0   k
           z11 = indexWord0 zs11 (k - m)
           z10 = indexWord0 zs0  (k - m)
           z12 = indexWord0 zs2  (k - m)
           z2  = indexWord0 zs2  (k - 2 * m)
-      writeWord zs k (z0 `xor` z11 `xor` z10 `xor` z12 `xor` z2)
+      writeWord zs (mulWordSize k) (z0 `xor` z11 `xor` z10 `xor` z12 `xor` z2)
     U.unsafeFreeze zs
   where
     lenXs = U.length xs
     lenYs = U.length ys
     lenZs = lenXs + lenYs - 1
 
-    m'    = ((lenXs `min` lenYs) + 1) `quot` 2
-    m     = m' - modWordSize m'
+    m    = (min lenXs lenYs + 1) `unsafeShiftR` (lgWordSize + 1)
+    m'   = mulWordSize m
 
-    xs0  = U.unsafeSlice 0 m xs
-    xs1  = U.unsafeSlice m (lenXs - m) xs
-    ys0  = U.unsafeSlice 0 m ys
-    ys1  = U.unsafeSlice m (lenYs - m) ys
+    xs0  = U.unsafeSlice 0 m' xs
+    xs1  = U.unsafeSlice m' (lenXs - m') xs
+    ys0  = U.unsafeSlice 0 m' ys
+    ys1  = U.unsafeSlice m' (lenYs - m') ys
 
     xs01 = xorBits xs0 xs1
     ys01 = xorBits ys0 ys1
@@ -205,17 +210,14 @@
     zs11 = karatsuba xs01 ys01
 
 indexWord0 :: U.Vector Bit -> Int -> Word
-indexWord0 bv i
-  | i <= - wordSize         = 0
-  | lenI <= 0               = 0
-  | i < 0, lenI >= wordSize = word0
-  | i < 0                   = word0 .&. loMask lenI
-  | lenI >= wordSize        = word
-  | otherwise               = word .&. loMask lenI
+indexWord0 bv i'
+  | i < 0 || lenI <= 0 = 0
+  | lenI >= wordSize   = word
+  | otherwise          = word .&. loMask lenI
   where
+    i     = mulWordSize i'
     lenI  = U.length bv - i
     word  = indexWord bv i
-    word0 = indexWord bv 0 `unsafeShiftL` (- i)
 
 mulBits :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
 mulBits xs ys
@@ -278,55 +280,25 @@
         0 -> go (n - wordSize)
         w -> n - countLeadingZeros w
 
-#if UseIntegerGmp
-
 bitsToByteArray :: U.Vector Bit -> ByteArray#
 bitsToByteArray xs = arr
   where
-    ys = if U.null xs then U.singleton 0 else cloneToWords xs
+    ys = if U.null xs then U.singleton (0 :: Word) else cloneToWords xs
     !(P.Vector _ _ (ByteArray arr)) = toPrimVector ys
 
-fromBigNat :: BigNat -> ByteArray
-fromBigNat = unsafeCoerce
--- fromBigNat (BN# arr) = ByteArray arr
-
-toBigNat :: ByteArray -> BigNat
-toBigNat = unsafeCoerce
--- toBigNat (ByteArray arr) = BN# arr
-
-bitsToInteger :: U.Vector Bit -> Integer
-bitsToInteger xs = bigNatToInteger (BN# (bitsToByteArray xs))
-
-#else
-
-integerToBits :: Integer -> U.Vector Bit
-integerToBits x = U.generate (bitLen x) (Bit . testBit x)
-
-bitLen :: Integer -> Int
-bitLen x
-  = fst
-  $ head
-  $ dropWhile (\(_, b) -> x >= b)
-  $ map (\a -> (a, 1 `shiftL` a))
-  $ map (1 `shiftL`)
-  $ [lgWordSize..]
-
-bitsToInteger :: U.Vector Bit -> Integer
-bitsToInteger = U.ifoldl' (\acc i (Bit b) -> if b then acc `setBit` i else acc) 0
-
-#endif
-
 -- | Execute the extended Euclidean algorithm.
 -- For polynomials @a@ and @b@, compute their unique greatest common divisor @g@
--- and the unique coefficient polynomial @s@ satisfying @as + bt = g@.
+-- and the unique coefficient polynomial @s@ satisfying \( a \cdot s + b \cdot t = g \).
 --
 -- >>> :set -XBinaryLiterals
 -- >>> gcdExt 0b101 0b0101
 -- (0b101,0b0)
 -- >>> gcdExt 0b11 0b111
 -- (0b1,0b10)
+--
+-- @since 1.0.2.0
 gcdExt :: F2Poly -> F2Poly -> (F2Poly, F2Poly)
-gcdExt = go 1 0
+gcdExt = go one zero
   where
     go s s' r r'
       | r' == 0   = (r, s)
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,146 +0,0 @@
-{-# LANGUAGE CPP              #-}
-{-# LANGUAGE MagicHash        #-}
-{-# LANGUAGE UnliftedFFITypes #-}
-
-#if UseLibGmp
-
-module Data.Bit.Gmp
-  ( mpnCom
-  , mpnLshift
-  , mpnRshift
-  , mpnScan0
-  , mpnScan1
-  , 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_lshift"
-  mpn_lshift :: MutableByteArray# s -> ByteArray# -> Int# -> Word# -> IO Word
-
-mpnLshift :: MutableByteArray s -> ByteArray -> Int -> Word -> ST s Word
-mpnLshift (MutableByteArray res#) (ByteArray arg#) (I# limbs#) (W# count#) =
-  unsafeIOToST (mpn_lshift res# arg# limbs# count#)
-{-# INLINE mpnLshift #-}
-
-foreign import ccall unsafe "__gmpn_rshift"
-  mpn_rshift :: MutableByteArray# s -> ByteArray# -> Int# -> Word# -> IO Word
-
-mpnRshift :: MutableByteArray s -> ByteArray -> Int -> Word -> ST s Word
-mpnRshift (MutableByteArray res#) (ByteArray arg#) (I# limbs#) (W# count#) =
-  unsafeIOToST (mpn_rshift res# arg# limbs# count#)
-{-# INLINE mpnRshift #-}
-
-foreign import ccall unsafe "__gmpn_scan0"
-  mpn_scan0 :: ByteArray# -> Word# -> IO Word
-
-mpnScan0 :: ByteArray -> Word -> Word
-mpnScan0 (ByteArray arg#) (W# bit#) =
-  unsafeDupablePerformIO (mpn_scan0 arg# bit#)
-{-# INLINE mpnScan0 #-}
-
-foreign import ccall unsafe "__gmpn_scan1"
-  mpn_scan1 :: ByteArray# -> Word# -> IO Word
-
-mpnScan1 :: ByteArray -> Word -> Word
-mpnScan1 (ByteArray arg#) (W# bit#) =
-  unsafeDupablePerformIO (mpn_scan1 arg# bit#)
-{-# INLINE mpnScan1 #-}
-
-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 #-}
-
-#else
-
-module Data.Bit.Gmp where
-
-#endif
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,12 +1,14 @@
-{-# LANGUAGE CPP                  #-}
-
 {-# LANGUAGE BangPatterns         #-}
+{-# LANGUAGE BinaryLiterals       #-}
+{-# LANGUAGE CPP                  #-}
 {-# LANGUAGE FlexibleContexts     #-}
+{-# LANGUAGE MagicHash            #-}
 {-# LANGUAGE RankNTypes           #-}
 {-# LANGUAGE FlexibleInstances    #-}
 {-# LANGUAGE UndecidableInstances #-}
 
-{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# OPTIONS_GHC -Wno-orphans #-}
+{- HLINT ignore "Unused LANGUAGE pragma" -}
 
 #ifndef BITVEC_THREADSAFE
 module Data.Bit.Immutable
@@ -21,7 +23,11 @@
   , castToWords8
   , cloneToWords8
 
+  , cloneFromByteString
+  , cloneToByteString
+
   , zipBits
+  , mapBits
   , invertBits
   , selectBits
   , excludeBits
@@ -33,11 +39,13 @@
   , listBits
   ) where
 
+#include "MachDeps.h"
+
 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
@@ -48,23 +56,22 @@
 #endif
 import Data.Bit.PdepPext
 import Data.Bit.Utils
+import qualified Data.ByteString.Internal as BS
 import Data.Primitive.ByteArray
 import qualified Data.Vector.Primitive as P
+import qualified Data.Vector.Storable as S
 import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Base as UB
 import qualified Data.Vector.Unboxed.Mutable as MU
 import Data.Word
-import Unsafe.Coerce
 
-#include "MachDeps.h"
-
-#if WORD_SIZE_IN_BITS == 64
-#define GMP_LIMB_SHIFT 3
-#elif WORD_SIZE_IN_BITS == 32
-#define GMP_LIMB_SHIFT 2
-#else
-#error unsupported WORD_SIZE_IN_BITS config
+#ifdef WORDS_BIGENDIAN
+import GHC.Exts
 #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 (.|.)
@@ -80,6 +87,27 @@
     | n < 0 || n >= U.length v = False
     | otherwise = unBit (U.unsafeIndex v n)
 
+  setBit v n
+    | n < 0 || n >= U.length v = v
+    | otherwise = runST $ do
+      u <- U.thaw v
+      MU.unsafeWrite u n (Bit True)
+      U.unsafeFreeze u
+
+  clearBit v n
+    | n < 0 || n >= U.length v = v
+    | otherwise = runST $ do
+      u <- U.thaw v
+      MU.unsafeWrite u n (Bit False)
+      U.unsafeFreeze u
+
+  complementBit v n
+    | n < 0 || n >= U.length v = v
+    | otherwise = runST $ do
+      u <- U.thaw v
+      unsafeFlipBit u n
+      U.unsafeFreeze u
+
   bit n
     | n < 0 = U.empty
     | otherwise = runST $ do
@@ -118,8 +146,11 @@
 -- to an unboxed vector of bits.
 -- Cf. 'Data.Bit.castFromWordsM'.
 --
--- >>> castFromWords (Data.Vector.Unboxed.singleton 123)
+-- >>> :set -XOverloadedLists
+-- >>> castFromWords [123]
 -- [1,1,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
+--
+-- @since 1.0.0.0
 castFromWords :: U.Vector Word -> U.Vector Bit
 castFromWords ws = BitVec (mulWordSize off) (mulWordSize len) arr
   where
@@ -127,11 +158,13 @@
 
 -- | Try to cast an unboxed vector of bits
 -- to an unboxed vector of words.
--- It succeeds if a vector of bits is aligned.
+-- It succeeds if the vector of bits is aligned.
 -- Use 'cloneToWords' otherwise.
 -- Cf. 'Data.Bit.castToWordsM'.
 --
--- prop> castToWords (castFromWords v) == Just v
+-- > castToWords (castFromWords v) == Just v
+--
+-- @since 1.0.0.0
 castToWords :: U.Vector Bit -> Maybe (U.Vector Word)
 castToWords (BitVec s n ws)
   | aligned s, aligned n =
@@ -145,135 +178,280 @@
 -- the last word will be zero-padded.
 -- Cf. 'Data.Bit.cloneToWordsM'.
 --
--- >>> cloneToWords (read "[1,1,0,1,1,1,1,0]")
+-- >>> :set -XOverloadedLists
+-- >>> cloneToWords [1,1,0,1,1,1,1]
 -- [123]
+--
+-- @since 1.0.0.0
 cloneToWords :: U.Vector Bit -> U.Vector Word
 cloneToWords v = runST $ do
   v' <- U.unsafeThaw v
   w  <- cloneToWordsM v'
   U.unsafeFreeze w
-{-# INLINE cloneToWords #-}
+{-# INLINABLE cloneToWords #-}
 
--- | Cast a unboxed vector of 'Word8'
+-- | Cast an unboxed vector of 'Word8'
 -- to an unboxed vector of bits.
--- This can be used in conjunction
--- with @bytestring-to-vector@ package
--- to convert from 'Data.ByteString.ByteString':
 --
--- >>> :set -XOverloadedStrings
--- >>> import Data.Vector.Storable.ByteString
--- >>> castFromWords8 (Data.Vector.convert (byteStringToVector "abc"))
--- [1,0,0,0,0,1,1,0,0,1,0,0,0,1,1,0,1,1,0,0,0,1,1,0]
+-- On big-endian architectures 'castFromWords8'
+-- resorts to copying instead of aliasing the underlying array.
+--
+-- >>> :set -XOverloadedLists
+-- >>> castFromWords8 [123]
+-- [1,1,0,1,1,1,1,0]
+--
+-- @since 1.0.3.0
 castFromWords8 :: U.Vector Word8 -> U.Vector Bit
 castFromWords8 ws = BitVec (off `shiftL` 3) (len `shiftL` 3) arr
   where
-    P.Vector off len arr = unsafeCoerce ws
+#ifdef WORDS_BIGENDIAN
+    UB.V_Word8 (P.Vector off' len arr') = ws
+    off = 0
+    arr = runST $ do
+      let lenWords = nWords $ len `shiftL` 3
+          len' = wordsToBytes lenWords
+      marr <- newByteArray len'
+      copyByteArray marr 0 arr' off' len
+      fillByteArray marr len (len' - len) 0
+      forM_ [0..lenWords - 1] $ \i -> do
+        W# w <- readByteArray marr i
+        writeByteArray marr i (W# (byteSwap# w))
+      unsafeFreezeByteArray marr
+#else
+    UB.V_Word8 (P.Vector off len arr) = ws
+#endif
 
 -- | Try to cast an unboxed vector of bits
 -- to an unboxed vector of 'Word8'.
--- It succeeds if a vector of bits is aligned.
+-- It succeeds if the vector of bits is aligned.
 -- Use 'Data.Bit.cloneToWords8' otherwise.
 --
--- prop> castToWords8 (castFromWords8 v) == Just v
+-- > castToWords8 (castFromWords8 v) == Just v
+--
+-- @since 1.0.3.0
 castToWords8 :: U.Vector Bit -> Maybe (U.Vector Word8)
+#ifdef WORDS_BIGENDIAN
+castToWords8 = const Nothing
+#else
 castToWords8 (BitVec s n ws)
-  | s .&. 7 == 0, n .&. 7 == 0 =
-    Just $ unsafeCoerce $ P.Vector (s `shiftR` 3) (n `shiftR` 3) ws
+  | s .&. 7 == 0, n .&. 7 == 0
+  = Just $ UB.V_Word8 $ P.Vector (s `shiftR` 3) (n `shiftR` 3) ws
   | otherwise = Nothing
+#endif
 
 -- | Clone an unboxed vector of bits
 -- to a new unboxed vector of 'Word8'.
--- If the bits don't completely fill the words,
+-- If the bits don't completely fill the bytes,
 -- the last 'Word8' will be zero-padded.
--- This can be used in conjunction
--- with @bytestring-to-vector@ package
--- to convert to 'Data.ByteString.ByteString':
 --
 -- >>> :set -XOverloadedLists
--- >>> import Data.Vector.Storable.ByteString
--- >>> vectorToByteString (Data.Vector.convert (Data.Bit.cloneToWords8 [1,0,0,0,0,1,1,0,0,1,0,0,0,1,1,0,1,1,0,0,0,1]))
--- "ab#"
+-- >>> cloneToWords8 [1,1,0,1,1,1,1]
+-- [123]
+--
+-- @since 1.0.3.0
 cloneToWords8 :: U.Vector Bit -> U.Vector Word8
 cloneToWords8 v = runST $ do
   v' <- U.unsafeThaw v
   w  <- cloneToWords8M v'
   U.unsafeFreeze w
-{-# INLINE cloneToWords8 #-}
+{-# INLINABLE cloneToWords8 #-}
 
+-- | Clone a 'BS.ByteString' to a new unboxed vector of bits.
+--
+-- >>> :set -XOverloadedStrings
+-- >>> cloneFromByteString "abc"
+-- [1,0,0,0,0,1,1,0,0,1,0,0,0,1,1,0,1,1,0,0,0,1,1,0]
+--
+-- @since 1.1.0.0
+cloneFromByteString :: BS.ByteString -> U.Vector Bit
+cloneFromByteString
+  = castFromWords8
+  . U.convert
+  . uncurry3 S.unsafeFromForeignPtr
+  . BS.toForeignPtr
+
+-- | Clone an unboxed vector of bits to a new 'BS.ByteString'.
+-- If the bits don't completely fill the bytes,
+-- the last character will be zero-padded.
+--
+-- >>> :set -XOverloadedLists
+-- >>> cloneToByteString [1,0,0,0,0,1,1,0,0,1,0,0,0,1,1,0,1,1,0,0,0,1]
+-- "ab#"
+--
+-- @since 1.1.0.0
+cloneToByteString :: U.Vector Bit -> BS.ByteString
+cloneToByteString
+  = uncurry3 BS.fromForeignPtr
+  . S.unsafeToForeignPtr
+  . U.convert
+  . cloneToWords8
+
+uncurry3 :: (a -> b -> c -> d) -> (a, b, c) -> d
+uncurry3 f (x, y, z) = f x y z
+
 -- | 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.
 --
+-- 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 (.&.) (read "[1,1,0]") (read "[0,1,1]") -- intersection
+-- >>> zipBits (.&.) [1,1,0] [0,1,1] -- intersection
 -- [0,1,0]
--- >>> zipBits (.|.) (read "[1,1,0]") (read "[0,1,1]") -- union
+-- >>> zipBits (.|.) [1,1,0] [0,1,1] -- union
 -- [1,1,1]
--- >>> zipBits (\x y -> x .&. complement y) (read "[1,1,0]") (read "[0,1,1]") -- difference
+-- >>> zipBits (\x y -> x .&. complement y) [1,1,0] [0,1,1] -- difference
 -- [1,0,0]
--- >>> zipBits xor (read "[1,1,0]") (read "[0,1,1]") -- symmetric difference
+-- >>> zipBits xor [1,1,0] [0,1,1] -- symmetric difference
 -- [1,0,1]
+--
+-- @since 1.0.0.0
 zipBits
   :: (forall a . Bits a => a -> a -> a)
   -> 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` GMP_LIMB_SHIFT
+        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
+  _ -> 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.
+--
+-- >>> :set -XOverloadedLists
+-- >>> import Data.Bits
+-- >>> mapBits complement [0,1,1]
+-- [1,0,0]
+--
+-- @since 1.1.0.0
+mapBits
+  :: (forall a . Bits a => a -> a)
+  -> U.Vector Bit
+  -> U.Vector Bit
+mapBits f = case (unBit (f (Bit False)), unBit (f (Bit True))) of
+  (False, False) -> (`U.replicate` Bit False) . U.length
+  (False, True)  -> id
+  (True, False)  -> invertBits
+  (True, True)   -> (`U.replicate` Bit True) . U.length
+{-# INLINE mapBits #-}
+
 -- | Invert (flip) all bits.
 --
--- >>> invertBits (read "[0,1,0,1,0]")
+-- >>> :set -XOverloadedLists
+-- >>> invertBits [0,1,0,1,0]
 -- [1,0,1,0,1]
+--
+-- @since 1.0.1.0
 invertBits
   :: 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` GMP_LIMB_SHIFT)
-  mpnCom brr arg w
+      b = wordsToBytes w
+  brr <- newByteArray b
+  ompCom brr arg b
   BitVec 0 l <$> unsafeFreezeByteArray brr
 #endif
 invertBits xs = runST $ do
@@ -283,35 +461,62 @@
     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 deposit instruction (PDEP).
+-- to the result. Similar to the
+-- [parallel bit extract instruction (PEXT)](https://en.wikipedia.org/wiki/X86_Bit_manipulation_instruction_set#Parallel_bit_deposit_and_extract).
 --
--- >>> selectBits (read "[0,1,0,1,1]") (read "[1,1,0,0,1]")
+-- 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]
 --
 -- Here is a reference (but slow) implementation:
 --
 -- > import qualified Data.Vector.Unboxed as U
--- > selectBits mask ws == U.map snd (U.filter (unBit . fst) (U.zip mask ws))
+-- > selectBits mask ws = U.map snd (U.filter (unBit . fst) (U.zip mask ws))
+--
+-- @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.
 --
--- >>> excludeBits (read "[0,1,0,1,1]") (read "[1,1,0,0,1]")
+-- 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]
 --
 -- Here is a reference (but slow) implementation:
 --
 -- > import qualified Data.Vector.Unboxed as U
--- > excludeBits mask ws == U.map snd (U.filter (not . unBit . fst) (U.zip mask ws))
+-- > excludeBits mask ws = U.map snd (U.filter (not . unBit . fst) (U.zip mask ws))
+--
+-- @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
@@ -319,12 +524,21 @@
 
 -- | Reverse the order of bits.
 --
--- >>> reverseBits (read "[1,1,0,1,0]")
+-- >>> :set -XOverloadedLists
+-- >>> reverseBits [1,1,0,1,0]
 -- [0,1,0,1,1]
 --
--- Consider using @vector-rotcev@ package
+-- Consider using the [vector-rotcev](https://hackage.haskell.org/package/vector-rotcev) package
 -- to reverse vectors in O(1) time.
+--
+-- @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
@@ -352,19 +566,27 @@
 -- with the specified value, if any.
 -- Similar to 'Data.Vector.Unboxed.elemIndex', but up to 64x faster.
 --
--- >>> bitIndex (Bit True) (read "[0,0,1,0,1]")
+-- >>> :set -XOverloadedLists
+-- >>> bitIndex 1 [0,0,1,0,1]
 -- Just 2
--- >>> bitIndex (Bit True) (read "[0,0,0,0,0]")
+-- >>> bitIndex 1 [0,0,0,0,0]
 -- Nothing
 --
--- prop> bitIndex bit == nthBitIndex bit 1
+-- > bitIndex bit == nthBitIndex bit 1
 --
 -- One can also use it to reduce a vector with disjunction or conjunction:
 --
--- >>> import Data.Maybe
--- >>> isAnyBitSet   = isJust    . bitIndex (Bit True)
--- >>> areAllBitsSet = isNothing . bitIndex (Bit False)
+-- > import Data.Maybe
+-- > isAnyBitSet   = isJust    . bitIndex 1
+-- > areAllBitsSet = isNothing . bitIndex 0
+--
+-- @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
@@ -422,30 +644,38 @@
     | n >= off + len = Nothing
     | otherwise = case ffs (indexByteArray arr n) of
       Nothing  -> go (n + 1)
-      r@Just{} -> r
+      Just r  -> Just $ mulWordSize (n - off) + r
 bitIndexInWords (Bit False) !off !len !arr = go off
  where
   go !n
     | n >= off + len = Nothing
     | otherwise = case ffs (complement (indexByteArray arr n)) of
-      Nothing  -> go (n + 1)
-      r@Just{} -> r
+      Nothing -> go (n + 1)
+      Just r  -> Just $ mulWordSize (n - off) + r
 
 -- | Return the index of the @n@-th bit in the vector
 -- with the specified value, if any.
 -- Here @n@ is 1-based and the index is 0-based.
 -- Non-positive @n@ results in an error.
 --
--- >>> nthBitIndex (Bit True) 2 (read "[0,1,0,1,1,1,0]")
+-- >>> :set -XOverloadedLists
+-- >>> nthBitIndex 1 2 [0,1,0,1,1,1,0] -- 2nd occurence of 1
 -- Just 3
--- >>> nthBitIndex (Bit True) 5 (read "[0,1,0,1,1,1,0]")
+-- >>> nthBitIndex 1 5 [0,1,0,1,1,1,0] -- 5th occurence of 1
 -- Nothing
 --
 -- One can use 'nthBitIndex' to implement
 -- to implement @select{0,1}@ queries
 -- for <https://en.wikipedia.org/wiki/Succinct_data_structure succinct dictionaries>.
+--
+-- @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
@@ -523,17 +753,20 @@
 
 -- | Return the number of set bits in a vector (population count, popcount).
 --
--- >>> countBits (read "[1,1,0,1,0,1]")
+-- >>> :set -XOverloadedLists
+-- >>> countBits [1,1,0,1,0,1]
 -- 4
 --
 -- One can combine 'countBits' with 'Data.Vector.Unboxed.take'
 -- to implement @rank{0,1}@ queries
 -- for <https://en.wikipedia.org/wiki/Succinct_data_structure succinct dictionaries>.
+--
+-- @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)
@@ -561,10 +794,13 @@
 countBitsInWords :: P.Vector Word -> Int
 countBitsInWords = P.foldl' (\acc word -> popCount word + acc) 0
 
--- | Return the indices of set bits in a vector.
+-- | Return 0-based indices of set bits in a vector.
 --
--- >>> listBits (read "[1,1,0,1,0,1]")
+-- >>> :set -XOverloadedLists
+-- >>> listBits [1,1,0,1,0,1]
 -- [0,1,3,5]
+--
+-- @since 0.1
 listBits :: U.Vector Bit -> [Int]
 listBits (BitVec _ 0 _)                      = []
 listBits (BitVec off len arr) | offBits == 0 = case modWordSize len of
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
@@ -1,15 +1,13 @@
-{-# LANGUAGE CPP                        #-}
-
 {-# LANGUAGE BangPatterns               #-}
-{-# LANGUAGE DeriveDataTypeable         #-}
+{-# LANGUAGE CPP                        #-}
 {-# LANGUAGE DeriveGeneric              #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
-{-# LANGUAGE LambdaCase                 #-}
 {-# LANGUAGE MagicHash                  #-}
 {-# LANGUAGE MultiParamTypeClasses      #-}
 {-# LANGUAGE TypeFamilies               #-}
 {-# LANGUAGE UnboxedTuples              #-}
 {-# LANGUAGE ViewPatterns               #-}
+{- HLINT ignore "Unused LANGUAGE pragma" -}
 
 #ifndef BITVEC_THREADSAFE
 module Data.Bit.Internal
@@ -24,24 +22,20 @@
   , writeWord
   , unsafeFlipBit
   , flipBit
-  , WithInternals(..)
   , modifyByteArray
   ) where
 
-#include "vector.h"
-
 import Control.DeepSeq
 import Control.Exception
-import Control.Monad
 import Control.Monad.Primitive
 import Control.Monad.ST
 import Data.Bits
 import Data.Bit.Utils
 import Data.Primitive.ByteArray
 import Data.Ratio
-import Data.Typeable
 import qualified Data.Vector.Generic as V
 import qualified Data.Vector.Generic.Mutable as MV
+import Data.Vector.Internal.Check (checkIndex, Checks(..))
 import qualified Data.Vector.Unboxed as U
 import GHC.Generics
 
@@ -51,27 +45,47 @@
 
 #ifndef BITVEC_THREADSAFE
 -- | A newtype wrapper with a custom instance
--- of "Data.Vector.Unboxed", which packs booleans
+-- for "Data.Vector.Unboxed", which packs booleans
 -- as efficient as possible (8 values per byte).
--- Vectors of `Bit` use 8x less memory
--- than vectors of 'Bool' (which stores one value per byte).
--- but random writes are up to 10% slower.
-newtype Bit = Bit { unBit :: Bool }
-  deriving (Bounded, Enum, Eq, Ord, FiniteBits, Bits, Typeable, Generic, NFData)
+-- Unboxed vectors of t'Bit' use 8x less memory
+-- than unboxed vectors of 'Bool' (which store one value per byte),
+-- but random writes are slightly slower.
+--
+-- @since 0.1
+newtype Bit = Bit {
+  unBit :: Bool -- ^ @since 0.2.0.0
+  } deriving
+  (Bounded, Enum, Eq, Ord
+  , FiniteBits -- ^ @since 0.2.0.0
+  , Bits
+  , Generic    -- ^ @since 1.0.1.0
+  , NFData     -- ^ @since 1.0.1.0
+  )
 #else
 -- | A newtype wrapper with a custom instance
--- of "Data.Vector.Unboxed", which packs booleans
+-- for "Data.Vector.Unboxed", which packs booleans
 -- as efficient as possible (8 values per byte).
--- Vectors of `Bit` use 8x less memory
--- than vectors of 'Bool' (which stores one value per byte).
--- but random writes are up to 20% slower.
-newtype Bit = Bit { unBit :: Bool }
-  deriving (Bounded, Enum, Eq, Ord, FiniteBits, Bits, Typeable, Generic, NFData)
+-- Unboxed vectors of t'Bit' use 8x less memory
+-- than unboxed vectors of 'Bool' (which store one value per byte),
+-- but random writes are slightly slower.
+--
+-- @since 1.0.0.0
+newtype Bit = Bit {
+  unBit :: Bool -- ^ @since 0.2.0.0
+  } deriving
+  (Bounded, Enum, Eq, Ord
+  , FiniteBits -- ^ @since 0.2.0.0
+  , Bits
+  , Generic    -- ^ @since 1.0.1.0
+  , NFData     -- ^ @since 1.0.1.0
+  )
 #endif
 
 -- | There is only one lawful 'Num' instance possible
 -- with '+' = 'xor' and
--- 'fromInteger' = 'Bit' . 'odd'.
+-- 'fromInteger' = v'Bit' . 'odd'.
+--
+-- @since 1.0.1.0
 instance Num Bit where
   Bit a * Bit b = Bit (a && b)
   Bit a + Bit b = Bit (a /= b)
@@ -81,31 +95,21 @@
   signum = id
   fromInteger = Bit . odd
 
+-- | @since 1.0.1.0
 instance Real Bit where
-  toRational (Bit False) = 0
-  toRational (Bit True)  = 1
+  toRational = fromIntegral
 
+-- | @since 1.0.1.0
 instance Integral Bit where
   quotRem _ (Bit False) = throw DivideByZero
   quotRem x (Bit True)  = (x, Bit False)
-  quot    _ (Bit False) = throw DivideByZero
-  quot    x (Bit True)  = x
-  rem     _ (Bit False) = throw DivideByZero
-  rem     _ (Bit True)  = Bit False
-
-  divMod = quotRem
-  div    = quot
-  mod    = rem
-
   toInteger (Bit False) = 0
   toInteger (Bit True)  = 1
 
+-- | @since 1.0.1.0
 instance Fractional Bit where
-  fromRational x = fromInteger (numerator x) `quot` fromInteger (denominator x)
-  _ / Bit False     = throw DivideByZero
-  x / Bit True      = x
-  recip (Bit False) = throw DivideByZero
-  recip (Bit True)  = Bit True
+  fromRational x = fromInteger (numerator x) / fromInteger (denominator x)
+  (/) = quot
 
 instance Show Bit where
   showsPrec _ (Bit False) = showString "0"
@@ -123,13 +127,6 @@
 data instance U.MVector s Bit = BitMVec !Int !Int !(MutableByteArray s)
 data instance U.Vector    Bit = BitVec  !Int !Int !ByteArray
 
-newtype WithInternals = WithInternals (U.Vector Bit)
-
-#if MIN_VERSION_primitive(0,6,3)
-instance Show WithInternals where
-  show (WithInternals v@(BitVec off len ba)) = show (off, len, ba, v)
-#endif
-
 readBit :: Int -> Word -> Bit
 readBit i w = Bit (w .&. (1 `unsafeShiftL` i) /= 0)
 {-# INLINE readBit #-}
@@ -138,10 +135,10 @@
 extendToWord (Bit False) = 0
 extendToWord (Bit True ) = complement 0
 
--- | read 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 result is padded with memory garbage.
+-- | Read 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 result is padded with memory garbage.
 indexWord :: U.Vector Bit -> Int -> Word
-indexWord !(BitVec _ 0 _) _ = 0
-indexWord !(BitVec off len' arr) !i' = word
+indexWord (BitVec _ 0 _) _ = 0
+indexWord (BitVec off len' arr) !i' = word
  where
   len    = off + len'
   i      = off + i'
@@ -150,19 +147,19 @@
   loWord = indexByteArray arr loIx
   hiWord = indexByteArray arr (loIx + 1)
 
-  word   = if nMod == 0
-    then loWord
-    else if loIx == divWordSize (len - 1)
-      then (loWord `unsafeShiftR` nMod)
-      else
-        (loWord `unsafeShiftR` nMod)
-          .|. (hiWord `unsafeShiftL` (wordSize - nMod))
+  word
+    | nMod == 0
+    = loWord
+    | loIx == divWordSize (len - 1)
+    = loWord `unsafeShiftR` nMod
+    | otherwise
+    = (loWord `unsafeShiftR` nMod) .|. (hiWord `unsafeShiftL` (wordSize - nMod))
 {-# INLINE indexWord #-}
 
--- | read 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 result is padded with memory garbage.
+-- | Read 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 result is padded with memory garbage.
 readWord :: PrimMonad m => U.MVector (PrimState m) Bit -> Int -> m Word
-readWord !(BitMVec _ 0 _) _ = pure 0
-readWord !(BitMVec off len' arr) !i' = do
+readWord (BitMVec _ 0 _) _ = pure 0
+readWord (BitMVec off len' arr) !i' = do
   let len  = off + len'
       i    = off + i'
       nMod = modWordSize i
@@ -178,9 +175,7 @@
         pure
           $   (loWord `unsafeShiftR` nMod)
           .|. (hiWord `unsafeShiftL` (wordSize - nMod))
-#if __GLASGOW_HASKELL__ >= 800
 {-# SPECIALIZE readWord :: U.MVector s Bit -> Int -> ST s Word #-}
-#endif
 {-# INLINE readWord #-}
 
 modifyByteArray
@@ -201,20 +196,13 @@
     let !(# state',  _ #) = fetchAndIntArray# mba ix (word2Int# msk) state  in
     let !(# state'', _ #) = fetchOrIntArray#  mba ix (word2Int# new) state' in
     (# state'', () #)
-
--- https://gitlab.haskell.org/ghc/ghc/issues/17334
-#if __GLASGOW_HASKELL__ == 808 && __GLASGOW_HASKELL_PATCHLEVEL1__ == 1
-{-# NOINLINE modifyByteArray #-}
-#else
 {-# INLINE modifyByteArray #-}
 #endif
 
-#endif
-
--- | 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.
+-- | 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 off len' arr) !i' !x
+writeWord (BitMVec _ 0 _) !_ !_ = pure ()
+writeWord (BitMVec off len' arr) !i' !x
   | iMod == 0
   = if len >= i + wordSize
     then writeByteArray arr iDiv x
@@ -240,9 +228,7 @@
     iMod   = modWordSize i
     iDiv   = divWordSize i
 
-#if __GLASGOW_HASKELL__ >= 800
 {-# SPECIALIZE writeWord :: U.MVector s Bit -> Int -> Word -> ST s () #-}
-#endif
 {-# INLINE writeWord #-}
 
 instance MV.MVector U.MVector Bit where
@@ -258,23 +244,17 @@
 
   {-# INLINE basicUnsafeReplicate #-}
   basicUnsafeReplicate n x
-    | n < 0 =  error
-    $  "Data.Bit.basicUnsafeReplicate: negative length: "
-    ++ show n
+    | n < 0 =  error $  "Data.Bit.basicUnsafeReplicate: negative length: " ++ show n
     | otherwise = do
       arr <- newByteArray (wordsToBytes $ nWords n)
       setByteArray arr 0 (nWords n) (extendToWord x :: Word)
       pure $ BitMVec 0 n arr
 
   {-# INLINE basicOverlaps #-}
-  basicOverlaps (BitMVec i' m' arr1) (BitMVec j' n' arr2) =
+  basicOverlaps (BitMVec i m arr1) (BitMVec j n arr2) =
     sameMutableByteArray arr1 arr2
       && (between i j (j + n) || between j i (i + m))
    where
-    i = divWordSize i'
-    m = nWords (i' + m') - i
-    j = divWordSize j'
-    n = nWords (j' + n') - j
     between x y z = x >= y && x < z
 
   {-# INLINE basicLength #-}
@@ -309,11 +289,7 @@
       in  (# state', () #)
 #endif
 
-  {-# INLINE basicClear #-}
-  basicClear _ = pure ()
-
   {-# INLINE basicSet #-}
-  basicSet (BitMVec _ 0 _) _ = pure ()
   basicSet (BitMVec off len arr) (extendToWord -> x) | offBits == 0 =
     case modWordSize len of
       0    -> setByteArray arr offWords lWords (x :: Word)
@@ -343,7 +319,6 @@
     lWords   = nWords (offBits + len)
 
   {-# INLINE basicUnsafeCopy #-}
-  basicUnsafeCopy _ (BitMVec _ 0 _) = pure ()
   basicUnsafeCopy (BitMVec offDst lenDst dst) (BitMVec offSrc _ src)
     | offDstBits == 0, offSrcBits == 0 = case modWordSize lenDst of
       0 -> copyMutableByteArray dst
@@ -407,10 +382,10 @@
         x <- readWord src i
         writeWord dst i x
         do_copy (i + wordSize)
-      | otherwise = return ()
+      | otherwise = pure ()
 
   {-# INLINE basicUnsafeMove #-}
-  basicUnsafeMove !dst !src@(BitMVec srcShift srcLen _)
+  basicUnsafeMove !dst src@(BitMVec srcShift srcLen _)
     | MV.basicOverlaps dst src = do
           -- Align shifts of src and srcCopy to speed up basicUnsafeCopy srcCopy src
       srcCopy <- MV.drop (modWordSize srcShift)
@@ -437,7 +412,7 @@
 #ifndef BITVEC_THREADSAFE
 
 -- | Flip the bit at the given position.
--- No bounds checks are performed.
+-- No bound checks are performed.
 -- Equivalent to 'flip' 'Data.Vector.Unboxed.Mutable.unsafeModify' 'Data.Bits.complement',
 -- but up to 2x faster.
 --
@@ -445,17 +420,26 @@
 -- either you modify it with 'id' (which is 'id' altogether)
 -- or with 'Data.Bits.complement' (which is 'unsafeFlipBit').
 --
--- >>> Data.Vector.Unboxed.modify (\v -> unsafeFlipBit v 1) (read "[1,1,1]")
--- [1,0,1]
+-- >>> :set -XOverloadedLists
+-- >>> Data.Vector.Unboxed.modify (`unsafeFlipBit` 2) [1,1,1,1]
+-- [1,1,0,1]
+--
+-- @since 1.0.0.0
 unsafeFlipBit :: PrimMonad m => U.MVector (PrimState m) Bit -> Int -> m ()
-unsafeFlipBit (BitMVec off _ arr) !i' = do
+unsafeFlipBit v i =
+  checkIndex Unsafe
+    i (MV.length v) $ basicFlipBit v i
+{-# INLINE unsafeFlipBit #-}
+
+basicFlipBit :: PrimMonad m => U.MVector (PrimState m) Bit -> Int -> m ()
+basicFlipBit (BitMVec off _ arr) !i' = do
   let i  = off + i'
       j  = divWordSize i
       k  = modWordSize i
       kk = 1 `unsafeShiftL` k :: Word
   word <- readByteArray arr j
   writeByteArray arr j (word `xor` kk)
-{-# INLINE unsafeFlipBit #-}
+{-# INLINE basicFlipBit #-}
 
 -- | Flip the bit at the given position.
 -- Equivalent to 'flip' 'Data.Vector.Unboxed.Mutable.modify' 'Data.Bits.complement',
@@ -465,17 +449,22 @@
 -- either you modify it with 'id' (which is 'id' altogether)
 -- or with 'Data.Bits.complement' (which is 'flipBit').
 --
--- >>> Data.Vector.Unboxed.modify (\v -> flipBit v 1) (read "[1,1,1]")
--- [1,0,1]
+-- >>> :set -XOverloadedLists
+-- >>> Data.Vector.Unboxed.modify (`flipBit` 2) [1,1,1,1]
+-- [1,1,0,1]
+--
+-- @since 1.0.0.0
 flipBit :: PrimMonad m => U.MVector (PrimState m) Bit -> Int -> m ()
 flipBit v i =
-  BOUNDS_CHECK(checkIndex) "flipBit" i (MV.length v) $ unsafeFlipBit v i
+  checkIndex Bounds
+    i (MV.length v) $
+      unsafeFlipBit v i
 {-# INLINE flipBit #-}
 
 #else
 
 -- | Flip the bit at the given position.
--- No bounds checks are performed.
+-- No bound checks are performed.
 -- Equivalent to 'flip' 'Data.Vector.Unboxed.Mutable.unsafeModify' 'Data.Bits.complement',
 -- but up to 33% faster and atomic.
 --
@@ -486,13 +475,19 @@
 -- >>> Data.Vector.Unboxed.modify (\v -> unsafeFlipBit v 1) (read "[1,1,1]")
 -- [1,0,1]
 unsafeFlipBit :: PrimMonad m => U.MVector (PrimState m) Bit -> Int -> m ()
-unsafeFlipBit (BitMVec off _ (MutableByteArray mba)) !i' = do
+unsafeFlipBit v i =
+  checkIndex Unsafe
+    i (MV.length v) $ basicFlipBit v i
+{-# INLINE unsafeFlipBit #-}
+
+basicFlipBit :: PrimMonad m => U.MVector (PrimState m) Bit -> Int -> m ()
+basicFlipBit (BitMVec off _ (MutableByteArray mba)) !i' = do
   let i       = off + i'
       !(I# j) = divWordSize i
       !(I# k) = 1 `unsafeShiftL` modWordSize i
   primitive $ \state ->
     let !(# state', _ #) = fetchXorIntArray# mba j k state in (# state', () #)
-{-# INLINE unsafeFlipBit #-}
+{-# INLINE basicFlipBit #-}
 
 -- | Flip the bit at the given position.
 -- Equivalent to 'flip' 'Data.Vector.Unboxed.Mutable.modify' 'Data.Bits.complement',
@@ -506,15 +501,15 @@
 -- [1,0,1]
 flipBit :: PrimMonad m => U.MVector (PrimState m) Bit -> Int -> m ()
 flipBit v i =
-  BOUNDS_CHECK(checkIndex) "flipBit" i (MV.length v) $ unsafeFlipBit v i
+  checkIndex Bounds
+    i (MV.length v) $ basicFlipBit v i
 {-# INLINE flipBit #-}
 
 #endif
 
 instance V.Vector U.Vector Bit where
-  basicUnsafeFreeze (BitMVec s n v) =
-    liftM (BitVec s n) (unsafeFreezeByteArray v)
-  basicUnsafeThaw (BitVec s n v) = liftM (BitMVec s n) (unsafeThawByteArray v)
+  basicUnsafeFreeze (BitMVec s n v) = BitVec s n <$> unsafeFreezeByteArray v
+  basicUnsafeThaw (BitVec s n v) = BitMVec s n <$> unsafeThawByteArray v
   basicLength (BitVec _ n _) = n
 
   basicUnsafeIndexM (BitVec off _ arr) !i' = do
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
@@ -1,9 +1,10 @@
-{-# LANGUAGE CPP              #-}
-
 {-# LANGUAGE BangPatterns        #-}
+{-# LANGUAGE CPP                 #-}
 {-# LANGUAGE FlexibleContexts    #-}
+{-# LANGUAGE MagicHash           #-}
 {-# LANGUAGE RankNTypes          #-}
 {-# LANGUAGE ScopedTypeVariables #-}
+{- HLINT ignore "Unused LANGUAGE pragma" -}
 
 #ifndef BITVEC_THREADSAFE
 module Data.Bit.Mutable
@@ -17,6 +18,7 @@
   , cloneToWords8M
 
   , zipInPlace
+  , mapInPlace
 
   , invertInPlace
   , selectBitsInPlace
@@ -25,6 +27,8 @@
   , reverseInPlace
   ) where
 
+#include "MachDeps.h"
+
 import Control.Monad
 import Control.Monad.Primitive
 import Control.Monad.ST
@@ -41,26 +45,35 @@
 import qualified Data.Vector.Unboxed.Mutable as MU
 import Data.Word
 
+#ifdef WORDS_BIGENDIAN
+import GHC.Exts
+#endif
+
 -- | Cast a vector of words to a vector of bits.
 -- Cf. 'Data.Bit.castFromWords'.
+--
+-- @since 1.0.0.0
 castFromWordsM :: MVector s Word -> MVector s Bit
 castFromWordsM (MU.MV_Word (P.MVector off len ws)) =
   BitMVec (mulWordSize off) (mulWordSize len) ws
 
 -- | Try to cast a vector of bits to a vector of words.
--- It succeeds if a vector of bits is aligned.
+-- It succeeds if the vector of bits is aligned.
 -- Use 'cloneToWordsM' otherwise.
 -- Cf. 'Data.Bit.castToWords'.
+--
+-- @since 1.0.0.0
 castToWordsM :: MVector s Bit -> Maybe (MVector s Word)
 castToWordsM (BitMVec s n ws)
-  | aligned s, aligned n = Just $ MU.MV_Word $ P.MVector (divWordSize s)
-                                                         (divWordSize n)
-                                                         ws
+  | aligned s, aligned n
+  = Just $ MU.MV_Word $ P.MVector (divWordSize s) (divWordSize n) ws
   | otherwise = Nothing
 
 -- | Clone a vector of bits to a new unboxed vector of words.
 -- If the bits don't completely fill the words, the last word will be zero-padded.
 -- Cf. 'Data.Bit.cloneToWords'.
+--
+-- @since 1.0.0.0
 cloneToWordsM
   :: PrimMonad m
   => MVector (PrimState m) Bit
@@ -72,7 +85,7 @@
   MU.unsafeCopy (MU.slice 0 lenBits w) v
   MU.set (MU.slice lenBits (mulWordSize lenWords - lenBits) w) (Bit False)
   pure $ MU.MV_Word $ P.MVector 0 lenWords arr
-{-# INLINE cloneToWordsM #-}
+{-# INLINABLE cloneToWordsM #-}
 
 -- | Clone a vector of bits to a new unboxed vector of 'Word8'.
 -- If the bits don't completely fill the words, the last 'Word8' will be zero-padded.
@@ -83,27 +96,46 @@
   -> m (MVector (PrimState m) Word8)
 cloneToWords8M v = do
   let lenBits  = MU.length v
-      lenWords = (lenBits + 7) `shiftR` 3
-  w@(BitMVec _ _ arr) <- MU.unsafeNew (lenWords `shiftL` 3)
-  MU.unsafeCopy (MU.slice 0 lenBits w) v
-  MU.set (MU.slice lenBits (lenWords `shiftL` 3 - lenBits) w) (Bit False)
-  pure $ MU.MV_Word8 $ P.MVector 0 lenWords arr
-{-# INLINE cloneToWords8M #-}
+      -- Take care about big-endian architectures: allocate full words!
+      actualLenBytes = (lenBits + 7) `shiftR` 3
+      roundedLenBytes = wordsToBytes (nWords lenBits)
+  ws@(BitMVec _ _ arr) <- MU.unsafeNew (roundedLenBytes `shiftL` 3)
+  MU.unsafeCopy (MU.slice 0 lenBits ws) v
+  MU.set (MU.slice lenBits (roundedLenBytes `shiftL` 3 - lenBits) ws) (Bit False)
 
--- | Zip two vectors with the given function.
--- rewriting contents of the second argument.
+#ifdef WORDS_BIGENDIAN
+  forM_ [0..nWords lenBits - 1] $ \i -> do
+    W# w <- readByteArray arr i
+    writeByteArray arr i (W# (byteSwap# w))
+#endif
+
+  pure $ MU.MV_Word8 $ P.MVector 0 actualLenBytes arr
+{-# INLINABLE cloneToWords8M #-}
+
+-- | Zip two vectors with the given function,
+-- rewriting the contents of the second argument.
 -- Cf. 'Data.Bit.zipBits'.
 --
+-- >>> :set -XOverloadedLists
 -- >>> import Data.Bits
--- >>> modify (zipInPlace (.&.) (read "[1,1,0]")) (read "[0,1,1]")
+-- >>> Data.Vector.Unboxed.modify (zipInPlace (.&.) [1,1,0]) [0,1,1]
 -- [0,1,0]
+-- >>> Data.Vector.Unboxed.modify (zipInPlace (\x y -> x .&. complement y) [1,1,0]) [0,1,1]
+-- [1,0,0]
 --
--- __Warning__: if the immutable vector is shorter than the mutable one,
--- it is a caller's responsibility to trim the result:
+-- __Warning__:
+-- If the immutable vector is longer than the mutable one,
+-- trailing bits will be ignored.
+-- If it's a mutable vector who is longer,
+-- trailing bits will be kept unchanged;
+-- it is caller's responsibility to trim the result:
 --
+-- >>> :set -XOverloadedLists
 -- >>> import Data.Bits
--- >>> modify (zipInPlace (.&.) (read "[1,1,0]")) (read "[0,1,1,1,1,1]")
+-- >>> Data.Vector.Unboxed.modify (zipInPlace (.&.) [1,1,0]) [0,1,1,1,1,1]
 -- [0,1,0,1,1,1] -- note trailing garbage
+--
+-- @since 1.0.0.0
 zipInPlace
   :: forall m.
      PrimMonad m
@@ -173,98 +205,142 @@
             writeByteArray ys i (f x y)
             loop (i + 1) accNew
 
-#if __GLASGOW_HASKELL__ >= 800
-{-# SPECIALIZE zipInPlace :: (forall a. Bits a => a -> a -> a) -> Vector Bit -> MVector s Bit -> ST s () #-}
-#endif
-{-# INLINE zipInPlace #-}
+{-# INLINABLE zipInPlace #-}
 
+-- | Apply a function to a mutable vector bitwise,
+-- rewriting its contents.
+-- Cf. 'Data.Bit.mapBits'.
+--
+-- >>> :set -XOverloadedLists
+-- >>> import Data.Bits
+-- >>> Data.Vector.Unboxed.modify (mapInPlace complement) [0,1,1]
+-- [1,0,0]
+--
+-- @since 1.1.0.0
+mapInPlace
+  :: PrimMonad m
+  => (forall a . Bits a => a -> a)
+  -> U.MVector (PrimState m) Bit
+  -> m ()
+mapInPlace f = case (unBit (f (Bit False)), unBit (f (Bit True))) of
+  (False, False) -> (`MU.set` Bit False)
+  (False, True)  -> const $ pure ()
+  (True, False)  -> invertInPlace
+  (True, True)   -> (`MU.set` Bit True)
+{-# INLINE mapInPlace #-}
+
 -- | Invert (flip) all bits in-place.
 --
--- >>> Data.Vector.Unboxed.modify invertInPlace (read "[0,1,0,1,0]")
+-- >>> :set -XOverloadedLists
+-- >>> Data.Vector.Unboxed.modify invertInPlace [0,1,0,1,0]
 -- [1,0,1,0,1]
+--
+-- @since 0.1
 invertInPlace :: PrimMonad m => U.MVector (PrimState m) Bit -> m ()
 invertInPlace xs = do
   let n = MU.length xs
   forM_ [0, wordSize .. n - 1] $ \i -> do
     x <- readWord xs i
     writeWord xs i (complement x)
-#if __GLASGOW_HASKELL__ >= 800
 {-# SPECIALIZE invertInPlace :: U.MVector s Bit -> ST s () #-}
-#endif
 
--- | Same as 'Data.Bit.selectBits', but deposit
--- selected bits in-place. Returns a number of selected bits.
--- It is caller's responsibility to trim the result to this number.
+-- | 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
+-- >>> runST $ do { vec <- U.unsafeThaw [1,1,0,0,1]; n <- selectBitsInPlace [0,1,0,1,1] vec; U.take n <$> U.unsafeFreeze vec }
+-- [1,0,1]
+--
+-- @since 0.1
 selectBitsInPlace
   :: PrimMonad m => U.Vector Bit -> U.MVector (PrimState m) Bit -> m Int
 selectBitsInPlace is xs = loop 0 0
  where
   !n = min (U.length is) (MU.length xs)
   loop !i !ct
-    | i >= n = return ct
+    | i >= n = pure ct
     | otherwise = do
       x <- readWord xs i
       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
--- excluded bits in-place. Returns a number of excluded bits.
--- It is caller's responsibility to trim the result to this number.
+-- | 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
+-- >>> runST $ do { vec <- U.unsafeThaw [1,1,0,0,1]; n <- excludeBitsInPlace [0,1,0,1,1] vec; U.take n <$> U.unsafeFreeze vec }
+-- [1,0]
+--
+-- @since 0.1
 excludeBitsInPlace
   :: PrimMonad m => U.Vector Bit -> U.MVector (PrimState m) Bit -> m Int
 excludeBitsInPlace is xs = loop 0 0
  where
   !n = min (U.length is) (MU.length xs)
   loop !i !ct
-    | i >= n = return ct
+    | i >= n = pure ct
     | otherwise = do
       x <- readWord xs i
       let !(nSet, x') =
             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.
 --
--- >>> Data.Vector.Unboxed.modify reverseInPlace (read "[1,1,0,1,0]")
+-- >>> :set -XOverloadedLists
+-- >>> Data.Vector.Unboxed.modify reverseInPlace [1,1,0,1,0]
 -- [0,1,0,1,1]
 --
--- Consider using @vector-rotcev@ package
+-- Consider using the [vector-rotcev](https://hackage.haskell.org/package/vector-rotcev) package
 -- to reverse vectors in O(1) time.
+--
+-- @since 0.1
 reverseInPlace :: PrimMonad m => U.MVector (PrimState m) Bit -> m ()
-reverseInPlace xs | len == 0  = pure ()
-                  | otherwise = loop 0
- where
-  len = MU.length xs
+reverseInPlace xs
+  | len == 0  = pure ()
+  | otherwise = loop 0
+  where
+    len = MU.length xs
 
-  loop !i
-    | i' <= j' = do
-      x <- readWord xs i
-      y <- readWord xs j'
+    loop !i
+      | i' <= j' = do
+        x <- readWord xs i
+        y <- readWord xs j'
 
-      writeWord xs i  (reverseWord y)
-      writeWord xs j' (reverseWord x)
+        writeWord xs i  (reverseWord y)
+        writeWord xs j' (reverseWord x)
 
-      loop i'
-    | i' < j = do
-      let w = (j - i) `shiftR` 1
-          k = j - w
-      x <- readWord xs i
-      y <- readWord xs k
+        loop i'
+      | i' < j = do
+        let w = (j - i) `shiftR` 1
+            k = j - w
+        x <- readWord xs i
+        y <- readWord xs k
 
-      writeWord xs i (meld w (reversePartialWord w y) x)
-      writeWord xs k (meld w (reversePartialWord w x) y)
+        writeWord xs i (meld w (reversePartialWord w y) x)
+        writeWord xs k (meld w (reversePartialWord w x) y)
 
-      loop i'
-    | otherwise = do
-      let w = j - i
-      x <- readWord xs i
-      writeWord xs i (meld w (reversePartialWord w x) x)
-   where
-    !j  = len - i
-    !i' = i + wordSize
-    !j' = j - wordSize
-#if __GLASGOW_HASKELL__ >= 800
+        loop i'
+      | otherwise = do
+        let w = j - i
+        x <- readWord xs i
+        writeWord xs i (meld w (reversePartialWord w x) x)
+     where
+      !j  = len - i
+      !i' = i + wordSize
+      !j' = j - wordSize
 {-# SPECIALIZE reverseInPlace :: U.MVector s Bit -> ST s () #-}
-#endif
diff --git a/src/Data/Bit/PdepPext.hs b/src/Data/Bit/PdepPext.hs
--- a/src/Data/Bit/PdepPext.hs
+++ b/src/Data/Bit/PdepPext.hs
@@ -1,12 +1,11 @@
 -- |
 -- Module:      Data.Bit.PdepPext
--- Copyright:   (c) 2020 Andrew Lelechenko
+-- Copyright:   (c) 2022 Andrew Lelechenko
 -- Licence:     BSD3
 --
 -- | Parallel bit deposit and extract instructions.
 -- https://en.wikipedia.org/wiki/Bit_Manipulation_Instruction_Sets#Parallel_bit_deposit_and_extract
 
-{-# LANGUAGE CPP          #-}
 {-# LANGUAGE MagicHash    #-}
 
 module Data.Bit.PdepPext
@@ -14,8 +13,6 @@
   , pext
   ) where
 
-#if MIN_VERSION_base(4,11,0)
-
 import GHC.Exts
 
 pdep :: Word -> Word -> Word
@@ -23,32 +20,3 @@
 
 pext :: Word -> Word -> Word
 pext (W# src#) (W# mask#) = W# (pext# src# mask#)
-
-#else
-
-import Data.Bits
-
-pdep :: Word -> Word -> Word
-pdep = go 0
-  where
-    go :: Word -> Word -> Word -> Word
-    go result _ 0 = result
-    go result src mask = go newResult newSrc newMask
-      where
-        lowest    = 1 `shiftL` countTrailingZeros mask
-        newResult = if src .&. 1 == 0 then result else result .|. lowest
-        newSrc    = src `shiftR` 1
-        newMask   = mask .&. complement lowest
-
-pext :: Word -> Word -> Word
-pext src mask = loop 0 0 0
-  where
-    loop i count acc
-      | i >= finiteBitSize (0 :: Word)
-      = acc
-      | testBit mask i
-      = loop (i + 1) (count + 1) (if testBit src i then setBit acc count else acc)
-      | otherwise
-      = loop (i + 1) count acc
-
-#endif
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/src/Data/Bit/Utils.hs b/src/Data/Bit/Utils.hs
--- a/src/Data/Bit/Utils.hs
+++ b/src/Data/Bit/Utils.hs
@@ -1,6 +1,4 @@
-{-# LANGUAGE BangPatterns               #-}
-{-# LANGUAGE CPP                        #-}
-{-# LANGUAGE MagicHash                  #-}
+{-# LANGUAGE MagicHash #-}
 
 module Data.Bit.Utils
   ( lgWordSize
@@ -25,29 +23,29 @@
   , toPrimVector
   ) where
 
-#include "MachDeps.h"
-
 import Data.Bits
 import qualified Data.Vector.Primitive as P
 import qualified Data.Vector.Unboxed as U
-#if __GLASGOW_HASKELL__ >= 810
+import qualified Data.Vector.Unboxed.Base as UB
 import GHC.Exts
-#endif
-import Unsafe.Coerce
 
 import Data.Bit.PdepPext
 
--- |The number of bits in a 'Word'.  A handy constant to have around when defining 'Word'-based bulk operations on bit vectors.
+-- | The number of bits in a 'Word'.  A handy constant to have around when defining 'Word'-based bulk operations on bit vectors.
 wordSize :: Int
 wordSize = finiteBitSize (0 :: Word)
 
-lgWordSize, wordSizeMask, wordSizeMaskC :: Int
+-- | The base 2 logarithm of 'wordSize'.
+lgWordSize :: Int
 lgWordSize = case wordSize of
   32 -> 5
   64 -> 6
   _  -> error "lgWordSize: unknown architecture"
 
+wordSizeMask :: Int
 wordSizeMask = wordSize - 1
+
+wordSizeMaskC :: Int
 wordSizeMaskC = complement wordSizeMask
 
 divWordSize :: Bits a => a -> a
@@ -87,11 +85,10 @@
 
 -- create a mask consisting of the lower n bits
 mask :: Int -> Word
-mask b = m
- where
-  m | b >= finiteBitSize m = complement 0
-    | b < 0                = 0
-    | otherwise            = bit b - 1
+mask b
+  | b >= wordSize = complement 0
+  | b < 0         = 0
+  | otherwise     = bit b - 1
 
 masked :: Int -> Word -> Word
 masked b x = x .&. mask b
@@ -101,35 +98,13 @@
 meld b lo hi = (lo .&. m) .|. (hi .&. complement m) where m = mask b
 {-# INLINE meld #-}
 
-#if __GLASGOW_HASKELL__ >= 810
 reverseWord :: Word -> Word
 reverseWord (W# w#) = W# (bitReverse# w#)
-#elif WORD_SIZE_IN_BITS == 64
-reverseWord :: Word -> Word
-reverseWord x0 = x6
- where
-  x1 = ((x0 .&. 0x5555555555555555) `shiftL`  1) .|. ((x0 .&. 0xAAAAAAAAAAAAAAAA) `shiftR`  1)
-  x2 = ((x1 .&. 0x3333333333333333) `shiftL`  2) .|. ((x1 .&. 0xCCCCCCCCCCCCCCCC) `shiftR`  2)
-  x3 = ((x2 .&. 0x0F0F0F0F0F0F0F0F) `shiftL`  4) .|. ((x2 .&. 0xF0F0F0F0F0F0F0F0) `shiftR`  4)
-  x4 = ((x3 .&. 0x00FF00FF00FF00FF) `shiftL`  8) .|. ((x3 .&. 0xFF00FF00FF00FF00) `shiftR`  8)
-  x5 = ((x4 .&. 0x0000FFFF0000FFFF) `shiftL` 16) .|. ((x4 .&. 0xFFFF0000FFFF0000) `shiftR` 16)
-  x6 = ((x5 .&. 0x00000000FFFFFFFF) `shiftL` 32) .|. ((x5 .&. 0xFFFFFFFF00000000) `shiftR` 32)
-#elif WORD_SIZE_IN_BITS == 32
-reverseWord :: Word -> Word
-reverseWord x0 = x5
- where
-  x1 = ((x0 .&. 0x55555555) `shiftL`  1) .|. ((x0 .&. 0xAAAAAAAA) `shiftR`  1)
-  x2 = ((x1 .&. 0x33333333) `shiftL`  2) .|. ((x1 .&. 0xCCCCCCCC) `shiftR`  2)
-  x3 = ((x2 .&. 0x0F0F0F0F) `shiftL`  4) .|. ((x2 .&. 0xF0F0F0F0) `shiftR`  4)
-  x4 = ((x3 .&. 0x00FF00FF) `shiftL`  8) .|. ((x3 .&. 0xFF00FF00) `shiftR`  8)
-  x5 = ((x4 .&. 0x0000FFFF) `shiftL` 16) .|. ((x4 .&. 0xFFFF0000) `shiftR` 16)
-#else
-#error unsupported WORD_SIZE_IN_BITS config
-#endif
 
 reversePartialWord :: Int -> Word -> Word
-reversePartialWord n w | n >= wordSize = reverseWord w
-                       | otherwise     = reverseWord w `shiftR` (wordSize - n)
+reversePartialWord n w
+  | n >= wordSize = reverseWord w
+  | otherwise     = reverseWord w `shiftR` (wordSize - n)
 
 ffs :: Word -> Maybe Int
 ffs 0 = Nothing
@@ -140,18 +115,22 @@
 selectWord msk src = (popCount msk, pext src msk)
 {-# INLINE selectWord #-}
 
-#if WORD_SIZE_IN_BITS == 64
-
 -- | Insert 0 between each consecutive bits of an input.
 -- xyzw --> (x0y0, z0w0)
 sparseBits :: Word -> (Word, Word)
-sparseBits w = (x, y)
+sparseBits = case wordSize of
+  32 -> sparseBits32
+  64 -> sparseBits64
+  _  -> error "sparseBits: unknown architecture"
+
+sparseBits64 :: Word -> (Word, Word)
+sparseBits64 w = (x, y)
   where
-    x = sparseBitsInternal (w .&. loMask 32)
-    y = sparseBitsInternal (w `shiftR` 32)
+    x = sparseBitsInternal64 (w .&. loMask 32)
+    y = sparseBitsInternal64 (w `shiftR` 32)
 
-sparseBitsInternal :: Word -> Word
-sparseBitsInternal x = x4
+sparseBitsInternal64 :: Word -> Word
+sparseBitsInternal64 x = x4
   where
     t  = (x  `xor` (x  `shiftR` 16)) .&. 0x00000000ffff0000
     x0 = x  `xor` (t  `xor` (t  `shiftL` 16));
@@ -165,18 +144,14 @@
     t3 = (x3 `xor` (x3 `shiftR` 1)) .&. 0x2222222222222222;
     x4 = x3 `xor` (t3 `xor` (t3 `shiftL` 1));
 
-#elif WORD_SIZE_IN_BITS == 32
-
--- | Insert 0 between each consecutive bits of an input.
--- xyzw --> (x0y0, z0w0)
-sparseBits :: Word -> (Word, Word)
-sparseBits w = (x, y)
+sparseBits32 :: Word -> (Word, Word)
+sparseBits32 w = (x, y)
   where
-    x = sparseBitsInternal (w .&. loMask 16)
-    y = sparseBitsInternal (w `shiftR` 16)
+    x = sparseBitsInternal32 (w .&. loMask 16)
+    y = sparseBitsInternal32 (w `shiftR` 16)
 
-sparseBitsInternal :: Word -> Word
-sparseBitsInternal x0 = x4
+sparseBitsInternal32 :: Word -> Word
+sparseBitsInternal32 x0 = x4
   where
     t0 = (x0 `xor` (x0 `shiftR` 8)) .&. 0x0000ff00;
     x1 = x0 `xor` (t0 `xor` (t0 `shiftL` 8));
@@ -187,10 +162,6 @@
     t3 = (x3 `xor` (x3 `shiftR` 1)) .&. 0x22222222;
     x4 = x3 `xor` (t3 `xor` (t3 `shiftL` 1));
 
-#else
-#error unsupported WORD_SIZE_IN_BITS config
-#endif
-
 loMask :: Int -> Word
 loMask n = 1 `unsafeShiftL` n - 1
 {-# INLINE loMask #-}
@@ -200,9 +171,9 @@
 {-# INLINE hiMask #-}
 
 fromPrimVector :: P.Vector Word -> U.Vector Word
-fromPrimVector = unsafeCoerce
+fromPrimVector = UB.V_Word
 {-# INLINE fromPrimVector #-}
 
 toPrimVector :: U.Vector Word -> P.Vector Word
-toPrimVector = unsafeCoerce
+toPrimVector (UB.V_Word ws) = ws
 {-# INLINE toPrimVector #-}
diff --git a/test/Main.hs b/test/Main.hs
--- a/test/Main.hs
+++ b/test/Main.hs
@@ -1,20 +1,20 @@
 {-# LANGUAGE CPP       #-}
-{-# LANGUAGE MagicHash #-}
 
 module Main where
 
+import Control.Exception
 import Data.Bit
-import Data.Bits
-import Data.Proxy
-import qualified Data.Vector.Unboxed as U
-import GHC.Exts
-import GHC.Integer.Logarithms
-import Test.QuickCheck.Classes
 import Test.Tasty
 import Test.Tasty.QuickCheck
 
+#ifdef MIN_VERSION_quickcheck_classes_base
+import Data.Proxy
+import Test.QuickCheck.Classes.Base
 import Support
+#endif
+
 import Tests.Conc (concTests)
+import Tests.F2Poly (f2polyTests)
 import Tests.MVector (mvectorTests)
 import qualified Tests.MVectorTS as TS (mvectorTests)
 import Tests.SetOps (setOpTests)
@@ -36,6 +36,7 @@
 lawsTests :: TestTree
 lawsTests = adjustOption (const $ QuickCheckTests 100)
   $ testGroup "Bit"
+#ifdef MIN_VERSION_quickcheck_classes_base
   $ map lawsToTest
   [ bitsLaws        (Proxy :: Proxy Bit)
   , eqLaws          (Proxy :: Proxy Bit)
@@ -43,78 +44,17 @@
   , boundedEnumLaws (Proxy :: Proxy Bit)
   , showLaws        (Proxy :: Proxy Bit)
   , showReadLaws    (Proxy :: Proxy Bit)
-#if MIN_VERSION_quickcheck_classes(0,6,3)
   , numLaws         (Proxy :: Proxy Bit)
-#endif
   , integralLaws    (Proxy :: Proxy Bit)
-  ]
-
-f2polyTests :: TestTree
-f2polyTests = testGroup "F2Poly"
-  [ testProperty "Addition"       prop_f2polyAdd
-  , testProperty "Multiplication" prop_f2polyMul
-  , testProperty "Square" prop_f2polySqr
-  , tenTimesLess $ testProperty "Multiplication long" prop_f2polyMulLong
-  , tenTimesLess $ testProperty "Square long" prop_f2polySqrLong
-  , testProperty "Remainder"      prop_f2polyRem
-  , testProperty "GCD"            prop_f2polyGCD
-  , tenTimesLess $ lawsToTest $
-    showLaws (Proxy :: Proxy F2Poly)
-#if MIN_VERSION_quickcheck_classes(0,6,3)
-  , lawsToTest $
-    numLaws (Proxy :: Proxy F2Poly)
+  ] ++
 #endif
-  , lawsToTest $
-    integralLaws (Proxy :: Proxy F2Poly)
+  [ testProperty "divideByZero" prop_bitDivideByZero
+  , testProperty "toRational"   prop_bitToRational
   ]
 
-prop_f2polyAdd :: F2Poly -> F2Poly -> Property
-prop_f2polyAdd x y = x + y === fromInteger (toInteger x `xor` toInteger y)
-
-prop_f2polyMul :: F2Poly -> F2Poly -> Property
-prop_f2polyMul x y = x * y === fromInteger (toInteger x `binMul` toInteger y)
-
-prop_f2polySqr :: F2Poly -> Property
-prop_f2polySqr x = x * x === fromInteger (toInteger x `binMul` toInteger x)
-
-prop_f2polyMulLong :: U.Vector Word -> U.Vector Word -> Property
-prop_f2polyMulLong xs ys = x * y === fromInteger (toInteger x `binMul` toInteger y)
-  where
-    x = toF2Poly $ castFromWords xs
-    y = toF2Poly $ castFromWords ys
-
-prop_f2polySqrLong :: U.Vector Word -> Property
-prop_f2polySqrLong xs = x * x === fromInteger (toInteger x `binMul` toInteger x)
-  where
-    x = toF2Poly $ castFromWords xs
-
-prop_f2polyRem :: F2Poly -> F2Poly -> Property
-prop_f2polyRem x y = y /= 0 ==> x `rem` y === fromInteger (toInteger x `binRem` toInteger y)
-
--- For polynomials @x@ and @y@, @gcdExt@ computes their unique greatest common
--- divisor @g@ and the unique coefficient polynomial @s@ satisfying @xs + yt = g@.
---
--- Thus it is sufficient to check @gcd == fst . gcdExt@ and @xs == g (mod y)@,
--- except if @y@ divides @x@, then @gcdExt x y@ is @(y, 0)@ and @xs `rem` y@ is zero,
--- so that it is then necessary to check @xs `rem` y == g `rem` y == 0@.
-prop_f2polyGCD :: F2Poly -> F2Poly -> Property
-prop_f2polyGCD x y = g === x `gcd` y .&&. (y /= 0 ==> (x * s) `rem` y === g `rem` y)
-  where
-    (g, s) = x `gcdExt` y
-
-binMul :: Integer -> Integer -> Integer
-binMul = go 0
-  where
-    go :: Integer -> Integer -> Integer -> Integer
-    go acc _ 0 = acc
-    go acc x y = go (if odd y then acc `xor` x else acc) (x `shiftL` 1) (y `shiftR` 1)
-
-binRem :: Integer -> Integer -> Integer
-binRem x y = go x
-  where
-    binLog n = I# (integerLog2# n)
-    ly = binLog y
+prop_bitToRational :: Bit -> Property
+prop_bitToRational x = fromRational (toRational x) === x
 
-    go z = if lz < ly then z else go (z `xor` (y `shiftL` (lz - ly)))
-      where
-        lz = binLog z
+prop_bitDivideByZero :: Bit -> Property
+prop_bitDivideByZero x =
+  ioProperty ((=== Left DivideByZero) <$> try (evaluate (x / 0)))
diff --git a/test/Support.hs b/test/Support.hs
--- a/test/Support.hs
+++ b/test/Support.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE CPP                 #-}
 {-# LANGUAGE FlexibleContexts    #-}
 {-# LANGUAGE FlexibleInstances   #-}
 {-# LANGUAGE RankNTypes          #-}
@@ -15,10 +16,16 @@
 import qualified Data.Vector.Generic.Mutable as M
 import qualified Data.Vector.Generic.New as N
 import qualified Data.Vector.Unboxed as U
-import Test.QuickCheck.Classes
 import Test.Tasty
 import Test.Tasty.QuickCheck
 
+#ifdef MIN_VERSION_quickcheck_classes_base
+import Test.QuickCheck.Classes.Base
+#endif
+#ifdef MIN_VERSION_quickcheck_classes
+import qualified Test.QuickCheck.Classes as QCC
+#endif
+
 instance Arbitrary Bit where
   arbitrary = Bit <$> arbitrary
   shrink    = fmap Bit . shrink . unBit
@@ -40,18 +47,35 @@
   function f = functionMap TS.unBit TS.Bit f
 
 instance (Arbitrary a, U.Unbox a) => Arbitrary (U.Vector a) where
-  arbitrary = (\v -> runST (N.run (v :: N.New U.Vector a) >>= U.freeze)) <$> arbitrary
-  shrink v = let len = U.length v in
-    [ U.slice s l v
-    | s <- [0 .. len - 1]
-    , l <- [0 .. len - s]
-    , (s, l) /= (0, len)
+  arbitrary = frequency
+    [ (10, U.fromList <$> arbitrary)
+    , (2 , U.drop <$> arbitrary <*> arbitrary)
+    , (2 , U.take <$> arbitrary <*> arbitrary)
+    , (2 , slice <$> arbitrary <*> arbitrary <*> arbitrary)
     ]
+   where
+    slice s n v = let (s', n') = trimSlice s n (U.length v) in U.slice s' n' v
+  shrink v = let len = U.length v in
+    [ U.take (len - s) v | s <- [1 .. len] ] ++
+    [ U.drop s         v | s <- [1 .. len] ] ++
+    [ v U.// [(i, x)] | i <- [0 .. len - 1], x <- shrink (v U.! i) ]
 
+instance {-# OVERLAPPING #-} Arbitrary (Large (U.Vector Bit)) where
+  arbitrary = Large . castFromWords <$> arbitrary
+  shrink (Large v) = Large <$> shrink v
+
+instance {-# OVERLAPPING #-} Arbitrary (Large (U.Vector TS.Bit)) where
+  arbitrary = Large . TS.castFromWords <$> arbitrary
+  shrink (Large v) = Large <$> shrink v
+
 instance Arbitrary F2Poly where
   arbitrary = toF2Poly <$> arbitrary
   shrink v = toF2Poly <$> shrink (unF2Poly v)
 
+instance {-# OVERLAPPING #-} Arbitrary (Large F2Poly) where
+  arbitrary = Large . toF2Poly . castFromWords <$> arbitrary
+  shrink (Large v) = Large . toF2Poly <$> shrink (unF2Poly v)
+
 instance (Show (v a), V.Vector v a) => Show (N.New v a) where
   showsPrec p = showsPrec p . V.new
 
@@ -62,19 +86,16 @@
 instance (V.Vector v a, Arbitrary a) => Arbitrary (N.New v a) where
   arbitrary = frequency
     [ (10, newFromList <$> arbitrary)
-    , (1 , N.drop <$> arbitrary <*> arbitrary)
-    , (1 , N.take <$> arbitrary <*> arbitrary)
-    , (1 , slice <$> arbitrary <*> arbitrary <*> arbitrary)
+    , (2 , N.drop <$> arbitrary <*> arbitrary)
+    , (2 , N.take <$> arbitrary <*> arbitrary)
+    , (2 , slice <$> arbitrary <*> arbitrary <*> arbitrary)
     ]
    where
     slice s n = N.apply
       $ \v -> let (s', n') = trimSlice s n (M.length v) in M.slice s' n' v
   shrink v =
-    [ N.slice s l v
-    | s <- [0 .. len - 1]
-    , l <- [0 .. len - s]
-    , (s, l) /= (0, len)
-    ]
+    [ N.take s v | s <- [0 .. len - 1] ] ++
+    [ N.drop s v | s <- [1 .. len] ]
     where len = runST (M.length <$> N.run v)
 
 trimSlice :: Integral a => a -> a -> a -> (a, a)
@@ -138,6 +159,18 @@
 tenTimesLess = adjustOption $
   \(QuickCheckTests n) -> QuickCheckTests (max 100 (n `div` 10))
 
+twoTimesMore :: TestTree -> TestTree
+twoTimesMore = adjustOption $
+  \(QuickCheckTests n) -> QuickCheckTests (n * 2)
+
+#ifdef MIN_VERSION_quickcheck_classes_base
 lawsToTest :: Laws -> TestTree
 lawsToTest (Laws name props) =
   testGroup name $ map (uncurry testProperty) props
+#endif
+
+#ifdef MIN_VERSION_quickcheck_classes
+lawsToTest' :: QCC.Laws -> TestTree
+lawsToTest' (QCC.Laws name props) =
+  testGroup name $ map (uncurry testProperty) props
+#endif
diff --git a/test/Tests/Conc.hs b/test/Tests/Conc.hs
--- a/test/Tests/Conc.hs
+++ b/test/Tests/Conc.hs
@@ -1,4 +1,6 @@
-module Tests.Conc where
+module Tests.Conc
+  ( concTests
+  ) where
 
 import Control.Concurrent
 import Control.Monad
@@ -8,13 +10,13 @@
 import qualified Data.Vector.Generic.Mutable as M
 import qualified Data.Vector.Unboxed as U
 import Test.Tasty
-import Test.Tasty.HUnit
+import Test.Tasty.QuickCheck
 
 concTests :: TestTree
 concTests = testGroup "Concurrency"
-  [ testCase "invertInPlace"  case_conc_invert
-  , testCase "reverseInPlace" case_conc_reverse
-  , testCase "zipInPlace"     case_conc_zip
+  [ testProperty "invertInPlace"  case_conc_invert
+  , testProperty "reverseInPlace" case_conc_reverse
+  , testProperty "zipInPlace"     case_conc_zip
   ]
 
 runConcurrently :: IO () -> IO () -> IO ()
@@ -26,8 +28,8 @@
   action2
   takeMVar m
 
-case_conc_invert :: IO ()
-case_conc_invert = replicateM_ 1000 $ do
+case_conc_invert :: Property
+case_conc_invert = once $ ioProperty $ replicateM_ 1000 $ do
   let len  = 64
       len' = 37
   vec <- M.replicate len (Bit True)
@@ -36,10 +38,10 @@
     (replicateM_ 1000 $ invertInPlace (M.take len' vec))
     (replicateM_ 1000 $ invertInPlace (M.drop len' vec))
   wec <- V.unsafeFreeze vec
-  assertEqual "should be equal" ref wec
+  pure $ ref === wec
 
-case_conc_reverse :: IO ()
-case_conc_reverse = replicateM_ 1000 $ do
+case_conc_reverse :: Property
+case_conc_reverse = once $ ioProperty $ replicateM_ 1000 $ do
   let len  = 128
       len' = 66
   vec <- M.new len
@@ -49,10 +51,10 @@
     (replicateM_ 1000 $ reverseInPlace (M.take len' vec))
     (replicateM_ 1000 $ reverseInPlace (M.drop len' vec))
   wec <- V.unsafeFreeze vec
-  assertEqual "should be equal" ref wec
+  pure $ ref === wec
 
-case_conc_zip :: IO ()
-case_conc_zip = replicateM_ 1000 $ do
+case_conc_zip :: Property
+case_conc_zip = once $ ioProperty $ replicateM_ 1000 $ do
   let len  = 128
       len' = 37
   vec <- M.replicate len (Bit True)
@@ -61,4 +63,4 @@
     (replicateM_ 1001 $ zipInPlace (const complement) ref (M.take len' vec))
     (replicateM_ 1001 $ zipInPlace (const complement) ref (M.drop len' vec))
   wec <- V.unsafeFreeze vec
-  assertEqual "should be equal" ref wec
+  pure $ ref === wec
diff --git a/test/Tests/F2Poly.hs b/test/Tests/F2Poly.hs
new file mode 100644
--- /dev/null
+++ b/test/Tests/F2Poly.hs
@@ -0,0 +1,118 @@
+{-# LANGUAGE CPP       #-}
+{-# LANGUAGE MagicHash #-}
+
+module Tests.F2Poly
+  ( f2polyTests
+  ) where
+
+import Control.Exception
+import Data.Bit
+import Data.Bits
+import Data.Ratio
+import GHC.Exts
+import GHC.Num.Integer
+import Test.Tasty
+import Test.Tasty.QuickCheck
+
+#ifdef MIN_VERSION_quickcheck_classes_base
+import Data.Proxy
+import Test.QuickCheck.Classes.Base
+#endif
+
+import Support
+
+f2polyTests :: TestTree
+f2polyTests = testGroup "F2Poly"
+  [ testProperty "Addition"            prop_f2polyAdd
+  , testProperty "Multiplication"      prop_f2polyMul
+  , testProperty "Square"              prop_f2polySqr
+  , tenTimesLess
+  $ testProperty "Multiplication long" prop_f2polyMulLong
+  , testProperty "Multiplication 1"    prop_f2polyMul1
+  , tenTimesLess
+  $ testProperty "Square long"         prop_f2polySqrLong
+  , testProperty "Remainder"           prop_f2polyRem
+  , testProperty "GCD"                 prop_f2polyGCD
+  , testProperty "Enum" $
+    \n -> let x = toEnum n in toEnum (fromEnum x) === (x :: F2Poly)
+#ifdef MIN_VERSION_quickcheck_classes_base
+  , tenTimesLess $ lawsToTest $
+    showLaws (Proxy :: Proxy F2Poly)
+  , lawsToTest $
+    numLaws (Proxy :: Proxy F2Poly)
+  , lawsToTest $
+    integralLaws (Proxy :: Proxy F2Poly)
+#endif
+  , testProperty "fromNegative" prop_f2polyFromNegative
+  , testProperty "divideByZero" prop_f2polyDivideByZero
+  , testProperty "toRational" prop_f2polyToRational
+  , testProperty "signum" $ \x -> x + signum x === (x + 1 :: F2Poly)
+  ]
+
+prop_f2polyAdd :: F2Poly -> F2Poly -> Property
+prop_f2polyAdd x y = x + y === fromInteger (toInteger x `xor` toInteger y)
+
+prop_f2polyMul :: F2Poly -> F2Poly -> Property
+prop_f2polyMul x y = x * y === fromInteger (toInteger x `binMul` toInteger y)
+
+prop_f2polySqr :: F2Poly -> Property
+prop_f2polySqr x = x * x === fromInteger (toInteger x `binMul` toInteger x)
+
+prop_f2polyMulLong :: Large F2Poly -> Large F2Poly -> Property
+prop_f2polyMulLong (Large x) (Large y) = prop_f2polyMul x y
+
+prop_f2polyMul1 :: Property
+prop_f2polyMul1 = prop_f2polyMul x y
+  where
+    x = fromInteger (1 `shiftL` 4358)
+    y = fromInteger (1 `shiftL` 4932 + 1 `shiftL` 2116)
+
+prop_f2polySqrLong :: Large F2Poly -> Property
+prop_f2polySqrLong (Large x) = prop_f2polySqr x
+
+prop_f2polyRem :: F2Poly -> F2Poly -> Property
+prop_f2polyRem x y = y /= 0 ==> x `rem` y === fromInteger (toInteger x `binRem` toInteger y)
+
+-- For polynomials @x@ and @y@, @gcdExt@ computes their unique greatest common
+-- divisor @g@ and the unique coefficient polynomial @s@ satisfying @xs + yt = g@.
+--
+-- Thus it is sufficient to check @gcd == fst . gcdExt@ and @xs == g (mod y)@,
+-- except if @y@ divides @x@, then @gcdExt x y@ is @(y, 0)@ and @xs `rem` y@ is zero,
+-- so that it is then necessary to check @xs `rem` y == g `rem` y == 0@.
+prop_f2polyGCD :: F2Poly -> F2Poly -> Property
+prop_f2polyGCD x y = g === x `gcd` y .&&. (y /= 0 ==> (x * s) `rem` y === g `rem` y)
+  where
+    (g, s) = x `gcdExt` y
+
+binMul :: Integer -> Integer -> Integer
+binMul = go 0
+  where
+    go :: Integer -> Integer -> Integer -> Integer
+    go acc _ 0 = acc
+    go acc x y = go (if odd y then acc `xor` x else acc) (x `shiftL` 1) (y `shiftR` 1)
+
+binRem :: Integer -> Integer -> Integer
+binRem x y = go x
+  where
+    binLog n = I# (word2Int# (integerLog2# n))
+    ly = binLog y
+
+    go 0 = 0
+    go z = if lz < ly then z else go (z `xor` (y `shiftL` (lz - ly)))
+      where
+        lz = binLog z
+
+prop_f2polyFromNegative :: Large Int -> Property
+prop_f2polyFromNegative (Large m) =
+  ioProperty ((=== Left Underflow) <$> try (evaluate (fromInteger neg :: F2Poly)))
+  where
+    neg = negate (1 + toInteger m * toInteger m)
+
+prop_f2polyToRational :: F2Poly -> Property
+prop_f2polyToRational x = denominator y === 1 .&&. fromInteger (numerator y) === x
+  where
+    y = toRational x
+
+prop_f2polyDivideByZero :: F2Poly -> Property
+prop_f2polyDivideByZero x =
+  ioProperty ((=== Left DivideByZero) <$> try (evaluate (x `quot` 0)))
diff --git a/test/Tests/MVector.hs b/test/Tests/MVector.hs
--- a/test/Tests/MVector.hs
+++ b/test/Tests/MVector.hs
@@ -1,13 +1,17 @@
 {-# LANGUAGE CPP #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 
+
 #ifndef BITVEC_THREADSAFE
-module Tests.MVector where
+module Tests.MVector (mvectorTests) where
 #else
-module Tests.MVectorTS where
+module Tests.MVectorTS (mvectorTests) where
 #endif
 
 import Support
 
+import Control.Exception
+import Control.Monad
 import Control.Monad.ST
 #ifndef BITVEC_THREADSAFE
 import Data.Bit
@@ -15,17 +19,19 @@
 import Data.Bit.ThreadSafe
 #endif
 import Data.Bits
-import Data.Proxy
 import qualified Data.Vector.Generic as V
-import qualified Data.Vector.Generic.Mutable as M (basicInitialize, basicSet)
+import qualified Data.Vector.Generic.Mutable as MG
 import qualified Data.Vector.Generic.New as N
-import qualified Data.Vector.Unboxed as B
+import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as M
-import Test.QuickCheck.Classes
 import Test.Tasty
-import Test.Tasty.HUnit
 import Test.Tasty.QuickCheck
 
+#ifdef MIN_VERSION_quickcheck_classes
+import Data.Proxy
+import Test.QuickCheck.Classes
+#endif
+
 mvectorTests :: TestTree
 mvectorTests = testGroup "Data.Vector.Unboxed.Mutable.Bit"
   [ testGroup "Data.Vector.Unboxed.Mutable functions"
@@ -35,160 +41,171 @@
     ]
   , testGroup "Read/write Words"
     [ tenTimesLess $
-      testProperty "cloneFromWords" prop_cloneFromWords_def
+      testProperty "castFromWords"  prop_castFromWords_def
     , testProperty "cloneToWords"   prop_cloneToWords_def
+    , tenTimesLess $
+      testProperty "castToWords_1"  prop_castToWords_1
+    , tenTimesLess $
+      testProperty "castToWords_2"  prop_castToWords_2
     ]
-  , lawsToTest $ muvectorLaws (Proxy :: Proxy Bit)
-  , testCase "basicInitialize 1" case_write_init_read1
-  , testCase "basicInitialize 2" case_write_init_read2
-  , testCase "basicInitialize 3" case_write_init_read3
-  , testCase "basicInitialize 4" case_write_init_read4
-  , testCase "basicSet 1"        case_write_set_read1
-  , testCase "basicSet 2"        case_write_set_read2
-  , testCase "basicSet 3"        case_write_set_read3
-  , testCase "basicSet 4"        case_write_set_read4
-  , testCase "basicSet 5"        case_set_read1
-  , testCase "basicSet 6"        case_set_read2
-  , testCase "basicSet 7"        case_set_read3
-  , testCase "basicSet 8"        case_set_read4
-  , testCase "basicUnsafeCopy1"  case_write_copy_read1
-  , testCase "basicUnsafeCopy2"  case_write_copy_read2
-  , testCase "basicUnsafeCopy3"  case_write_copy_read3
-  , testCase "basicUnsafeCopy4"  case_write_copy_read4
-  , testCase "basicUnsafeCopy5"  case_write_copy_read5
+#ifdef MIN_VERSION_quickcheck_classes
+  , lawsToTest' $ muvectorLaws (Proxy :: Proxy Bit)
+#endif
+  , testProperty "basicInitialize 1" case_write_init_read1
+  , testProperty "basicInitialize 2" case_write_init_read2
+  , testProperty "basicInitialize 3" case_write_init_read3
+  , testProperty "basicInitialize 4" case_write_init_read4
+  , testProperty "basicSet 1"        case_write_set_read1
+  , testProperty "basicSet 2"        case_write_set_read2
+  , testProperty "basicSet 3"        case_write_set_read3
+  , testProperty "basicSet 4"        case_write_set_read4
+  , testProperty "basicSet 5"        case_set_read1
+  , testProperty "basicSet 6"        case_set_read2
+  , testProperty "basicSet 7"        case_set_read3
+  , testProperty "basicSet 8"        case_set_read4
+  , testProperty "basicUnsafeCopy1"  case_write_copy_read1
+  , testProperty "basicUnsafeCopy2"  case_write_copy_read2
+  , testProperty "basicUnsafeCopy3"  case_write_copy_read3
+  , testProperty "basicUnsafeCopy4"  case_write_copy_read4
+  , testProperty "basicUnsafeCopy5"  case_write_copy_read5
   , tenTimesLess $
     testProperty "flipBit" prop_flipBit
+  , testProperty "new negative"       prop_new_neg
+  , testProperty "replicate negative" prop_replicate_neg
+  , testProperty "move preserves data around" prop_move_around
+  , adjustOption (\(QuickCheckMaxRatio n) -> QuickCheckMaxRatio (max 1000 n)) $
+    testProperty "copy preserves data around" prop_copy_around
   ]
 
-prop_flipBit :: B.Vector Bit -> NonNegative Int -> Property
-prop_flipBit xs (NonNegative k) = B.length xs > 0 ==> ys === ys'
+prop_flipBit :: U.Vector Bit -> NonNegative Int -> Property
+prop_flipBit xs (NonNegative k) = U.length xs > 0 ==> ys === ys'
   where
-    k'  = k `mod` B.length xs
-    ys  = B.modify (\v -> M.modify v complement k') xs
-    ys' = B.modify (\v -> flipBit v k') xs
+    k'  = k `mod` U.length xs
+    ys  = U.modify (\v -> M.modify v complement k') xs
+    ys' = U.modify (`flipBit` k') xs
 
-case_write_init_read1 :: IO ()
-case_write_init_read1 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_write_init_read1 :: Property
+case_write_init_read1 = (=== Bit True) $ runST $ do
   arr <- M.new 2
   M.write arr 0 (Bit True)
-  M.basicInitialize (M.slice 1 1 arr)
+  MG.basicInitialize (M.slice 1 1 arr)
   M.read arr 0
 
-case_write_init_read2 :: IO ()
-case_write_init_read2 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_write_init_read2 :: Property
+case_write_init_read2 = (=== Bit True) $ runST $ do
   arr <- M.new 2
   M.write arr 1 (Bit True)
-  M.basicInitialize (M.slice 0 1 arr)
+  MG.basicInitialize (M.slice 0 1 arr)
   M.read arr 1
 
-case_write_init_read3 :: IO ()
+case_write_init_read3 :: Property
 case_write_init_read3 =
-  assertEqual "should be equal" (Bit True, Bit True) $ runST $ do
+  (=== (Bit True, Bit True)) $ runST $ do
     arr <- M.new 2
     M.write arr 0 (Bit True)
     M.write arr 1 (Bit True)
-    M.basicInitialize (M.slice 1 0 arr)
+    MG.basicInitialize (M.slice 1 0 arr)
     (,) <$> M.read arr 0 <*> M.read arr 1
 
-case_write_init_read4 :: IO ()
+case_write_init_read4 :: Property
 case_write_init_read4 =
-  assertEqual "should be equal" (Bit True, Bit True) $ runST $ do
+  (=== (Bit True, Bit True)) $ runST $ do
     arr <- M.new 3
     M.write arr 0 (Bit True)
     M.write arr 2 (Bit True)
-    M.basicInitialize (M.slice 1 1 arr)
+    MG.basicInitialize (M.slice 1 1 arr)
     (,) <$> M.read arr 0 <*> M.read arr 2
 
-case_write_set_read1 :: IO ()
-case_write_set_read1 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_write_set_read1 :: Property
+case_write_set_read1 = (=== Bit True) $ runST $ do
   arr <- M.new 2
   M.write arr 0 (Bit True)
-  M.basicSet (M.slice 1 1 arr) (Bit False)
+  MG.basicSet (M.slice 1 1 arr) (Bit False)
   M.read arr 0
 
-case_write_set_read2 :: IO ()
-case_write_set_read2 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_write_set_read2 :: Property
+case_write_set_read2 = (=== Bit True) $ runST $ do
   arr <- M.new 2
   M.write arr 1 (Bit True)
-  M.basicSet (M.slice 0 1 arr) (Bit False)
+  MG.basicSet (M.slice 0 1 arr) (Bit False)
   M.read arr 1
 
-case_write_set_read3 :: IO ()
+case_write_set_read3 :: Property
 case_write_set_read3 =
-  assertEqual "should be equal" (Bit True, Bit True) $ runST $ do
+  (=== (Bit True, Bit True)) $ runST $ do
     arr <- M.new 2
     M.write arr 0 (Bit True)
     M.write arr 1 (Bit True)
-    M.basicSet (M.slice 1 0 arr) (Bit False)
+    MG.basicSet (M.slice 1 0 arr) (Bit False)
     (,) <$> M.read arr 0 <*> M.read arr 1
 
-case_write_set_read4 :: IO ()
+case_write_set_read4 :: Property
 case_write_set_read4 =
-  assertEqual "should be equal" (Bit True, Bit True) $ runST $ do
+  (=== (Bit True, Bit True)) $ runST $ do
     arr <- M.new 3
     M.write arr 0 (Bit True)
     M.write arr 2 (Bit True)
-    M.basicSet (M.slice 1 1 arr) (Bit False)
+    MG.basicSet (M.slice 1 1 arr) (Bit False)
     (,) <$> M.read arr 0 <*> M.read arr 2
 
-case_set_read1 :: IO ()
-case_set_read1 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_set_read1 :: Property
+case_set_read1 = (=== Bit True) $ runST $ do
   arr <- M.new 1
-  M.basicSet arr (Bit True)
+  MG.basicSet arr (Bit True)
   M.read arr 0
 
-case_set_read2 :: IO ()
-case_set_read2 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_set_read2 :: Property
+case_set_read2 = (=== Bit True) $ runST $ do
   arr <- M.new 2
-  M.basicSet (M.slice 1 1 arr) (Bit True)
+  MG.basicSet (M.slice 1 1 arr) (Bit True)
   M.read arr 1
 
-case_set_read3 :: IO ()
-case_set_read3 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_set_read3 :: Property
+case_set_read3 = (=== Bit True) $ runST $ do
   arr <- M.new 192
-  M.basicSet (M.slice 71 121 arr) (Bit True)
+  MG.basicSet (M.slice 71 121 arr) (Bit True)
   M.read arr 145
 
-case_set_read4 :: IO ()
-case_set_read4 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_set_read4 :: Property
+case_set_read4 = (=== Bit True) $ runST $ do
   arr <- M.slice 27 38 <$> M.new 65
-  M.basicSet arr (Bit True)
+  MG.basicSet arr (Bit True)
   M.read arr 21
 
-case_write_copy_read1 :: IO ()
-case_write_copy_read1 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_write_copy_read1 :: Property
+case_write_copy_read1 = (=== Bit True) $ runST $ do
   src <- M.slice 37 28 <$> M.new 65
   M.write src 27 (Bit True)
   dst <- M.slice 37 28 <$> M.new 65
   M.copy dst src
   M.read dst 27
 
-case_write_copy_read2 :: IO ()
-case_write_copy_read2 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_write_copy_read2 :: Property
+case_write_copy_read2 = (=== Bit True) $ runST $ do
   src <- M.slice 32 33 <$> M.new 65
   M.write src 0 (Bit True)
   dst <- M.slice 32 33 <$> M.new 65
   M.copy dst src
   M.read dst 0
 
-case_write_copy_read3 :: IO ()
-case_write_copy_read3 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_write_copy_read3 :: Property
+case_write_copy_read3 = (=== Bit True) $ runST $ do
   src <- M.slice 1 1 <$> M.new 2
   M.write src 0 (Bit True)
   dst <- M.slice 1 1 <$> M.new 2
   M.copy dst src
   M.read dst 0
 
-case_write_copy_read4 :: IO ()
-case_write_copy_read4 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_write_copy_read4 :: Property
+case_write_copy_read4 = (=== Bit True) $ runST $ do
   src <- M.slice 12 52 <$> M.new 64
   M.write src 22 (Bit True)
   dst <- M.slice 12 52 <$> M.new 64
   M.copy dst src
   M.read dst 22
 
-case_write_copy_read5 :: IO ()
-case_write_copy_read5 = assertEqual "should be equal" (Bit True) $ runST $ do
+case_write_copy_read5 :: Property
+case_write_copy_read5 = (=== Bit True) $ runST $ do
   src <- M.slice 48 80 <$> M.new 128
   M.write src 46 (Bit True)
   dst <- M.slice 48 80 <$> M.new 128
@@ -198,34 +215,103 @@
 prop_slice_def
   :: NonNegative Int
   -> NonNegative Int
-  -> N.New B.Vector Bit
+  -> N.New U.Vector Bit
   -> Property
 prop_slice_def (NonNegative s) (NonNegative n) xs =
   l > 0 ==> runST $ do
     let xs' = V.new xs
     xs1 <- N.run xs
     xs2 <- V.unsafeFreeze (M.slice s' n' xs1)
-    return (B.toList xs2 === sliceList s' n' (B.toList xs'))
+    return (U.toList xs2 === sliceList s' n' (U.toList xs'))
   where
     l = V.length (V.new xs)
     s' = s `mod` l
     n' = n `mod` (l - s')
 
-prop_grow_def :: B.Vector Bit -> NonNegative Int -> Bool
+prop_grow_def :: U.Vector Bit -> NonNegative Int -> Bool
 prop_grow_def xs (NonNegative m) = runST $ do
-  let n = B.length xs
-  v0  <- B.thaw xs
+  let n = U.length xs
+  v0  <- U.thaw xs
   v1  <- M.grow v0 m
-  fv0 <- B.freeze v0
-  fv1 <- B.freeze v1
-  return (fv0 == B.take n fv1)
+  fv0 <- U.freeze v0
+  fv1 <- U.freeze v1
+  return (fv0 == U.take n fv1)
 
-prop_cloneFromWords_def :: N.New B.Vector Word -> Property
-prop_cloneFromWords_def ws =
-  runST (N.run ws >>= pure . castFromWordsM >>= V.unsafeFreeze)
+prop_castFromWords_def :: N.New U.Vector Word -> Property
+prop_castFromWords_def ws =
+  runST (N.run ws >>= V.unsafeFreeze . castFromWordsM)
     === castFromWords (V.new ws)
 
-prop_cloneToWords_def :: N.New B.Vector Bit -> Property
+prop_cloneToWords_def :: N.New U.Vector Bit -> Property
 prop_cloneToWords_def xs =
   runST (N.run xs >>= cloneToWordsM >>= V.unsafeFreeze)
     === cloneToWords (V.new xs)
+
+prop_castToWords_1 :: N.New U.Vector Word -> Property
+prop_castToWords_1 xs = runST $ do
+  vs <- N.run xs
+  vs' <- cloneToWordsM (castFromWordsM vs)
+  case castToWordsM (castFromWordsM vs) of
+    Nothing -> pure $ property False
+    Just vs'' -> do
+      ws'  <- V.unsafeFreeze vs'
+      ws'' <- V.unsafeFreeze vs''
+      pure $ ws' === ws''
+
+prop_castToWords_2 :: N.New U.Vector Bit -> Property
+prop_castToWords_2 xs = runST $ do
+  vs <- N.run xs
+  case castToWordsM vs of
+    Nothing  -> pure $ property True
+    Just ws -> do
+      ws' <- V.unsafeFreeze (castFromWordsM ws)
+      ws'' <- V.unsafeFreeze vs
+      pure $ ws' === ws''
+
+prop_replicate_neg :: Positive Int -> Bit -> Property
+prop_replicate_neg (Positive n) x = ioProperty $ do
+  ret <- try (evaluate (runST $ MG.basicUnsafeReplicate (-n) x >>= U.unsafeFreeze))
+  pure $ property $ case ret of
+    Left (_ :: SomeException) -> True
+    _ -> False
+
+prop_new_neg :: Positive Int -> Property
+prop_new_neg (Positive n) = ioProperty $ do
+  ret <- try (evaluate (runST $ MG.basicUnsafeNew (-n) >>= U.unsafeFreeze :: U.Vector Bit))
+  pure $ property $ case ret of
+    Left (_ :: SomeException) -> True
+    _ -> False
+
+prop_move_around :: Int -> Int -> Int -> U.Vector Bit -> Property
+prop_move_around srcFrom' sliceLen' dstFrom' xs = ioProperty $ do
+  let l = U.length xs
+  when (l < 1) discard
+  let sliceLen = sliceLen' `mod` l
+      srcFrom = srcFrom' `mod` (l - sliceLen)
+      dstFrom = dstFrom' `mod` (l - sliceLen)
+  ys <- V.thaw xs
+  let src = M.slice srcFrom sliceLen ys
+      dst = M.slice dstFrom sliceLen ys
+  M.move dst src
+  xs' <- V.unsafeFreeze ys
+  let slicePrefix = V.slice 0 dstFrom
+      sliceSuffix = V.slice (dstFrom + sliceLen) (l - dstFrom - sliceLen)
+  pure $ slicePrefix xs === slicePrefix xs' .&&. sliceSuffix xs === sliceSuffix xs'
+
+prop_copy_around :: Bool -> Int -> Int -> Int -> U.Vector Bit -> Property
+prop_copy_around b srcFrom' sliceLen' dstFrom' xs = ioProperty $ do
+  let l = U.length xs
+  when (l < 4) discard
+  let sliceLen = 1 + sliceLen' `mod` ((l - 2) `quot` 2)
+      firstFrom = srcFrom' `mod` (l - 2 * sliceLen)
+      secondFrom = firstFrom + sliceLen + dstFrom' `mod` (l - firstFrom - 2 * sliceLen)
+      (srcFrom, dstFrom) = if b then (firstFrom, secondFrom) else (secondFrom, firstFrom)
+  ys <- V.thaw xs
+  let src = M.slice srcFrom sliceLen ys
+      dst = M.slice dstFrom sliceLen ys
+  when (M.overlaps src dst) discard
+  M.copy dst src
+  xs' <- V.unsafeFreeze ys
+  let slicePrefix = V.slice 0 dstFrom
+      sliceSuffix = V.slice (dstFrom + sliceLen) (l - dstFrom - sliceLen)
+  pure $ slicePrefix xs === slicePrefix xs' .&&. sliceSuffix xs === sliceSuffix xs'
diff --git a/test/Tests/SetOps.hs b/test/Tests/SetOps.hs
--- a/test/Tests/SetOps.hs
+++ b/test/Tests/SetOps.hs
@@ -2,12 +2,12 @@
 {-# LANGUAGE RankNTypes #-}
 
 #ifndef BITVEC_THREADSAFE
-module Tests.SetOps where
+module Tests.SetOps (setOpTests) where
 #else
-module Tests.SetOpsTS where
+module Tests.SetOpsTS (setOpTests) where
 #endif
 
-import Support ()
+import Support (twoTimesMore)
 
 import Control.Monad
 import Control.Monad.ST
@@ -19,34 +19,88 @@
 import Test.Tasty.QuickCheck hiding ((.&.))
 
 setOpTests :: TestTree
-setOpTests = testGroup
-  "Set operations"
-  [ testProperty "generalize"       prop_generalize
-  , testProperty "zipBits"          prop_zipBits
-  , testProperty "zipInPlace"       prop_zipInPlace
-  , testProperty "invertBits"       prop_invertBits
-  , testProperty "invertBitsWords"  prop_invertBitsWords
-  , testProperty "invertInPlace"    prop_invertInPlace
-  , testProperty "invertInPlace middle" prop_invertInPlace_middle
-  , testProperty "reverseBits"      prop_reverseBits
-  , testProperty "reverseBitsWords" prop_reverseBitsWords
-  , testProperty "reverseInPlace"   prop_reverseInPlace
-  , testProperty "reverseInPlace middle" prop_reverseInPlace_middle
-  , testProperty "selectBits"       prop_selectBits_def
-  , testProperty "excludeBits"      prop_excludeBits_def
-  , testProperty "countBits"        prop_countBits_def
+setOpTests = testGroup "Set operations"
+  [ testProperty "generalize1"              prop_generalize1
+  , testProperty "generalize2"              prop_generalize2
+  , twoTimesMore
+  $ testProperty "zipBits"                  prop_zipBits
+  , testProperty "zipInPlace"               prop_zipInPlace
+  , testProperty "zipInPlace around"        prop_zipInPlace_around
+
+  , testProperty "mapBits"                  prop_mapBits
+  , testProperty "mapInPlace"               prop_mapInPlace
+
+  , testProperty "union"                    prop_union_def
+  , testProperty "intersection"             prop_intersection_def
+  , testProperty "difference"               prop_difference_def
+  , testProperty "symDiff"                  prop_symDiff_def
+
+  , mkGroup "invertBits" prop_invertBits
+
+  , testProperty "invertInPlace"            prop_invertInPlace
+  , testProperty "invertInPlace around"     prop_invertInPlace_around
+  , testProperty "invertInPlaceWords"       prop_invertInPlaceWords
+  , testProperty "invertInPlace middle"     prop_invertInPlace_middle
+  , testProperty "invertInPlaceLong middle" prop_invertInPlaceLong_middle
+
+  , mkGroup "reverseBits" prop_reverseBits
+
+  , testProperty "reverseInPlace"            prop_reverseInPlace
+  , testProperty "reverseInPlaceWords"       prop_reverseInPlaceWords
+  , testProperty "reverseInPlace middle"     prop_reverseInPlace_middle
+  , testProperty "reverseInPlaceLong middle" prop_reverseInPlaceLong_middle
+
+  , mkGroup2 "selectBits"  prop_selectBits_def
+  , mkGroup2 "excludeBits" prop_excludeBits_def
+
+  , mkGroup "countBits" prop_countBits_def
   ]
 
-prop_generalize :: Fun (Bit, Bit) Bit -> Bit -> Bit -> Property
-prop_generalize fun x y = curry (applyFun fun) x y === generalize (curry (applyFun fun)) x y
+mkGroup :: String -> (U.Vector Bit -> Property) -> TestTree
+mkGroup name prop = testGroup name
+  [ testProperty "simple" prop
+  , testProperty "simple_long" (prop . getLarge)
+  , 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 from) (NonNegative len) (NonNegative excess) =
+      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 =
+  applyFun fun x === generalize1 (applyFun fun) x
+
+prop_generalize2 :: Fun (Bit, Bit) Bit -> Bit -> Bit -> Property
+prop_generalize2 fun x y =
+  curry (applyFun fun) x y === generalize2 (curry (applyFun fun)) x y
+
 prop_union_def :: U.Vector Bit -> U.Vector Bit -> Property
 prop_union_def xs ys =
-  zipBits (.|.) xs ys === U.zipWith (.|.) xs ys
+  xs .|. ys === U.zipWith (.|.) xs ys
 
 prop_intersection_def :: U.Vector Bit -> U.Vector Bit -> Property
 prop_intersection_def xs ys =
-  zipBits (.&.) xs ys === U.zipWith (.&.) xs ys
+  xs .&. ys === U.zipWith (.&.) xs ys
 
 prop_difference_def :: U.Vector Bit -> U.Vector Bit -> Property
 prop_difference_def xs ys =
@@ -56,34 +110,62 @@
 
 prop_symDiff_def :: U.Vector Bit -> U.Vector Bit -> Property
 prop_symDiff_def xs ys =
-  zipBits xor xs ys === U.zipWith xor xs ys
+  xs `xor` ys === U.zipWith xor xs ys
 
 prop_zipBits :: Fun (Bit, Bit) Bit -> U.Vector Bit -> U.Vector Bit -> Property
 prop_zipBits fun xs ys =
-  U.zipWith f xs ys === zipBits (generalize f) xs ys
+  U.zipWith f xs ys === zipBits (generalize2 f) xs ys
   where
     f = curry $ applyFun fun
 
 prop_zipInPlace :: Fun (Bit, Bit) Bit -> U.Vector Bit -> U.Vector Bit -> Property
 prop_zipInPlace fun xs ys =
-  U.zipWith f xs ys === U.take (min (U.length xs) (U.length ys)) (U.modify (zipInPlace (generalize f) xs) ys)
+  U.zipWith f xs ys === U.take (min (U.length xs) (U.length ys)) (U.modify (zipInPlace (generalize2 f) xs) ys)
   where
     f = curry $ applyFun fun
 
+prop_zipInPlace_around :: Fun (Bit, Bit) Bit -> U.Vector Bit -> U.Vector Bit -> Property
+prop_zipInPlace_around fun ts xs = ioProperty $ do
+  let l = U.length xs
+      f = curry $ applyFun fun
+  when (l < 2) discard
+  ys <- U.thaw xs
+  let zs = MU.slice 1 (l - 2) ys
+  zipInPlace (generalize2 f) ts zs
+  hd <- MU.read ys 0
+  lst <- MU.read ys (MU.length ys - 1)
+  pure $ hd === U.head xs .&&. lst === U.last xs
+
+prop_mapBits :: Fun Bit Bit -> U.Vector Bit -> Property
+prop_mapBits fun xs =
+  U.map (applyFun fun) xs === mapBits (generalize1 (applyFun fun)) xs
+
+prop_mapInPlace :: Fun Bit Bit -> U.Vector Bit -> Property
+prop_mapInPlace fun xs =
+  U.map (applyFun fun) xs === U.modify (mapInPlace (generalize1 (applyFun fun))) xs
+
 prop_invertBits :: U.Vector Bit -> Property
 prop_invertBits xs =
-  U.map complement xs === invertBits xs
-
-prop_invertBitsWords :: U.Vector Word -> Property
-prop_invertBitsWords ws =
-  U.map complement xs === invertBits xs
-  where
-    xs = castFromWords ws
+  U.map complement xs === complement xs
 
 prop_invertInPlace :: U.Vector Bit -> Property
 prop_invertInPlace xs =
   U.map complement xs === U.modify invertInPlace xs
 
+prop_invertInPlace_around :: U.Vector Bit -> Property
+prop_invertInPlace_around xs = ioProperty $ do
+  let l = U.length xs
+  when (l < 2) discard
+  ys <- U.thaw xs
+  let zs = MU.slice 1 (l - 2) ys
+  invertInPlace zs
+  hd <- MU.read ys 0
+  lst <- MU.read ys (MU.length ys - 1)
+  pure $ hd === U.head xs .&&. lst === U.last xs
+
+prop_invertInPlaceWords :: Large (U.Vector Bit) -> Property
+prop_invertInPlaceWords = prop_invertInPlace . getLarge
+
 prop_invertInPlace_middle :: NonNegative Int -> NonNegative Int -> NonNegative Int -> Property
 prop_invertInPlace_middle (NonNegative from) (NonNegative len) (NonNegative excess) = runST $ do
   let totalLen = from + len + excess
@@ -96,26 +178,29 @@
   invertInPlace middle
   wec <- U.unsafeFreeze vec
 
-  let refLeft  = U.take from ref
-      wecLeft  = U.take from wec
-      refRight = U.drop (from + len) ref
-      wecRight = U.drop (from + len) wec
-  pure $ refLeft === wecLeft .&&. refRight === wecRight
+  let refLeft   = U.take from ref
+      wecLeft   = U.take from wec
+      refRight  = U.drop (from + len) ref
+      wecRight  = U.drop (from + len) wec
+      refMiddle = U.map complement (U.take len (U.drop from ref))
+      wecMiddle = U.take len (U.drop from wec)
+  pure $ refLeft === wecLeft .&&. refRight === wecRight .&&. refMiddle === wecMiddle
 
+prop_invertInPlaceLong_middle :: NonNegative Int -> NonNegative Int -> NonNegative Int -> Property
+prop_invertInPlaceLong_middle (NonNegative x) (NonNegative y) (NonNegative z) =
+  prop_invertInPlace_middle (NonNegative $ x * 31) (NonNegative $ y * 37) (NonNegative $ z * 29)
+
 prop_reverseBits :: U.Vector Bit -> Property
 prop_reverseBits xs =
   U.reverse xs === reverseBits xs
 
-prop_reverseBitsWords :: U.Vector Word -> Property
-prop_reverseBitsWords ws =
-  U.reverse xs === reverseBits xs
-  where
-    xs = castFromWords ws
-
 prop_reverseInPlace :: U.Vector Bit -> Property
 prop_reverseInPlace xs =
   U.reverse xs === U.modify reverseInPlace xs
 
+prop_reverseInPlaceWords :: Large (U.Vector Bit) -> Property
+prop_reverseInPlaceWords = prop_reverseInPlace . getLarge
+
 prop_reverseInPlace_middle :: NonNegative Int -> NonNegative Int -> NonNegative Int -> Property
 prop_reverseInPlace_middle (NonNegative from) (NonNegative len) (NonNegative excess) = runST $ do
   let totalLen = from + len + excess
@@ -128,12 +213,18 @@
   reverseInPlace middle
   wec <- U.unsafeFreeze vec
 
-  let refLeft  = U.take from ref
-      wecLeft  = U.take from wec
-      refRight = U.drop (from + len) ref
-      wecRight = U.drop (from + len) wec
-  pure $ refLeft === wecLeft .&&. refRight === wecRight
+  let refLeft   = U.take from ref
+      wecLeft   = U.take from wec
+      refRight  = U.drop (from + len) ref
+      wecRight  = U.drop (from + len) wec
+      refMiddle = U.reverse (U.take len (U.drop from ref))
+      wecMiddle = U.take len (U.drop from wec)
+  pure $ refLeft === wecLeft .&&. refRight === wecRight .&&. refMiddle === wecMiddle
 
+prop_reverseInPlaceLong_middle :: NonNegative Int -> NonNegative Int -> NonNegative Int -> Property
+prop_reverseInPlaceLong_middle (NonNegative x) (NonNegative y) (NonNegative z) =
+  prop_reverseInPlace_middle (NonNegative $ x * 31) (NonNegative $ y * 37) (NonNegative $ z * 29)
+
 select :: U.Unbox a => U.Vector Bit -> U.Vector a -> U.Vector a
 select mask ws = U.map snd (U.filter (unBit . fst) (U.zip mask ws))
 
@@ -151,8 +242,15 @@
 
 -------------------------------------------------------------------------------
 
-generalize :: (Bit -> Bit -> Bit) -> (forall a. Bits a => a -> a -> a)
-generalize f = case (f (Bit False) (Bit False), f (Bit False) (Bit True), f (Bit True) (Bit False), f (Bit True) (Bit True)) of
+generalize1 :: (Bit -> Bit) -> (forall a. Bits a => a -> a)
+generalize1 f = case (f (Bit False), f (Bit True)) of
+  (Bit False, Bit False) -> const zeroBits
+  (Bit False, Bit True)  -> id
+  (Bit True,  Bit False) -> complement
+  (Bit True,  Bit True)  -> const $ complement zeroBits
+
+generalize2 :: (Bit -> Bit -> Bit) -> (forall a. Bits a => a -> a -> a)
+generalize2 f = case (f (Bit False) (Bit False), f (Bit False) (Bit True), f (Bit True) (Bit False), f (Bit True) (Bit True)) of
   (Bit False, Bit False, Bit False, Bit False) -> \_ _ -> zeroBits
   (Bit False, Bit False, Bit False, Bit True)  -> \x y -> x .&. y
   (Bit False, Bit False, Bit True,  Bit False) -> \x y -> x .&. complement y
diff --git a/test/Tests/Vector.hs b/test/Tests/Vector.hs
--- a/test/Tests/Vector.hs
+++ b/test/Tests/Vector.hs
@@ -1,100 +1,201 @@
-module Tests.Vector where
+{-# LANGUAGE CPP #-}
 
+module Tests.Vector
+  ( vectorTests
+  ) where
+
 import Support
 
 import Prelude hiding (and, or)
+import Control.Exception
 import Data.Bit
 import Data.Bits
-import Data.List hiding (and, or)
-import qualified Data.Vector.Unboxed as U hiding (reverse, and, or, any, all, findIndex)
+import Data.List (findIndex)
+import qualified Data.Vector.Primitive as P
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Base as UB
 import Data.Word
 import Test.Tasty
-import Test.Tasty.HUnit
-import Test.Tasty.QuickCheck
+import Test.Tasty.QuickCheck (Property, NonNegative(..), Positive(..), testProperty, Large(..), (===), property, once, (==>), ioProperty, (.&&.), counterexample)
 
+#include "MachDeps.h"
+
 vectorTests :: TestTree
 vectorTests = testGroup "Data.Vector.Unboxed.Bit"
   [ testGroup "Data.Vector.Unboxed functions"
     [ testProperty "toList . fromList == id" prop_toList_fromList
-    , testProperty "fromList . toList == id" prop_fromList_toList
+    , mkGroup      "fromList . toList == id" prop_fromList_toList
     , testProperty "slice"                   prop_slice_def
     ]
   , tenTimesLess $
     testProperty "cloneFromWords" prop_cloneFromWords_def
-  , testProperty "cloneToWords"   prop_cloneToWords_def
+  , mkGroup      "cloneToWords"   prop_cloneToWords_def
   , tenTimesLess $
+    testProperty "castToWords_1"   prop_castToWords_1
+  , tenTimesLess $
+    testProperty "castToWords_2"   prop_castToWords_2
+  , tenTimesLess $
     testProperty "cloneFromWords8" prop_cloneFromWords8_def
-  , testProperty "cloneToWords8"   prop_cloneToWords8_def
-  , testProperty "reverse"        prop_reverse_def
-  , testProperty "countBits"      prop_countBits_def
-  , testProperty "listBits"       prop_listBits_def
-  , testGroup "Boolean operations"
-    [ testProperty "and" prop_and_def
-    , testProperty "or" prop_or_def
+  , mkGroup      "cloneToWords8"   prop_cloneToWords8_def
+  , tenTimesLess $
+    testProperty "castToWords8_1"  prop_castToWords8_1
+  , tenTimesLess $
+    testProperty "castToWords8_2"  prop_castToWords8_2
+  , testProperty "cloneToByteString" prop_cloneToByteString
+  , mkGroup "reverse"        prop_reverse_def
+  , testGroup "countBits"
+    [ testProperty "special case 1" case_countBits_1
+    , mkGroup "matches definition"  prop_countBits_def
     ]
-  , testGroup "Search operations"
-    [ testProperty "first" prop_first_def
+  , testGroup "listBits"
+    [ testProperty "special case 1" case_listBits_1
+    , testProperty "special case 2" case_listBits_2
+    , mkGroup "matches definition"  prop_listBits_def
     ]
+  , mkGroup "and"            prop_and_def
+  , mkGroup "or"             prop_or_def
+  , testGroup "bitIndex"
+    [ testProperty "special case 1" case_bitIndex_1
+    , testProperty "special case 2" case_bitIndex_2
+    , testProperty "special case 3" case_bitIndex_3
+    , testProperty "special case 4" case_bitIndex_4
+    , testProperty "special case 5" case_bitIndex_5
+    , testProperty "special case 6" case_bitIndex_6
+    , testProperty "special case 7" case_bitIndex_7
+    , mkGroup "True"               (prop_bitIndex_1 (Bit True))
+    , mkGroup "False"              (prop_bitIndex_1 (Bit False))
+    ]
   , testGroup "nthBitIndex"
-    [ testCase "special case 1" case_nthBit_1
-    , testProperty "matches bitIndex True"              prop_nthBit_1
-    , testProperty "matches bitIndex False"             prop_nthBit_2
+    [ testProperty "special case 1"                     case_nthBit_1
+    , testProperty "special case 2"                     case_nthBit_2
+    , testProperty "special case 3"                     case_nthBit_3
+    , testProperty "special case 4"                     case_nthBit_4
+    , testProperty "special case 5"                     case_nthBit_5
+    , testProperty "special case 6"                     case_nthBit_6
+    , testProperty "special case 7"                     case_nthBit_7
+    , mkGroup      "matches bitIndex True"              prop_nthBit_1
+    , mkGroup      "matches bitIndex False"             prop_nthBit_2
     , testProperty "matches sequence of bitIndex True"  prop_nthBit_3
     , testProperty "matches sequence of bitIndex False" prop_nthBit_4
     , testProperty "matches countBits"                  prop_nthBit_5
+    , testProperty "negative argument"                  prop_nthBit_6
     ]
   , testGroup "Bits instance"
     [ testProperty "rotate is reversible" prop_rotate
     , testProperty "bit"                  prop_bit
     , testProperty "shiftL"               prop_shiftL
     , testProperty "shiftR"               prop_shiftR
+    , testProperty "zeroBits"             prop_zeroBits
+    , testProperty "bitSize"              prop_bitSize
+    , testProperty "isSigned"             prop_isSigned
+    , testProperty "setBit"               prop_setBit
+    , testProperty "clearBit"             prop_clearBit
+    , testProperty "complementBit"        prop_complementBit
     ]
   ]
 
-prop_toList_fromList :: [Bit] -> Bool
-prop_toList_fromList xs = U.toList (U.fromList xs) == xs
+mkGroup :: String -> (U.Vector Bit -> Property) -> TestTree
+mkGroup name prop = testGroup name
+  [ testProperty "simple" prop
+  , testProperty "simple_long" (prop . getLarge)
+  , 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 from) (NonNegative len) (NonNegative excess) =
+      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)
 
-prop_fromList_toList :: U.Vector Bit -> Bool
-prop_fromList_toList xs = U.fromList (U.toList xs) == xs
+prop_toList_fromList :: [Bit] -> Property
+prop_toList_fromList xs = U.toList (U.fromList xs) === xs
 
-prop_slice_def :: Int -> Int -> U.Vector Bit -> Bool
-prop_slice_def s n xs = sliceList s' n' (U.toList xs)
-  == U.toList (U.slice s' n' xs)
-  where (s', n') = trimSlice s n (U.length xs)
+prop_fromList_toList :: U.Vector Bit -> Property
+prop_fromList_toList xs = U.fromList (U.toList xs) === xs
 
+prop_slice_def :: Int -> Int -> U.Vector Bit -> Property
+prop_slice_def s n xs =
+  sliceList s' n' (U.toList xs) === U.toList (U.slice s' n' xs)
+  where
+    (s', n') = trimSlice s n (U.length xs)
+
 prop_cloneFromWords_def :: U.Vector Word -> Property
 prop_cloneFromWords_def ws =
   U.toList (castFromWords ws) === concatMap wordToBitList (U.toList ws)
 
-prop_cloneToWords_def :: U.Vector Bit -> Bool
-prop_cloneToWords_def xs = U.toList (cloneToWords xs) == loop (U.toList xs)
+prop_cloneToWords_def :: U.Vector Bit -> Property
+prop_cloneToWords_def xs = U.toList (cloneToWords xs) === loop (U.toList xs)
  where
   loop [] = []
   loop bs = case packBitsToWord bs of
     (w, bs') -> w : loop bs'
 
+prop_castToWords_1 :: U.Vector Word -> Property
+prop_castToWords_1 ws =
+  Just ws === castToWords (castFromWords ws)
+
+prop_castToWords_2 :: U.Vector Bit -> Property
+prop_castToWords_2 xs = case castToWords xs of
+  Nothing -> property True
+  Just ws -> castFromWords ws === xs
+
 prop_cloneFromWords8_def :: U.Vector Word8 -> Property
-prop_cloneFromWords8_def ws =
-  U.toList (castFromWords8 ws) === concatMap wordToBitList (U.toList ws)
+prop_cloneFromWords8_def ws
+  = counterexample ("offset = " ++ show off ++ " len = " ++ show len)
+  $ U.toList (castFromWords8 ws) === concatMap wordToBitList (U.toList ws)
+  where
+    UB.V_Word8 (P.Vector off len _) = ws
 
-prop_cloneToWords8_def :: U.Vector Bit -> Bool
-prop_cloneToWords8_def xs = U.toList (cloneToWords8 xs) == loop (U.toList xs)
- where
-  loop [] = []
-  loop bs = case packBitsToWord bs of
-    (w, bs') -> w : loop bs'
+prop_cloneToWords8_def :: U.Vector Bit -> Property
+prop_cloneToWords8_def xs@(BitVec off len _)
+  = counterexample ("offset = " ++ show off ++ " len = " ++ show len)
+  $ U.toList (cloneToWords8 xs) === loop (U.toList xs)
+  where
+    loop [] = []
+    loop bs = case packBitsToWord bs of
+      (w, bs') -> w : loop bs'
 
-prop_reverse_def :: U.Vector Bit -> Bool
+prop_castToWords8_1 :: U.Vector Word8 -> Property
+#ifdef WORDS_BIGENDIAN
+prop_castToWords8_1 ws = Nothing === castToWords8 (castFromWords8 ws)
+#else
+prop_castToWords8_1 ws
+  = counterexample ("offset = " ++ show off ++ " len = " ++ show len)
+  $ Just ws === castToWords8 (castFromWords8 ws)
+  where
+    UB.V_Word8 (P.Vector off len _) = ws
+#endif
+
+prop_castToWords8_2 :: U.Vector Bit -> Property
+prop_castToWords8_2 xs = case castToWords8 xs of
+  Nothing -> property True
+  Just ws -> castFromWords8 ws === xs
+
+prop_reverse_def :: U.Vector Bit -> Property
 prop_reverse_def xs =
-  reverse (U.toList xs) == U.toList (U.modify reverseInPlace xs)
+  reverse (U.toList xs) === U.toList (U.modify reverseInPlace xs)
 
-prop_countBits_def :: U.Vector Bit -> Bool
-prop_countBits_def xs = countBits xs == length (filter unBit (U.toList xs))
+prop_countBits_def :: U.Vector Bit -> Property
+prop_countBits_def xs = countBits xs === length (filter unBit (U.toList xs))
 
+case_countBits_1 :: Property
+case_countBits_1 = once $
+  countBits (U.drop 64 (U.replicate 128 (Bit False))) === 0
+
 prop_listBits_def :: U.Vector Bit -> Property
 prop_listBits_def xs =
   listBits xs === [ i | (i, x) <- zip [0 ..] (U.toList xs), unBit x ]
 
+case_listBits_1 :: Property
+case_listBits_1 = once $
+  listBits (U.drop 24 (U.replicate 64 (Bit False))) === []
+
+case_listBits_2 :: Property
+case_listBits_2 = once $
+  listBits (U.drop 24 (U.replicate 128 (Bit True))) === [0..103]
+
 and :: U.Vector Bit -> Bool
 and xs = case bitIndex (Bit False) xs of
   Nothing -> True
@@ -111,9 +212,37 @@
 prop_or_def :: U.Vector Bit -> Property
 prop_or_def xs = or xs === any unBit (U.toList xs)
 
-prop_first_def :: Bit -> U.Vector Bit -> Bool
-prop_first_def b xs = bitIndex b xs == findIndex (b ==) (U.toList xs)
+case_bitIndex_1 :: Property
+case_bitIndex_1 = once $
+  bitIndex (Bit True) (U.generate 128 (Bit . (== 64))) === Just 64
 
+case_bitIndex_2 :: Property
+case_bitIndex_2 = once $
+  bitIndex (Bit False) (U.generate 128 (Bit . (/= 64))) === Just 64
+
+case_bitIndex_3 :: Property
+case_bitIndex_3 = once $
+  bitIndex (Bit True) (U.drop 63 (U.generate 128 (Bit . (== 64)))) === Just 1
+
+case_bitIndex_4 :: Property
+case_bitIndex_4 = once $
+  bitIndex (Bit False) (U.drop 63 (U.generate 128 (Bit . (/= 64)))) === Just 1
+
+case_bitIndex_5 :: Property
+case_bitIndex_5 = once $
+  bitIndex (Bit False) (U.drop 63 (U.replicate 65 (Bit True))) === Nothing
+
+case_bitIndex_6 :: Property
+case_bitIndex_6 = once $
+  bitIndex (Bit False) (U.drop 63 (U.generate 66 (Bit . (== 63)))) === Just 1
+
+case_bitIndex_7 :: Property
+case_bitIndex_7 = once $
+  bitIndex (Bit False) (U.drop 1023 (U.generate 1097 (Bit . (/= 1086)))) === Just 63
+
+prop_bitIndex_1 :: Bit -> U.Vector Bit -> Property
+prop_bitIndex_1 b xs = bitIndex b xs === findIndex (b ==) (U.toList xs)
+
 prop_nthBit_1 :: U.Vector Bit -> Property
 prop_nthBit_1 xs = bitIndex (Bit True) xs === nthBitIndex (Bit True) 1 xs
 
@@ -147,18 +276,50 @@
     count = countBits xs
     n' = n `mod` count + 1
 
-case_nthBit_1 :: IO ()
-case_nthBit_1 =
-  assertEqual "should be equal" Nothing
-    $ nthBitIndex (Bit True) 1
-    $ U.slice 61 4
-    $ U.replicate 100 (Bit False)
+prop_nthBit_6 :: NonNegative Int -> U.Vector Bit -> Property
+prop_nthBit_6 (NonNegative n) xs = ioProperty $ do
+  ret <- try (evaluate (nthBitIndex (Bit True) (-n) xs))
+  pure $ property $ case ret of
+    Left ErrorCall{} -> True
+    _ -> False
 
+case_nthBit_1 :: Property
+case_nthBit_1 = once $
+  nthBitIndex (Bit True) 1 (U.slice 61 4 (U.replicate 100 (Bit False))) === Nothing
+
+case_nthBit_2 :: Property
+case_nthBit_2 = once $
+  nthBitIndex (Bit False) 1 (U.slice 61 4 (U.replicate 100 (Bit True))) === Nothing
+
+case_nthBit_3 :: Property
+case_nthBit_3 = once $
+  nthBitIndex (Bit True) 1 (U.drop 63 (U.generate 128 (Bit . (== 64)))) === Just 1
+
+case_nthBit_4 :: Property
+case_nthBit_4 = once $
+  nthBitIndex (Bit False) 1 (U.drop 63 (U.generate 128 (Bit . (/= 64)))) === Just 1
+
+case_nthBit_5 :: Property
+case_nthBit_5 = once $
+  nthBitIndex (Bit False) 1 (U.drop 63 (U.replicate 65 (Bit True))) === Nothing
+
+case_nthBit_6 :: Property
+case_nthBit_6 = once $
+  nthBitIndex (Bit False) 1 (U.drop 63 (U.generate 66 (Bit . (== 63)))) === Just 1
+
+case_nthBit_7 :: Property
+case_nthBit_7 = once $
+  nthBitIndex (Bit False) 1 (U.drop 1023 (U.generate 1097 (Bit . (/= 1086)))) === Just 63
+
 prop_rotate :: Int -> U.Vector Bit -> Property
 prop_rotate n v = v === (v `rotate` n) `rotate` (-n)
 
-prop_bit :: NonNegative Int -> Property
-prop_bit (NonNegative n) = testBit v n .&&. popCount v === 1 .&&. U.length v == n + 1
+prop_bit :: Int -> Property
+prop_bit n
+  | n >= 0
+  = testBit v n .&&. popCount v === 1 .&&. U.length v === n + 1
+  | otherwise
+  = not (testBit v n) .&&. popCount v === 0 .&&. U.length v === 0
   where
     v :: U.Vector Bit
     v = bit n
@@ -172,3 +333,29 @@
 prop_shiftR (NonNegative n) v = U.drop n v === U.drop n u .&&. popCount (U.take n u) === 0
   where
     u = (v `shiftR` n) `shiftL` n
+
+prop_zeroBits :: Property
+prop_zeroBits = once $
+  U.length (zeroBits :: U.Vector Bit) === 0
+
+prop_bitSize :: U.Vector Bit -> Property
+prop_bitSize v = bitSizeMaybe v === Nothing
+
+prop_isSigned :: U.Vector Bit -> Property
+prop_isSigned v = isSigned v === False
+
+prop_setBit :: Int -> U.Vector Bit -> Property
+prop_setBit n v = v `setBit` n === U.imap ((.|.) . Bit . (== n)) v
+
+prop_clearBit :: Int -> U.Vector Bit -> Property
+prop_clearBit n v = v `clearBit` n === U.imap ((.&.) . Bit . (/= n)) v
+
+prop_complementBit :: Int -> U.Vector Bit -> Property
+prop_complementBit n v = v `complementBit` n === U.imap (xor . Bit . (== n)) v
+
+prop_cloneToByteString :: U.Vector Bit -> Property
+prop_cloneToByteString v@(BitVec off len _)
+  = counterexample ("offset = " ++ show off ++ " len = " ++ show len)
+  $ cloneToByteString (cloneFromByteString bs) === bs
+  where
+    bs = cloneToByteString v