diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 # 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)
 
-A newtype over `Bool` with a better `Vector` instance.
+A newtype over `Bool` with a better `Vector` instance: 8x less memory, up to 1000x faster.
 
 The [`vector`](https://hackage.haskell.org/package/vector)
 package represents unboxed arrays of `Bool`
@@ -12,7 +12,7 @@
 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 64x faster than `Vector Bool`.
+`Vector Bit` is up to 1000x faster than `Vector Bool`.
 
 ## Thread safety
 
@@ -23,16 +23,6 @@
 * `Data.Bit.ThreadSafe` is slower (up to 20%),
   but writes and flips are thread-safe.
 
-## Similar packages
-
-* [`bv`](https://hackage.haskell.org/package/bv) and
-  [`bv-little`](https://hackage.haskell.org/package/bv-little)
-  do not offer mutable vectors.
-
-* [`array`](https://hackage.haskell.org/package/array)
-  is memory-efficient for `Bool`, but lacks
-  a handy `Vector` interface and is not thread-safe.
-
 ## Quick start
 
 Consider the following (very naive) implementation of
@@ -103,8 +93,9 @@
 25
 ```
 
-And vice-versa, query an address of the _n_-th set bit
+And vice versa, query an address of the _n_-th set bit
 (which corresponds to the _n_-th prime number here):
+
 ```haskell
 > nthBitIndex (Bit True) 10 eratosthenes
 Just 29
@@ -116,3 +107,71 @@
 because the former is thread-unsafe with regards to writes.
 There is a moderate performance penalty (up to 20%)
 for using the thread-safe interface.
+
+## 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.
+
+For example, consider three possible representations of a set of `Word16`:
+
+* As an `IntSet` with a readily available `union` function.
+* 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.
+
+## 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,
+operating on `UPoly Bit`,
+this package provides even faster arithmetic routines
+exposed via `F2Poly` data type and its instances.
+
+```haskell
+> :set -XBinaryLiterals
+> -- (1 + x) (1 + x + x^2) = 1 + x^3 (mod 2)
+> 0b11 * 0b111 :: F2Poly
+F2Poly {unF2Poly = [1,0,0,1]}
+```
+
+Use `fromInteger` / `toInteger` to convert binary polynomials
+from `Integer` to `F2Poly` and back.
+
+## 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.
+
+  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.
+
+* Flag `bmi2`, disabled by default, experimental.
+
+  Depend on `bits-extra` package and use it for `nthBitIndex`.
+  This is supposed to be faster, but have not been properly polished yet.
+
+## Similar packages
+
+* [`bv`](https://hackage.haskell.org/package/bv) and
+  [`bv-little`](https://hackage.haskell.org/package/bv-little)
+  do not offer mutable vectors.
+
+* [`array`](https://hackage.haskell.org/package/array)
+  is memory-efficient for `Bool`, but lacks
+  a handy `Vector` interface and is not thread-safe.
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -5,20 +5,24 @@
 import Bench.BitIndex
 import Bench.Intersection
 import Bench.Invert
+import Bench.Product
 import Bench.RandomFlip
 import Bench.RandomRead
 import Bench.RandomWrite
 import Bench.Reverse
+import Bench.Sum
 import Bench.Union
 
 main :: IO ()
 main = defaultMain
-  [ bgroup "bitIndex"     $ map benchBitIndex      [5..10]
-  , bgroup "invert"       $ map benchInvert        [5..10]
-  , bgroup "intersection" $ map benchIntersection  [5..10]
-  , bgroup "randomWrite"  $ map benchRandomWrite   [5..10]
-  , bgroup "randomFlip"   $ map benchRandomFlip    [5..10]
-  , bgroup "randomRead"   $ map benchRandomRead    [5..10]
-  , bgroup "reverse"      $ map benchReverse       [5..10]
-  , bgroup "union"        $ map benchUnion         [5..10]
+  [ bgroup "bitIndex"     $ map benchBitIndex     [5..14]
+  , bgroup "invert"       $ map benchInvert       [5..14]
+  , bgroup "intersection" $ map benchIntersection [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 "reverse"      $ map benchReverse      [5..14]
+  , bgroup "sum"          $ map benchSum          [5..14]
+  , bgroup "union"        $ map benchUnion        [5..14]
   ]
diff --git a/bench/Bench/Intersection.hs b/bench/Bench/Intersection.hs
--- a/bench/Bench/Intersection.hs
+++ b/bench/Bench/Intersection.hs
@@ -40,12 +40,12 @@
 
 benchIntersection :: Int -> Benchmark
 benchIntersection k = bgroup (show (1 `shiftL` k :: Int))
-  [ bench "Bit/zipBits"    $ nf (intersectionBit    (randomVec Bit k))    (randomVec2 Bit k)
-  , bench "Bit/zipWith"    $ nf (intersectionBit'   (randomVec Bit k))    (randomVec2 Bit k)
-  , bench "Bit.TS/zipBits" $ nf (intersectionBitTS  (randomVec TS.Bit k)) (randomVec2 TS.Bit k)
-  , bench "Bit.TS/zipWith" $ nf (intersectionBitTS' (randomVec TS.Bit k)) (randomVec2 TS.Bit k)
-  , bench "Vector"         $ nf (intersectionVector (randomVec id k))     (randomVec2 id k)
-  , bench "IntSet"         $ nf (intersectionIntSet (randomSet k))        (randomSet2 k)
+  [ 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)
   ]
 
 intersectionBit :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
diff --git a/bench/Bench/Invert.hs b/bench/Bench/Invert.hs
--- a/bench/Bench/Invert.hs
+++ b/bench/Bench/Invert.hs
@@ -24,21 +24,21 @@
 
 benchInvert :: Int -> Benchmark
 benchInvert k = bgroup (show (1 `shiftL` k :: Int))
-  [ bench "Bit/invertInPlace"     $ nf invertBit    (randomVec Bit k)
+  [ bench "Bit/invertBits"        $ nf invertBit    (randomVec Bit k)
   , bench "Bit/map-complement"    $ nf invertBit'   (randomVec Bit k)
-  , bench "Bit.TS/invertInPlace"  $ nf invertBitTS  (randomVec TS.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)
   ]
 
 invertBit :: U.Vector Bit -> U.Vector Bit
-invertBit = U.modify invertInPlace
+invertBit = invertBits
 
 invertBit' :: U.Vector Bit -> U.Vector Bit
 invertBit' = U.map complement
 
 invertBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit
-invertBitTS = U.modify TS.invertInPlace
+invertBitTS = TS.invertBits
 
 invertBitTS' :: U.Vector TS.Bit -> U.Vector TS.Bit
 invertBitTS' = U.map complement
diff --git a/bench/Bench/Product.hs b/bench/Bench/Product.hs
new file mode 100644
--- /dev/null
+++ b/bench/Bench/Product.hs
@@ -0,0 +1,44 @@
+module Bench.Product
+  ( benchProduct
+  ) 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
+
+benchProduct :: Int -> Benchmark
+benchProduct k = bgroup (show (1 `shiftL` k :: Int))
+  [ bench "Bit/product"         $ nf (\x -> productBit    (randomVec Bit k) x)    (randomVec2 Bit k)
+  , bench "Bit/productShort"    $ nf (\x -> productBit    (randomVec Bit k) x)    (U.take 32 $ randomVec2 Bit k)
+  , bench "Bit/square"          $ nf (\x -> productBit    (randomVec Bit k) x)    (randomVec Bit k)
+  , bench "Bit.TS/product"      $ nf (\x -> productBitTS  (randomVec TS.Bit k) x) (randomVec2 TS.Bit k)
+  , bench "Bit.TS/productShort" $ nf (\x -> productBitTS  (randomVec TS.Bit k) x) (U.take 32 $ randomVec2 TS.Bit k)
+  , bench "Bit.TS/square"       $ nf (\x -> productBitTS  (randomVec TS.Bit k) x) (randomVec TS.Bit k)
+  ]
+
+productBit :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
+productBit xs ys = unF2Poly (toF2Poly xs * toF2Poly ys)
+
+productBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit -> U.Vector TS.Bit
+productBitTS xs ys = TS.unF2Poly (TS.toF2Poly xs * TS.toF2Poly ys)
diff --git a/bench/Bench/Reverse.hs b/bench/Bench/Reverse.hs
--- a/bench/Bench/Reverse.hs
+++ b/bench/Bench/Reverse.hs
@@ -24,21 +24,21 @@
 
 benchReverse :: Int -> Benchmark
 benchReverse k = bgroup (show (1 `shiftL` k :: Int))
-  [ bench "Bit/reverseInPlace"    $ nf reverseBit    (randomVec Bit k)
-  , bench "Bit/reverse"           $ nf reverseBit'   (randomVec Bit k)
-  , bench "Bit.TS/reverseInPlace" $ nf reverseBitTS  (randomVec TS.Bit k)
-  , bench "Bit.TS/reverse"        $ nf reverseBitTS' (randomVec TS.Bit k)
-  , bench "Vector"                $ nf reverseVector (randomVec id k)
+  [ 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)
   ]
 
 reverseBit :: U.Vector Bit -> U.Vector Bit
-reverseBit = U.modify reverseInPlace
+reverseBit = reverseBits
 
 reverseBit' :: U.Vector Bit -> U.Vector Bit
 reverseBit' = U.reverse
 
 reverseBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit
-reverseBitTS = U.modify TS.reverseInPlace
+reverseBitTS = TS.reverseBits
 
 reverseBitTS' :: U.Vector TS.Bit -> U.Vector TS.Bit
 reverseBitTS' = U.reverse
diff --git a/bench/Bench/Sum.hs b/bench/Bench/Sum.hs
new file mode 100644
--- /dev/null
+++ b/bench/Bench/Sum.hs
@@ -0,0 +1,42 @@
+module Bench.Sum
+  ( benchSum
+  ) 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
+
+benchSum :: Int -> Benchmark
+benchSum k = bgroup (show (1 `shiftL` k :: Int))
+  [ bench "Bit/add"    $ nf (\x -> sumBit    (randomVec Bit k) x)    (randomVec2 Bit k)
+  , bench "Bit/sum"    $ nf sum [(1 :: F2Poly) .. fromInteger (1 `shiftL` k)]
+  , bench "Bit.TS/add" $ nf (\x -> sumBitTS  (randomVec TS.Bit k) x) (randomVec2 TS.Bit k)
+  , bench "Bit.TS/sum" $ nf sum [(1 :: TS.F2Poly) .. fromInteger (1 `shiftL` k)]
+  ]
+
+sumBit :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
+sumBit xs ys = unF2Poly (toF2Poly xs + toF2Poly ys)
+
+sumBitTS :: U.Vector TS.Bit -> U.Vector TS.Bit -> U.Vector TS.Bit
+sumBitTS xs ys = TS.unF2Poly (TS.toF2Poly xs + TS.toF2Poly ys)
diff --git a/bench/Bench/Union.hs b/bench/Bench/Union.hs
--- a/bench/Bench/Union.hs
+++ b/bench/Bench/Union.hs
@@ -40,12 +40,12 @@
 
 benchUnion :: Int -> Benchmark
 benchUnion k = bgroup (show (1 `shiftL` k :: Int))
-  [ bench "Bit/zipBits"    $ nf (unionBit    (randomVec Bit k))    (randomVec2 Bit k)
-  , bench "Bit/zipWith"    $ nf (unionBit'   (randomVec Bit k))    (randomVec2 Bit k)
-  , bench "Bit.TS/zipBits" $ nf (unionBitTS  (randomVec TS.Bit k)) (randomVec2 TS.Bit k)
-  , bench "Bit.TS/zipWith" $ nf (unionBitTS' (randomVec TS.Bit k)) (randomVec2 TS.Bit k)
-  , bench "Vector"         $ nf (unionVector (randomVec id k))     (randomVec2 id k)
-  , bench "IntSet"         $ nf (unionIntSet (randomSet k))        (randomSet2 k)
+  [ 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)
   ]
 
 unionBit :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
diff --git a/bitvec.cabal b/bitvec.cabal
--- a/bitvec.cabal
+++ b/bitvec.cabal
@@ -1,5 +1,5 @@
 name: bitvec
-version: 1.0.0.1
+version: 1.0.1.0
 cabal-version: >=1.10
 build-type: Simple
 license: BSD3
@@ -9,10 +9,11 @@
 homepage: https://github.com/Bodigrim/bitvec
 synopsis: Space-efficient bit vectors
 description:
-  A newtype over 'Bool' with a better 'Vector' instance.
+  A newtype over 'Bool' with a better 'Vector' instance: 8x less memory, up to 1000x faster.
   .
-  The [vector](https://hackage.haskell.org/package/vector)
+  The <https://hackage.haskell.org/package/vector vector>
   package represents unboxed arrays of 'Bool'
+  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.__
@@ -20,7 +21,7 @@
   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 64x faster than 'Vector Bool'.
+  'Vector' 'Bit' is up to 1000x faster than 'Vector' 'Bool'.
   .
   === Thread safety
   .
@@ -33,11 +34,11 @@
   .
   === Similar packages
   .
-  * [bv](https://hackage.haskell.org/package/bv) and
-    [bv-little](https://hackage.haskell.org/package/bv-little)
+  * <https://hackage.haskell.org/package/bv bv> and
+    <https://hackage.haskell.org/package/bv-little bv-little>
     do not offer mutable vectors.
   .
-  * [array](https://hackage.haskell.org/package/array)
+  * <https://hackage.haskell.org/package/array array>
     is memory-efficient for 'Bool', but lacks
     a handy 'Vector' interface and is not thread-safe.
 
@@ -56,15 +57,23 @@
 
 flag bmi2
   description: Enable bmi2 instruction set
-  manual: False
   default: False
 
+flag integer-gmp
+  description: Use integer-gmp package for binary polynomials
+  default: True
+
+flag libgmp
+  description: Link against GMP library
+  default: True
+
 library
   exposed-modules:
     Data.Bit
     Data.Bit.ThreadSafe
   build-depends:
     base >=4.8 && <5,
+    deepseq,
     ghc-prim,
     primitive >=0.5,
     vector >=0.11
@@ -77,6 +86,9 @@
   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
@@ -87,9 +99,15 @@
     Data.Bit.Utils
   ghc-options: -O2 -Wall
   include-dirs: src
-  if (flag(bmi2)) && (impl(ghc >=8.4.1))
+  if flag(bmi2) && impl(ghc >=8.4.1)
     ghc-options: -mbmi2 -msse4.2
     cpp-options: -DBMI2_ENABLED
+  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
 
 test-suite bitvec-tests
   type: exitcode-stdio-1.0
@@ -133,9 +151,11 @@
     Bench.BitIndex
     Bench.Invert
     Bench.Intersection
+    Bench.Product
     Bench.RandomFlip
     Bench.RandomRead
     Bench.RandomWrite
     Bench.Reverse
+    Bench.Sum
     Bench.Union
   ghc-options: -O2 -Wall
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,3 +1,10 @@
+# 1.0.1.0
+
+* Implement arithmetic of binary polynomials.
+* Add `invertBits` and `reverseBits` functions.
+* Add `Num`, `Real`, `Integral`, `Fractional` and `NFData` instances.
+* Performance improvements.
+
 # 1.0.0.1
 
 * Performance improvements.
diff --git a/src/Data/Bit.hs b/src/Data/Bit.hs
--- a/src/Data/Bit.hs
+++ b/src/Data/Bit.hs
@@ -33,6 +33,8 @@
 
   -- * Immutable operations
   , zipBits
+  , invertBits
+  , reverseBits
   , bitIndex
   , nthBitIndex
   , countBits
@@ -46,20 +48,27 @@
   , cloneToWordsM
 
   -- * Mutable operations
-  , invertInPlace
   , zipInPlace
+  , invertInPlace
+  , reverseInPlace
   , selectBitsInPlace
   , excludeBitsInPlace
-  , reverseInPlace
+
+  -- * Binary polynomials
+  , F2Poly
+  , unF2Poly
+  , toF2Poly
   ) where
 
 import Prelude hiding (and, or)
 
 #ifndef BITVEC_THREADSAFE
+import Data.Bit.F2Poly
 import Data.Bit.Immutable
 import Data.Bit.Internal
 import Data.Bit.Mutable
 #else
+import Data.Bit.F2PolyTS
 import Data.Bit.ImmutableTS
 import Data.Bit.InternalTS
 import Data.Bit.MutableTS
diff --git a/src/Data/Bit/F2Poly.hs b/src/Data/Bit/F2Poly.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Bit/F2Poly.hs
@@ -0,0 +1,366 @@
+{-# LANGUAGE CPP                        #-}
+
+{-# LANGUAGE BangPatterns               #-}
+{-# LANGUAGE DeriveDataTypeable         #-}
+{-# LANGUAGE DeriveGeneric              #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE LambdaCase                 #-}
+{-# LANGUAGE MagicHash                  #-}
+{-# LANGUAGE RankNTypes                 #-}
+
+#ifndef BITVEC_THREADSAFE
+module Data.Bit.F2Poly
+#else
+module Data.Bit.F2PolyTS
+#endif
+  ( F2Poly
+  , unF2Poly
+  , toF2Poly
+  ) where
+
+import Control.DeepSeq
+import Control.Exception
+import Control.Monad
+import Control.Monad.ST
+#ifndef BITVEC_THREADSAFE
+import Data.Bit.Immutable
+import Data.Bit.Internal
+import Data.Bit.Mutable
+#else
+import Data.Bit.ImmutableTS
+import Data.Bit.InternalTS
+import Data.Bit.MutableTS
+#endif
+import Data.Bit.Utils
+import Data.Bits
+import Data.Coerce
+import Data.List hiding (dropWhileEnd)
+import Data.Typeable
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as MU
+import GHC.Generics
+
+#if UseIntegerGmp
+import Data.Primitive.ByteArray
+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'.
+--
+-- Polynomials are stored normalized, without leading zero coefficients.
+--
+-- '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)
+-- >>> 0b11 * 0b111 :: F2Poly
+-- F2Poly {unF2Poly = [1,0,0,1]}
+newtype F2Poly = F2Poly {
+  unF2Poly :: U.Vector Bit
+  -- ^ Convert 'F2Poly' to a vector of coefficients
+  -- (first element corresponds to a constant term).
+  }
+  deriving (Eq, Ord, Show, Typeable, Generic, NFData)
+
+-- | Make 'F2Poly' from a list of coefficients
+-- (first element corresponds to a constant term).
+toF2Poly :: U.Vector Bit -> F2Poly
+toF2Poly xs = F2Poly $ dropWhileEnd $ castFromWords $ cloneToWords xs
+
+-- | Addition and multiplication are evaluated modulo 2.
+--
+-- 'abs' = 'id' and 'signum' = 'const' 1.
+--
+-- 'fromInteger' converts a binary polynomial, encoded as 'Integer',
+-- to 'F2Poly' encoding.
+instance Num F2Poly where
+  (+) = coerce xorBits
+  (-) = coerce xorBits
+  negate = id
+  abs    = id
+  signum = const (F2Poly (U.singleton (Bit True)))
+  (*) = 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
+
+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
+
+instance Real F2Poly where
+  toRational = fromIntegral
+
+-- | 'toInteger' converts a binary polynomial, encoded as 'F2Poly',
+-- to 'Integer' encoding.
+instance Integral F2Poly where
+  toInteger = bitsToInteger . unF2Poly
+  quotRem (F2Poly xs) (F2Poly ys) = (F2Poly (dropWhileEnd qs), F2Poly (dropWhileEnd rs))
+    where
+      (qs, rs) = quotRemBits xs ys
+  rem = coerce ((dropWhileEnd .) . remBits)
+  divMod = quotRem
+  mod = rem
+
+xorBits
+  :: U.Vector Bit
+  -> U.Vector Bit
+  -> U.Vector Bit
+#if UseIntegerGmp
+-- GMP has platform-dependent ASM implementations for mpn_xor_n,
+-- which are impossible to beat by native Haskell.
+xorBits (BitVec _ 0 _) ys = ys
+xorBits xs (BitVec _ 0 _) = xs
+xorBits (BitVec 0 lx xarr) (BitVec 0 ly 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
+xorBits xs ys = dropWhileEnd $ runST $ do
+  let lx = U.length xs
+      ly = U.length ys
+      (shorterLen, longerLen, longer) = if lx >= ly then (ly, lx, xs) else (lx, ly, ys)
+  zs <- MU.replicate longerLen (Bit False)
+  forM_ [0, wordSize .. shorterLen - 1] $ \i ->
+    writeWord zs i (indexWord xs i `xor` indexWord ys i)
+  U.unsafeCopy (MU.drop shorterLen zs) (U.drop shorterLen longer)
+  U.unsafeFreeze zs
+
+-- | Must be >= wordSize.
+karatsubaThreshold :: Int
+karatsubaThreshold = 4096
+
+karatsuba :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
+karatsuba xs ys
+  | 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
+      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)
+    U.unsafeFreeze zs
+  where
+    lenXs = U.length xs
+    lenYs = U.length ys
+    lenZs = lenXs + lenYs - 1
+
+    m'    = ((lenXs `min` lenYs) + 1) `quot` 2
+    m     = if karatsubaThreshold < wordSize then m' else m' - modWordSize m'
+
+    xs0  = U.slice 0 m xs
+    xs1  = U.slice m (lenXs - m) xs
+    ys0  = U.slice 0 m ys
+    ys1  = U.slice m (lenYs - m) ys
+
+    xs01 = xorBits xs0 xs1
+    ys01 = xorBits ys0 ys1
+    zs0  = karatsuba xs0 ys0
+    zs2  = karatsuba xs1 ys1
+    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
+  where
+    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
+  | lenXs == 0 || lenYs == 0 = U.empty
+  | lenXs <= wordSize && lenYs <= wordSize = mulShortShort x0 y0
+  | lenYs <= wordSize                      = mulLongShort  xs y0
+  | lenXs <= wordSize                      = mulLongShort  ys x0
+  | otherwise = runST $ do
+    zs <- MU.replicate lenZs (Bit False)
+    forM_ [0 .. lenYs - 1] $ \k ->
+      MU.unsafeWrite zs k
+        (zipAndCountParityBits xs (U.unsafeSlice (lenYs - 1 - k) (k + 1) rys))
+    forM_ [lenYs .. lenZs - 1] $ \k ->
+      MU.unsafeWrite zs k
+        (zipAndCountParityBits (U.unsafeSlice (k - (lenYs - 1)) (lenXs + lenYs + 1 - k) xs) rys)
+    U.unsafeFreeze zs
+  where
+    lenXs = U.length xs
+    lenYs = U.length ys
+    lenZs = lenXs + lenYs - 1
+    rys   = reverseBits ys
+    x0 = indexWord xs 0 .&. loMask lenXs
+    y0 = indexWord ys 0 .&. loMask lenYs
+
+mulShortShort :: Word -> Word -> U.Vector Bit
+mulShortShort xs ys = runST $ do
+  zs <- MU.replicate lenZs (Bit False)
+  forM_ [0 .. lenYs - 1] $ \k -> do
+    let yk = rys `shiftR` (lenYs - 1 - k)
+        l  = (k + 1) `min` lenXs
+    MU.unsafeWrite zs k (fromIntegral $ popCount $ xs .&. yk .&. loMask l)
+  forM_ [lenYs .. lenZs - 1] $ \k -> do
+    let xk = xs `shiftR` (k - (lenYs - 1))
+        l  = (lenXs + lenYs + 1 - k) `min` lenYs
+    MU.unsafeWrite zs k (fromIntegral $ popCount $ xk .&. rys .&. loMask l)
+  U.unsafeFreeze zs
+  where
+    clzXs = countLeadingZeros xs
+    lenXs = wordSize - clzXs
+    clzYs = countLeadingZeros ys
+    lenYs = wordSize - clzYs
+    lenZs = lenXs + lenYs - 1
+    rys   = reverseWord (ys `shiftL` clzYs)
+
+mulLongShort :: U.Vector Bit -> Word -> U.Vector Bit
+mulLongShort xs ys = runST $ do
+  zs <- MU.replicate lenZs (Bit False)
+  forM_ [0 .. lenYs - 1] $ \k -> do
+    let yk = rys `shiftR` (lenYs - 1 - k)
+        l  = (k + 1) `min` lenXs
+    MU.unsafeWrite zs k (fromIntegral $ popCount $ x0 .&. yk .&. loMask l)
+  forM_ [lenYs .. lenZs - 1] $ \k -> do
+    let xk = indexWord xs (k - (lenYs - 1))
+        l  = (lenXs + lenYs + 1 - k) `min` lenYs
+    MU.unsafeWrite zs k (fromIntegral $ popCount $ xk .&. rys .&. loMask l)
+  U.unsafeFreeze zs
+  where
+    lenXs = U.length xs
+    clzYs = countLeadingZeros ys
+    lenYs = wordSize - clzYs
+    lenZs = lenXs + lenYs - 1
+    rys   = reverseWord (ys `shiftL` clzYs)
+    x0    = indexWord xs 0
+
+zipAndCountParityBits :: U.Vector Bit -> U.Vector Bit -> Bit
+zipAndCountParityBits xs ys
+  | nMod == 0 = fromIntegral $ popCnt
+  | otherwise = fromIntegral $ popCnt `xor` lastPopCnt
+  where
+    n = min (U.length xs) (U.length ys)
+    nMod = modWordSize n
+    ff i = indexWord xs i .&. indexWord ys i
+    popCnt = foldl' (\acc i -> acc `xor` popCount (ff i)) 0 [0, wordSize .. n - nMod - 1]
+    lastPopCnt = popCount (ff (n - nMod) .&. loMask nMod)
+
+sqrBits :: U.Vector Bit -> U.Vector Bit
+sqrBits xs = runST $ do
+    let lenXs = U.length xs
+    zs <- MU.replicate (lenXs `shiftL` 1) (Bit False)
+    forM_ [0, wordSize .. lenXs - 1] $ \i -> do
+      let (z0, z1) = sparseBits (indexWord xs i)
+      writeWord zs (i `shiftL` 1) z0
+      writeWord zs (i `shiftL` 1 + wordSize) z1
+    U.unsafeFreeze zs
+
+quotRemBits :: U.Vector Bit -> U.Vector Bit -> (U.Vector Bit, U.Vector Bit)
+quotRemBits xs ys
+  | U.null ys = throw DivideByZero
+  | U.length xs < U.length ys = (U.empty, xs)
+  | otherwise = runST $ do
+    let lenXs = U.length xs
+        lenYs = U.length ys
+        lenQs = lenXs - lenYs + 1
+    qs <- MU.replicate lenQs (Bit False)
+    rs <- MU.replicate lenXs (Bit False)
+    U.unsafeCopy rs xs
+    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do
+      Bit r <- MU.unsafeRead rs (lenYs - 1 + i)
+      when r $ do
+        MU.unsafeWrite qs i (Bit True)
+        zipInPlace xor ys (MU.drop i rs)
+    let rs' = MU.unsafeSlice 0 lenYs rs
+    (,) <$> U.unsafeFreeze qs <*> U.unsafeFreeze rs'
+
+remBits :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
+remBits xs ys
+  | U.null ys = throw DivideByZero
+  | U.length xs < U.length ys = xs
+  | otherwise = runST $ do
+    let lenXs = U.length xs
+        lenYs = U.length ys
+        lenQs = lenXs - lenYs + 1
+    rs <- MU.replicate lenXs (Bit False)
+    U.unsafeCopy rs xs
+    forM_ [lenQs - 1, lenQs - 2 .. 0] $ \i -> do
+      Bit r <- MU.unsafeRead rs (lenYs - 1 + i)
+      when r $ do
+        zipInPlace xor ys (MU.drop i rs)
+    let rs' = MU.unsafeSlice 0 lenYs rs
+    U.unsafeFreeze rs'
+
+dropWhileEnd
+  :: U.Vector Bit
+  -> U.Vector Bit
+dropWhileEnd xs = U.unsafeSlice 0 (go (U.length xs)) xs
+  where
+    go n
+      | n < wordSize = wordSize - countLeadingZeros (indexWord xs 0 .&. loMask n)
+      | otherwise    = case indexWord xs (n - wordSize) of
+        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
+    !(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`)
+  $ [0..]
+
+bitsToInteger :: U.Vector Bit -> Integer
+bitsToInteger = U.ifoldl' (\acc i (Bit b) -> if b then acc `setBit` i else acc) 0
+
+#endif
diff --git a/src/Data/Bit/F2PolyTS.hs b/src/Data/Bit/F2PolyTS.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Bit/F2PolyTS.hs
@@ -0,0 +1,4 @@
+{-# LANGUAGE CPP #-}
+
+#define BITVEC_THREADSAFE
+#include "Data/Bit/F2Poly.hs"
diff --git a/src/Data/Bit/Gmp.hs b/src/Data/Bit/Gmp.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Bit/Gmp.hs
@@ -0,0 +1,146 @@
+{-# 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
@@ -14,18 +14,21 @@
   , cloneToWords
 
   , zipBits
-
+  , invertBits
   , selectBits
   , excludeBits
-  , bitIndex
+  , reverseBits
 
+  , bitIndex
   , nthBitIndex
   , countBits
   , listBits
   ) where
 
+import Control.Monad
 import Control.Monad.ST
 import Data.Bits
+import Data.Bit.Gmp
 #ifndef BITVEC_THREADSAFE
 import Data.Bit.Internal
 import Data.Bit.Mutable
@@ -39,8 +42,17 @@
 import qualified Data.Vector.Primitive as P
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
-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
+#endif
+
 -- | Cast a vector of words to a vector of bits.
 -- Cf. 'Data.Bit.castFromWordsM'.
 --
@@ -48,7 +60,8 @@
 -- [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]
 castFromWords :: U.Vector Word -> U.Vector Bit
 castFromWords ws = BitVec (mulWordSize off) (mulWordSize len) arr
-  where P.Vector off len arr = unsafeCoerce ws
+  where
+    P.Vector off len arr = toPrimVector ws
 
 -- | Try to cast a vector of bits to a vector of words.
 -- It succeeds if a vector of bits is aligned.
@@ -58,11 +71,11 @@
 -- prop> castToWords (castFromWords v) == Just v
 castToWords :: U.Vector Bit -> Maybe (U.Vector Word)
 castToWords (BitVec s n ws)
-  | aligned s, aligned n = Just $ unsafeCoerce $ P.Vector (divWordSize s)
-                                                          (divWordSize n)
-                                                          ws
+  | aligned s, aligned n =
+    Just $ fromPrimVector $ P.Vector (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.cloneToWordsM'.
@@ -77,8 +90,13 @@
 {-# INLINE cloneToWords #-}
 
 -- | Zip two vectors with the given function.
--- Similar to 'Data.Vector.Unboxed.zipWith', but up to 16x faster.
+-- Similar to 'Data.Vector.Unboxed.zipWith',
+-- but up to 1000x (!) faster.
 --
+-- For sufficiently dense sets, represented as bitmaps,
+-- 'zipBits' is up to 32x faster than
+-- 'Data.IntSet.union', 'Data.IntSet.intersection', etc.
+--
 -- >>> import Data.Bits
 -- >>> zipBits (.&.) (read "[1,1,0]") (read "[0,1,1]") -- intersection
 -- [0,1,0]
@@ -93,11 +111,67 @@
   -> U.Vector Bit
   -> U.Vector Bit
   -> U.Vector Bit
-zipBits f xs ys | U.length xs >= U.length ys = zs
-                | otherwise                  = U.slice 0 (U.length xs) zs
-  where zs = U.modify (zipInPlace f xs) ys
+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
+        w = nWords l
+        b = w `shiftL` GMP_LIMB_SHIFT
+    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)
+    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
 {-# INLINE zipBits #-}
 
+-- | Invert (flip) all bits.
+--
+-- >>> invertBits (read "[0,1,0,1,0]")
+-- [1,0,1,0,1]
+invertBits
+  :: U.Vector Bit
+  -> U.Vector Bit
+invertBits (BitVec _ 0 _) = U.empty
+#if UseLibGmp
+invertBits (BitVec 0 l arg) = runST $ do
+  let w = nWords l
+  brr <- newByteArray (w `shiftL` GMP_LIMB_SHIFT)
+  mpnCom brr arg w
+  BitVec 0 l <$> unsafeFreezeByteArray brr
+#endif
+invertBits xs = runST $ do
+  let n = U.length xs
+  ys <- MU.new n
+  forM_ [0, wordSize .. n - 1] $ \i ->
+    writeWord ys i (complement (indexWord xs i))
+  U.unsafeFreeze ys
+
 -- | For each set bit of the first argument, deposit
 -- the corresponding bit of the second argument
 -- to the result. Similar to the parallel deposit instruction (PDEP).
@@ -132,6 +206,26 @@
   n   <- excludeBitsInPlace is xs1
   U.unsafeFreeze (MU.take n xs1)
 
+-- | Reverse the order of bits.
+--
+-- >>> reverseBits (read "[1,1,0,1,0]")
+-- [0,1,0,1,1]
+reverseBits :: U.Vector Bit -> U.Vector Bit
+reverseBits xs = runST $ do
+  let n    = U.length xs
+  ys <- MU.new n
+
+  forM_ [0, wordSize .. n - wordSize] $ \i ->
+    writeWord ys (n - i - wordSize) (reverseWord (indexWord xs i))
+
+  let nMod = modWordSize n
+  when (nMod /= 0) $ do
+    let x = indexWord xs (mulWordSize (divWordSize n))
+    y <- readWord ys 0
+    writeWord ys 0 (meld nMod (reversePartialWord nMod x) y)
+
+  U.unsafeFreeze ys
+
 clipLoBits :: Bit -> Int -> Word -> Word
 clipLoBits (Bit True ) k w = w `unsafeShiftR` k
 clipLoBits (Bit False) k w = (w `unsafeShiftR` k) .|. hiMask (wordSize - k)
@@ -320,6 +414,10 @@
 -- for <https://en.wikipedia.org/wiki/Succinct_data_structure succinct dictionaries>.
 countBits :: U.Vector Bit -> Int
 countBits (BitVec _ 0 _)                      = 0
+#if UseLibGmp
+countBits (BitVec 0 len arr) | modWordSize len == 0 =
+  fromIntegral (mpnPopcount arr (divWordSize len))
+#endif
 countBits (BitVec off len arr) | offBits == 0 = case modWordSize len of
   0    -> countBitsInWords (P.Vector offWords lWords arr)
   nMod -> countBitsInWords (P.Vector offWords (lWords - 1) arr)
diff --git a/src/Data/Bit/Internal.hs b/src/Data/Bit/Internal.hs
--- a/src/Data/Bit/Internal.hs
+++ b/src/Data/Bit/Internal.hs
@@ -2,6 +2,7 @@
 
 {-# LANGUAGE BangPatterns               #-}
 {-# LANGUAGE DeriveDataTypeable         #-}
+{-# LANGUAGE DeriveGeneric              #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE LambdaCase                 #-}
 {-# LANGUAGE MagicHash                  #-}
@@ -28,15 +29,20 @@
 
 #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 qualified Data.Vector.Unboxed as U
+import GHC.Generics
 
 #ifdef BITVEC_THREADSAFE
 import GHC.Exts
@@ -50,7 +56,7 @@
 -- 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)
+  deriving (Bounded, Enum, Eq, Ord, FiniteBits, Bits, Typeable, Generic, NFData)
 #else
 -- | A newtype wrapper with a custom instance
 -- of "Data.Vector.Unboxed", which packs booleans
@@ -59,9 +65,47 @@
 -- 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)
+  deriving (Bounded, Enum, Eq, Ord, FiniteBits, Bits, Typeable, Generic, NFData)
 #endif
 
+-- | There is only one lawful 'Num' instance possible
+-- with '+' = 'xor' and
+-- 'fromInteger' = 'Bit' . 'odd'.
+instance Num Bit where
+  Bit a * Bit b = Bit (a && b)
+  Bit a + Bit b = Bit (a /= b)
+  Bit a - Bit b = Bit (a /= b)
+  negate = id
+  abs    = id
+  signum = id
+  fromInteger = Bit . odd
+
+instance Real Bit where
+  toRational (Bit False) = 0
+  toRational (Bit True)  = 1
+
+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
+
+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
+
 instance Show Bit where
   showsPrec _ (Bit False) = showString "0"
   showsPrec _ (Bit True ) = showString "1"
@@ -93,15 +137,14 @@
 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 zero-padded.
+-- | 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 off len' arr) i' = word .&. msk
+indexWord (BitVec off len' arr) i' = word
  where
   len    = off + len'
   i      = off + i'
   nMod   = modWordSize i
   loIx   = divWordSize i
-  msk    = if len - i >= wordSize then complement 0 else loMask (len - i)
   loWord = indexByteArray arr loIx
   hiWord = indexByteArray arr (loIx + 1)
 
@@ -113,17 +156,16 @@
         (loWord `unsafeShiftR` nMod)
           .|. (hiWord `unsafeShiftL` (wordSize - nMod))
 
--- | 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 zero-padded.
+-- | 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 off len' arr) i' = do
   let len  = off + len'
       i    = off + i'
       nMod = modWordSize i
       loIx = divWordSize i
-      msk  = if len - i >= wordSize then complement 0 else loMask (len - i)
   loWord <- readByteArray arr loIx
 
-  word   <- if nMod == 0
+  if nMod == 0
     then pure loWord
     else if loIx == divWordSize (len - 1)
       then pure (loWord `unsafeShiftR` nMod)
@@ -132,8 +174,9 @@
         pure
           $   (loWord `unsafeShiftR` nMod)
           .|. (hiWord `unsafeShiftL` (wordSize - nMod))
-
-  pure $ word .&. msk
+#if __GLASGOW_HASKELL__ >= 800
+{-# SPECIALIZE readWord :: U.MVector s Bit -> Int -> ST s Word #-}
+#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.
 writeWord :: PrimMonad m => U.MVector (PrimState m) Bit -> Int -> Word -> m ()
@@ -173,6 +216,9 @@
         writeByteArray arr (loIx + 1)
           $   (hiWord .&. hiMask nMod)
           .|. (x `unsafeShiftR` (wordSize - nMod))
+#if __GLASGOW_HASKELL__ >= 800
+{-# SPECIALIZE writeWord :: U.MVector s Bit -> Int -> Word -> ST s () #-}
+#endif
 
 instance MV.MVector U.MVector Bit where
   {-# INLINE basicInitialize #-}
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
@@ -22,7 +22,9 @@
   , reverseInPlace
   ) where
 
+import Control.Monad
 import Control.Monad.Primitive
+import Control.Monad.ST
 #ifndef BITVEC_THREADSAFE
 import Data.Bit.Internal
 #else
@@ -85,37 +87,27 @@
   -> Vector Bit
   -> MVector (PrimState m) Bit
   -> m ()
-zipInPlace f xs ys = loop 0
- where
-  !n = min (U.length xs) (MU.length ys)
-  loop !i
-    | i >= n = pure ()
-    | otherwise = do
-      let x = indexWord xs i
-      y <- readWord ys i
-      writeWord ys i (f x y)
-      loop (i + wordSize)
+zipInPlace f xs ys = do
+  let n = min (U.length xs) (MU.length ys)
+  forM_ [0, wordSize .. n - 1] $ \i -> do
+    let x = indexWord xs i
+    y <- readWord ys i
+    writeWord ys i (f x y)
 {-# INLINE zipInPlace #-}
 
 -- | Invert (flip) all bits in-place.
 --
--- Combine with 'Data.Vector.Unboxed.modify'
--- or simply resort to 'Data.Vector.Unboxed.map' 'Data.Bits.complement'
--- to operate on immutable vectors.
---
 -- >>> Data.Vector.Unboxed.modify invertInPlace (read "[0,1,0,1,0]")
 -- [1,0,1,0,1]
 invertInPlace :: PrimMonad m => U.MVector (PrimState m) Bit -> m ()
-invertInPlace xs = loop 0
- where
-  !n = MU.length xs
-  loop !i
-    | i >= n = pure ()
-    | otherwise = do
-      x <- readWord xs i
-      writeWord xs i (complement x)
-      loop (i + wordSize)
-{-# INLINE invertInPlace #-}
+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.
@@ -152,10 +144,6 @@
 
 -- | Reverse the order of bits in-place.
 --
--- Combine with 'Data.Vector.Unboxed.modify'
--- or simply resort to 'Data.Vector.Unboxed.reverse'
--- to operate on immutable vectors.
---
 -- >>> Data.Vector.Unboxed.modify reverseInPlace (read "[1,1,0,1,0]")
 -- [0,1,0,1,1]
 reverseInPlace :: PrimMonad m => U.MVector (PrimState m) Bit -> m ()
@@ -191,3 +179,6 @@
     !j  = len - i
     !i' = i + wordSize
     !j' = j - wordSize
+#if __GLASGOW_HASKELL__ >= 800
+{-# SPECIALIZE reverseInPlace :: U.MVector s Bit -> ST s () #-}
+#endif
diff --git a/src/Data/Bit/Select1.hs b/src/Data/Bit/Select1.hs
--- a/src/Data/Bit/Select1.hs
+++ b/src/Data/Bit/Select1.hs
@@ -141,7 +141,9 @@
 select1 :: Word -> Int -> Int
 #if WORD_SIZE_IN_BITS == 64
 select1 w i = fromIntegral $ select1Word64 (fromIntegral w) (fromIntegral i)
-#else
+#elif WORD_SIZE_IN_BITS == 32
 select1 w i = fromIntegral $ select1Word32 (fromIntegral w) (fromIntegral i)
+#else
+#error unsupported WORD_SIZE_IN_BITS config
 #endif
 {-# INLINE select1 #-}
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
@@ -18,11 +18,17 @@
   , ffs
   , loMask
   , hiMask
+  , sparseBits
+  , fromPrimVector
+  , toPrimVector
   ) where
 
 #include "MachDeps.h"
 
 import Data.Bits
+import qualified Data.Vector.Primitive as P
+import qualified Data.Vector.Unboxed as U
+import Unsafe.Coerce
 
 -- |The number of bits in a 'Word'.  A handy constant to have around when defining 'Word'-based bulk operations on bit vectors.
 wordSize :: Int
@@ -97,15 +103,17 @@
   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)
-#else
+#elif WORD_SIZE_IN_BITS == 32
 reverseWord :: Word -> Word
 reverseWord x0 = x5
  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)
+  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
@@ -127,8 +135,65 @@
                          (if testBit x i then setBit y ct else y)
     | otherwise = loop (i + 1) ct y
 
+#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)
+  where
+    x = sparseBitsInternal (w .&. loMask 32)
+    y = sparseBitsInternal (w `shiftR` 32)
+
+sparseBitsInternal :: Word -> Word
+sparseBitsInternal x = x4
+  where
+    t  = (x  `xor` (x  `shiftR` 16)) .&. 0x00000000ffff0000
+    x0 = x  `xor` (t  `xor` (t  `shiftL` 16));
+
+    t0 = (x0 `xor` (x0 `shiftR` 8)) .&. 0x0000ff000000ff00;
+    x1 = x0 `xor` (t0 `xor` (t0 `shiftL` 8));
+    t1 = (x1 `xor` (x1 `shiftR` 4)) .&. 0x00f000f000f000f0;
+    x2 = x1 `xor` (t1 `xor` (t1 `shiftL` 4));
+    t2 = (x2 `xor` (x2 `shiftR` 2)) .&. 0x0c0c0c0c0c0c0c0c;
+    x3 = x2 `xor` (t2 `xor` (t2 `shiftL` 2));
+    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)
+  where
+    x = sparseBitsInternal (w .&. loMask 16)
+    y = sparseBitsInternal (w `shiftR` 16)
+
+sparseBitsInternal :: Word -> Word
+sparseBitsInternal x0 = x4
+  where
+    t0 = (x0 `xor` (x0 `shiftR` 8)) .&. 0x0000ff00;
+    x1 = x0 `xor` (t0 `xor` (t0 `shiftL` 8));
+    t1 = (x1 `xor` (x1 `shiftR` 4)) .&. 0x00f000f0;
+    x2 = x1 `xor` (t1 `xor` (t1 `shiftL` 4));
+    t2 = (x2 `xor` (x2 `shiftR` 2)) .&. 0x0c0c0c0c;
+    x3 = x2 `xor` (t2 `xor` (t2 `shiftL` 2));
+    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 `shiftL` n - 1
 
 hiMask :: Int -> Word
 hiMask n = complement (1 `shiftL` n - 1)
+
+fromPrimVector :: P.Vector Word -> U.Vector Word
+fromPrimVector = unsafeCoerce
+
+toPrimVector :: U.Vector Word -> P.Vector Word
+toPrimVector = unsafeCoerce
diff --git a/test/Main.hs b/test/Main.hs
--- a/test/Main.hs
+++ b/test/Main.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE CPP #-}
+
 module Main where
 
 import Data.Bit
@@ -14,11 +16,32 @@
 main :: IO ()
 main = defaultMain $ testGroup
   "All"
-  [showReadTests, mvectorTests, TS.mvectorTests, setOpTests, vectorTests]
+  [lawsTests, f2polyTests, mvectorTests, TS.mvectorTests, setOpTests, vectorTests]
 
-showReadTests :: TestTree
-showReadTests =
-  testGroup "Show/Read"
-    $ map (uncurry testProperty)
-    $ lawsProperties
-    $ showReadLaws (Proxy :: Proxy Bit)
+lawsTests :: TestTree
+lawsTests = testGroup "Laws"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [ bitsLaws        (Proxy :: Proxy Bit)
+  , eqLaws          (Proxy :: Proxy Bit)
+  , ordLaws         (Proxy :: Proxy Bit)
+  , 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"
+  $ map (uncurry testProperty)
+  $ concatMap lawsProperties
+  [ showLaws        (Proxy :: Proxy F2Poly)
+#if MIN_VERSION_quickcheck_classes(0,6,3)
+  , numLaws         (Proxy :: Proxy F2Poly)
+#endif
+  , integralLaws    (Proxy :: Proxy F2Poly)
+  ]
+
diff --git a/test/Support.hs b/test/Support.hs
--- a/test/Support.hs
+++ b/test/Support.hs
@@ -38,7 +38,17 @@
   function f = functionMap TS.unBit TS.Bit f
 
 instance (Arbitrary a, U.Unbox a) => Arbitrary (U.Vector a) where
-  arbitrary = V.new <$> arbitrary
+  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)
+    ]
+
+instance Arbitrary F2Poly where
+  arbitrary = toF2Poly <$> arbitrary
+  shrink v = toF2Poly <$> shrink (unF2Poly v)
 
 instance (Show (v a), V.Vector v a) => Show (N.New v a) where
   showsPrec p = showsPrec p . V.new
diff --git a/test/Tests/MVector.hs b/test/Tests/MVector.hs
--- a/test/Tests/MVector.hs
+++ b/test/Tests/MVector.hs
@@ -37,7 +37,6 @@
     [ testProperty "cloneFromWords" prop_cloneFromWords_def
     , testProperty "cloneToWords"   prop_cloneToWords_def
     ]
-  , testProperty "reverseInPlace" prop_reverseInPlace_def
   , testGroup "MVector laws"
   $ map (uncurry testProperty)
   $ lawsProperties
@@ -223,9 +222,3 @@
 prop_cloneToWords_def xs =
   runST (N.run xs >>= cloneToWordsM >>= V.unsafeFreeze)
     === cloneToWords (V.new xs)
-
-prop_reverseInPlace_def :: N.New B.Vector Bit -> Property
-prop_reverseInPlace_def xs =
-  runST (N.run xs >>= \v -> reverseInPlace v >> V.unsafeFreeze v)
-    === B.reverse (V.new xs)
-
diff --git a/test/Tests/SetOps.hs b/test/Tests/SetOps.hs
--- a/test/Tests/SetOps.hs
+++ b/test/Tests/SetOps.hs
@@ -1,10 +1,11 @@
+{-# LANGUAGE RankNTypes #-}
+
 module Tests.SetOps where
 
 import Support ()
 
 import Data.Bit
 import Data.Bits
-import Data.List.NonEmpty (NonEmpty(..))
 import qualified Data.Vector.Unboxed as U
 import Test.Tasty
 import Test.Tasty.QuickCheck hiding ((.&.))
@@ -12,78 +13,76 @@
 setOpTests :: TestTree
 setOpTests = testGroup
   "Set operations"
-  [ testProperty "union"         prop_union_def
-  , testProperty "intersection"  prop_intersection_def
-  , testProperty "difference"    prop_difference_def
-  , testProperty "symDiff"       prop_symDiff_def
-  -- , testProperty "unions"        prop_unions_def
-  -- , testProperty "intersections" prop_unions_def
-  , testProperty "invert"        prop_invert_def
-  , testProperty "select"        prop_select_def
-  , testProperty "exclude"       prop_exclude_def
-  , testProperty "selectBits"    prop_selectBits_def
-  , testProperty "excludeBits"   prop_excludeBits_def
-  , testProperty "countBits"     prop_countBits_def
+  [ 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 "reverseBits"     prop_reverseBits
+  , testProperty "reverseInPlace"  prop_reverseInPlace
+  , testProperty "select"          prop_select_def
+  , testProperty "exclude"         prop_exclude_def
+  , testProperty "selectBits"      prop_selectBits_def
+  , testProperty "excludeBits"     prop_excludeBits_def
+  , testProperty "countBits"       prop_countBits_def
   ]
 
-union :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
-union = zipBits (.|.)
+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
 
 prop_union_def :: U.Vector Bit -> U.Vector Bit -> Property
 prop_union_def xs ys =
-  U.toList (union xs ys) === zipWith (.|.) (U.toList xs) (U.toList ys)
-
-intersection :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
-intersection = zipBits (.&.)
+  zipBits (.|.) xs ys === U.zipWith (.|.) xs ys
 
 prop_intersection_def :: U.Vector Bit -> U.Vector Bit -> Property
 prop_intersection_def xs ys =
-  U.toList (intersection xs ys) === zipWith (.&.) (U.toList xs) (U.toList ys)
-
-difference :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
-difference = zipBits (\a b -> a .&. complement b)
+  zipBits (.&.) xs ys === U.zipWith (.&.) xs ys
 
 prop_difference_def :: U.Vector Bit -> U.Vector Bit -> Property
-prop_difference_def xs ys = U.toList (difference xs ys)
-  === zipWith diff (U.toList xs) (U.toList ys)
-  where diff x y = x .&. complement y
-
-symDiff :: U.Vector Bit -> U.Vector Bit -> U.Vector Bit
-symDiff = zipBits xor
+prop_difference_def xs ys =
+  zipBits diff xs ys === U.zipWith diff xs ys
+  where
+    diff x y = x .&. complement y
 
 prop_symDiff_def :: U.Vector Bit -> U.Vector Bit -> Property
 prop_symDiff_def xs ys =
-  U.toList (symDiff xs ys) === zipWith xor (U.toList xs) (U.toList ys)
+  zipBits xor xs ys === U.zipWith xor xs ys
 
-unions :: NonEmpty (U.Vector Bit) -> U.Vector Bit
-unions (x :| xs) = U.slice 0 l $ U.modify (go xs) x
- where
-  l = minimum $ fmap U.length (x :| xs)
-  go []       _   = pure ()
-  go (y : ys) acc = do
-    zipInPlace (.|.) y acc
-    go ys acc
+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
+  where
+    f = curry $ applyFun fun
 
-prop_unions_def :: U.Vector Bit -> [U.Vector Bit] -> Property
-prop_unions_def xs xss = unions (xs :| xss) === foldr union xs xss
+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)
+  where
+    f = curry $ applyFun fun
 
-intersections :: NonEmpty (U.Vector Bit) -> U.Vector Bit
-intersections (x :| xs) = U.slice 0 l $ U.modify (go xs) x
- where
-  l = minimum $ fmap U.length (x :| xs)
-  go []       _   = pure ()
-  go (y : ys) acc = do
-    zipInPlace (.&.) y acc
-    go ys acc
+prop_invertBits :: U.Vector Bit -> Property
+prop_invertBits xs =
+  U.map complement xs === invertBits xs
 
-prop_intersections_def :: U.Vector Bit -> [U.Vector Bit] -> Property
-prop_intersections_def xs xss =
-  intersections (xs :| xss) === foldr intersection xs xss
+prop_invertBitsWords :: U.Vector Word -> Property
+prop_invertBitsWords ws =
+  U.map complement xs === invertBits xs
+  where
+    xs = castFromWords ws
 
-prop_invert_def :: U.Vector Bit -> Bool
-prop_invert_def xs =
-  U.toList (U.modify invertInPlace xs) == map complement (U.toList xs)
+prop_invertInPlace :: U.Vector Bit -> Property
+prop_invertInPlace xs =
+  U.map complement xs === U.modify invertInPlace xs
 
+prop_reverseBits :: U.Vector Bit -> Property
+prop_reverseBits xs =
+  U.reverse xs === reverseBits xs
+
+prop_reverseInPlace :: U.Vector Bit -> Property
+prop_reverseInPlace xs =
+  U.reverse xs === U.modify reverseInPlace xs
+
 select :: U.Unbox a => U.Vector Bit -> U.Vector a -> [a]
 select mask ws = U.toList (U.map snd (U.filter (unBit . fst) (U.zip mask ws)))
 
@@ -107,3 +106,27 @@
 
 prop_countBits_def :: U.Vector Bit -> Bool
 prop_countBits_def xs = countBits xs == U.length (selectBits xs xs)
+
+-------------------------------------------------------------------------------
+
+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
+  (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
+  (Bit False, Bit False, Bit True,  Bit True)  -> \x _ -> x
+
+  (Bit False, Bit True,  Bit False, Bit False) -> \x y -> complement x .&. y
+  (Bit False, Bit True,  Bit False, Bit True)  -> \_ y -> y
+  (Bit False, Bit True,  Bit True,  Bit False) -> \x y -> x `xor` y
+  (Bit False, Bit True,  Bit True,  Bit True)  -> \x y -> x .|. y
+
+  (Bit True,  Bit False, Bit False, Bit False) -> \x y -> complement (x .|. y)
+  (Bit True,  Bit False, Bit False, Bit True)  -> \x y -> complement (x `xor` y)
+  (Bit True,  Bit False, Bit True,  Bit False) -> \_ y -> complement y
+  (Bit True,  Bit False, Bit True,  Bit True)  -> \x y -> x .|. complement y
+
+  (Bit True,  Bit True,  Bit False, Bit False) -> \x _ -> complement x
+  (Bit True,  Bit True,  Bit False, Bit True)  -> \x y -> complement x .|. y
+  (Bit True,  Bit True,  Bit True,  Bit False) -> \x y -> complement (x .&. y)
+  (Bit True,  Bit True,  Bit True,  Bit True)  -> \_ _ -> complement zeroBits