diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -133,6 +133,10 @@
 using `chimera` and `memoize` packages.
 
 ```haskell
+#!/usr/bin/env cabal
+{- cabal:
+build-depends: base, chimera, memoize, time
+-}
 {-# LANGUAGE TypeApplications #-}
 import Data.Chimera
 import Data.Function.Memoize
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -1,85 +1,12 @@
-{-# LANGUAGE CPP #-}
-
 module Main where
 
-import Control.Monad.State (evalState, put, get)
-import Data.Bits
-import Data.Chimera
 import Test.Tasty.Bench
-import Test.Tasty.Patterns.Printer
-import System.Random
 
-#ifdef MIN_VERSION_ral
-import qualified Data.RAList as RAL
-#endif
-
-sizes :: Num a => [a]
-sizes = [7, 8, 9, 10]
+import Memoize
+import Read
 
 main :: IO ()
-main = defaultMain $ (: []) $ mapLeafBenchmarks addCompare $ bgroup "read"
-  [ bgroup chimeraBenchName (map benchReadChimera sizes)
-  , bgroup "List"           (map benchReadList    sizes)
-#ifdef MIN_VERSION_ral
-  , bgroup "RAL"            (map benchReadRAL     sizes)
-#endif
+main = defaultMain
+  [ readBenchmark
+  , memoizeBenchmark
   ]
-
-chimeraBenchName :: String
-chimeraBenchName = "Chimera"
-
-addCompare :: ([String] -> Benchmark -> Benchmark)
-addCompare (size : name : path)
-  | name /= chimeraBenchName
-  = bcompare (printAwkExpr (locateBenchmark (size : chimeraBenchName : path)))
-addCompare _ = id
-
-randomChimera :: UChimera Int
-randomChimera = flip evalState (mkStdGen 42) $ tabulateM $ const $ do
-  g <- get
-  let (x, g') = random g
-  put g'
-  pure x
-
-randomList :: [Int]
-randomList = randoms (mkStdGen 42)
-
-#ifdef MIN_VERSION_ral
-randomRAL :: RAL.RAList Int
-randomRAL = RAL.fromList $ take (1 `shiftL` (maximum sizes)) $ randoms (mkStdGen 42)
-#endif
-
-randomIndicesWord :: [Word]
-randomIndicesWord = randoms (mkStdGen 42)
-
-randomIndicesInt :: [Int]
-randomIndicesInt = randoms (mkStdGen 42)
-
-benchReadChimera :: Int -> Benchmark
-benchReadChimera k
-  = bench (show n)
-  $ nf (sum . map (index randomChimera))
-  $ map (.&. (n - 1))
-  $ take (fromIntegral n) randomIndicesWord
-  where
-    n = 1 `shiftL` k
-
-benchReadList :: Int -> Benchmark
-benchReadList k
-  = bench (show n)
-  $ nf (sum . map (randomList !!))
-  $ map (.&. (n - 1))
-  $ take n randomIndicesInt
-  where
-    n = 1 `shiftL` k
-
-#ifdef MIN_VERSION_ral
-benchReadRAL :: Int -> Benchmark
-benchReadRAL k
-  = bench (show n)
-  $ nf (sum . map (randomRAL RAL.!))
-  $ map (.&. (n - 1))
-  $ take n randomIndicesInt
-  where
-    n = 1 `shiftL` k
-#endif
diff --git a/bench/Memoize.hs b/bench/Memoize.hs
new file mode 100644
--- /dev/null
+++ b/bench/Memoize.hs
@@ -0,0 +1,44 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE TypeApplications #-}
+
+module Memoize
+  ( memoizeBenchmark
+  ) where
+
+import Data.Bits
+import Data.Chimera
+import Data.Foldable
+import Data.Function
+import qualified Data.Vector.Unboxed as U
+import Test.Tasty.Bench
+
+memoizeBenchmark :: Benchmark
+memoizeBenchmark = bgroup "memoize"
+  [ bgroup "memoizeFix" $ memoizeFixBenchmark memoizeFix
+  , bgroup "memoizeFix unboxed" $ memoizeFixBenchmark (index . tabulateFix @U.Vector)
+  , bgroup "fix memoize" $ memoizeFixBenchmark (fix . (memoize .))
+  ]
+
+memoizeFixBenchmark :: (forall a. U.Unbox a => ((Word -> a) -> Word -> a) -> Word -> a) -> [Benchmark]
+memoizeFixBenchmark fixer =
+  [ bench "isOdd" $ nf (\f -> let isOdd = fixer f in
+      foldl' (\acc n -> xor acc (isOdd n)) False [0..10000]) isOddF
+  , bench "isPrime" $ nf (\f -> let isPrime = fixer f in
+      foldl' (\acc n -> xor acc (isPrime n)) False [0..10000]) isPrimeF
+  , bench "fibo" $ nf (\f -> let fibo = fixer f in
+      foldl' (\acc n -> acc + fibo n) 0 [0..10000]) fiboF
+  , bench "collatz" $ nf (\f -> let collatz = fixer f in
+      foldl' (\acc n -> acc + collatz n) 0 [0..1000]) collatzF
+  ]
+
+isOddF :: (Word -> Bool) -> Word -> Bool
+isOddF f n = n /= 0 && not (f (n - 1))
+
+isPrimeF :: (Word -> Bool) -> Word -> Bool
+isPrimeF f n = n > 1 && and [ n `rem` d /= 0 | d <- [2 .. floor (sqrt (fromIntegral n :: Double))], f d]
+
+fiboF :: (Word -> Word) -> Word -> Word
+fiboF f n = if n < 2 then fromIntegral n else f (n - 1) + f (n - 2)
+
+collatzF :: (Word -> Word) -> Word -> Word
+collatzF f n = if n <= 1 then 0 else 1 + f (if even n then n `quot` 2 else 3 * n + 1)
diff --git a/bench/Read.hs b/bench/Read.hs
new file mode 100644
--- /dev/null
+++ b/bench/Read.hs
@@ -0,0 +1,87 @@
+{-# LANGUAGE CPP #-}
+
+module Read
+  ( readBenchmark
+  ) where
+
+import Control.Monad.State (evalState, put, get)
+import Data.Bits
+import Data.Chimera
+import Test.Tasty.Bench
+import Test.Tasty.Patterns.Printer
+import System.Random
+
+#ifdef MIN_VERSION_ral
+import qualified Data.RAList as RAL
+#endif
+
+sizes :: Num a => [a]
+sizes = [7, 8, 9, 10]
+
+readBenchmark :: Benchmark
+readBenchmark = mapLeafBenchmarks addCompare $ bgroup "read"
+  [ bgroup chimeraBenchName (map benchReadChimera sizes)
+  , bgroup "List"           (map benchReadList    sizes)
+#ifdef MIN_VERSION_ral
+  , bgroup "RAL"            (map benchReadRAL     sizes)
+#endif
+  ]
+
+chimeraBenchName :: String
+chimeraBenchName = "Chimera"
+
+addCompare :: ([String] -> Benchmark -> Benchmark)
+addCompare (size : name : path)
+  | name /= chimeraBenchName
+  = bcompare (printAwkExpr (locateBenchmark (size : chimeraBenchName : path)))
+addCompare _ = id
+
+randomChimera :: UChimera Int
+randomChimera = flip evalState (mkStdGen 42) $ tabulateM $ const $ do
+  g <- get
+  let (x, g') = random g
+  put g'
+  pure x
+
+randomList :: [Int]
+randomList = randoms (mkStdGen 42)
+
+#ifdef MIN_VERSION_ral
+randomRAL :: RAL.RAList Int
+randomRAL = RAL.fromList $ take (1 `shiftL` (maximum sizes)) $ randoms (mkStdGen 42)
+#endif
+
+randomIndicesWord :: [Word]
+randomIndicesWord = randoms (mkStdGen 42)
+
+randomIndicesInt :: [Int]
+randomIndicesInt = randoms (mkStdGen 42)
+
+benchReadChimera :: Int -> Benchmark
+benchReadChimera k
+  = bench (show n)
+  $ nf (sum . map (index randomChimera))
+  $ map (.&. (n - 1))
+  $ take (fromIntegral n) randomIndicesWord
+  where
+    n = 1 `shiftL` k
+
+benchReadList :: Int -> Benchmark
+benchReadList k
+  = bench (show n)
+  $ nf (sum . map (randomList !!))
+  $ map (.&. (n - 1))
+  $ take n randomIndicesInt
+  where
+    n = 1 `shiftL` k
+
+#ifdef MIN_VERSION_ral
+benchReadRAL :: Int -> Benchmark
+benchReadRAL k
+  = bench (show n)
+  $ nf (sum . map (randomRAL RAL.!))
+  $ map (.&. (n - 1))
+  $ take n randomIndicesInt
+  where
+    n = 1 `shiftL` k
+#endif
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,5 +1,19 @@
+# 0.4.1.0
+
+* Fix divergence of `fromInfinite` and `fromListWithDef` on infinite inputs.
+
+# 0.4.0.0
+
+* Remove instances `Foldable` and `Traversable`, they are too dangerous to diverge.
+* Add `HalfWord` and `ThirdWord` types,
+  change types of `toZCurve`, `fromZCurve`, `toZCurve3`, `fromZCurve3` accordingly.
+* Add `throughZCurveFix` and `throughZCurveFix3`.
+* Add `imapSubvectors`.
+* Add `prependVector`.
+
 # 0.3.4.0
 
+* Breaking change: remove deprecated `zipSubvectors`, use `zipWithSubvectors`.
 * Add `foldr` catamorphism and `fromInfinite` / `toInfinite` conversions.
 * Add `iterateWithIndex` and `iterateWithIndexM`.
 
diff --git a/chimera.cabal b/chimera.cabal
--- a/chimera.cabal
+++ b/chimera.cabal
@@ -1,15 +1,14 @@
-cabal-version:      2.0
+cabal-version:      2.2
 name:               chimera
-version:            0.3.4.0
-license:            BSD3
+version:            0.4.1.0
+license:            BSD-3-Clause
 license-file:       LICENSE
 copyright:          2017-2019 Bodigrim
 maintainer:         andrew.lelechenko@gmail.com
 author:             Bodigrim
 tested-with:
-    ghc ==9.8.1 ghc ==9.6.3 ghc ==9.4.7 ghc ==9.2.8 ghc ==9.0.2
+    ghc ==9.8.1 ghc ==9.6.3 ghc ==9.4.8 ghc ==9.2.8 ghc ==9.0.2
     ghc ==8.10.7 ghc ==8.8.4 ghc ==8.6.5 ghc ==8.4.4 ghc ==8.2.2
-    ghc ==8.0.2
 
 homepage:           https://github.com/Bodigrim/chimera#readme
 synopsis:
@@ -63,17 +62,19 @@
 
     hs-source-dirs:   src
     other-modules:
-        Data.Chimera.FromIntegral
         Data.Chimera.Compat
+        Data.Chimera.FromIntegral
+        Data.Chimera.Internal
+        Data.Chimera.Memoize
 
     default-language: Haskell2010
     ghc-options:      -Wall -Wcompat
     build-depends:
-        base >=4.9 && <5,
+        base >=4.10 && <5,
         infinite-list <0.2,
         primitive <0.10,
         transformers <0.7,
-        vector <0.14
+        vector <0.14,
 
     if arch(aarch64)
         c-sources:    cbits/aarch64.c
@@ -83,7 +84,7 @@
         build-depends:
             adjunctions <4.5,
             distributive <0.7,
-            mtl <2.4
+            mtl <2.4,
 
 test-suite chimera-test
     type:             exitcode-stdio-1.0
@@ -94,17 +95,22 @@
     build-depends:
         base >=4.5 && <5,
         chimera,
+        infinite-list,
         QuickCheck >=2.10 && <2.15,
         tasty <1.6,
         tasty-hunit <0.11,
         tasty-quickcheck <0.11,
         tasty-smallcheck <0.9,
-        vector
+        vector,
 
 benchmark chimera-bench
     type:             exitcode-stdio-1.0
     main-is:          Bench.hs
     hs-source-dirs:   bench
+    other-modules:
+        Memoize
+        Read
+
     default-language: Haskell2010
     ghc-options:      -Wall -Wcompat
     build-depends:
@@ -113,4 +119,5 @@
         mtl,
         random <1.3,
         tasty >=1.4.2,
-        tasty-bench >=0.3.2 && <0.4
+        tasty-bench >=0.3.2 && <0.4,
+        vector,
diff --git a/src/Data/Chimera.hs b/src/Data/Chimera.hs
--- a/src/Data/Chimera.hs
+++ b/src/Data/Chimera.hs
@@ -1,15 +1,3 @@
-{-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE CPP #-}
-{-# LANGUAGE DeriveTraversable #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE LambdaCase #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TupleSections #-}
-{-# LANGUAGE TypeApplications #-}
-{-# LANGUAGE TypeFamilies #-}
-
 -- |
 -- Module:      Data.Chimera
 -- Copyright:   (c) 2018-2019 Andrew Lelechenko
@@ -41,6 +29,7 @@
 
   -- * Manipulation
   interleave,
+  prependVector,
 
   -- * Elimination
   index,
@@ -60,40 +49,17 @@
   -- * Subvectors
   -- $subvectors
   mapSubvectors,
+  imapSubvectors,
   traverseSubvectors,
   zipWithSubvectors,
   zipWithMSubvectors,
   sliceSubvectors,
 ) where
 
-import Control.Applicative
-import Control.Monad.Fix
-import Control.Monad.Trans.Class
-import qualified Control.Monad.Trans.State.Lazy as LazyState
-import Control.Monad.Zip
-import Data.Bits
-import qualified Data.Foldable as F
-import Data.Functor.Identity
-import Data.List.Infinite (Infinite (..))
-import qualified Data.List.Infinite as Inf
-import qualified Data.Primitive.Array as A
-import qualified Data.Vector as V
-import qualified Data.Vector.Generic as G
-import qualified Data.Vector.Unboxed as U
-import GHC.Exts (fromListN)
 import Prelude hiding (Applicative (..), and, cycle, div, drop, foldr, fromIntegral, iterate, not, or, (*), (^))
 
-#ifdef MIN_VERSION_mtl
-import Control.Monad.Reader (MonadReader, ask, local)
-#endif
-#ifdef MIN_VERSION_distributive
-import Data.Distributive
-#ifdef MIN_VERSION_adjunctions
-import qualified Data.Functor.Rep as Rep
-#endif
-#endif
-
-import Data.Chimera.FromIntegral
+import Data.Chimera.Internal
+import Data.Chimera.Memoize
 
 -- $monadic
 -- Be careful: the stream is infinite, so
@@ -138,658 +104,3 @@
 -- > ch1 = tabulate (Bit . f1)
 -- > ch2 = tabulate (Bit . f2)
 -- > ch3 = zipWithSubvectors (zipBits (.&.)) ch1 ch2
-
--- | Lazy infinite streams with elements from @a@,
--- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
--- Use 'tabulate', 'tabulateFix', etc. to create a stream
--- and 'index' to access its arbitrary elements
--- in constant time.
---
--- @since 0.2.0.0
-newtype Chimera v a = Chimera {unChimera :: A.Array (v a)}
-  deriving
-    ( Functor
-      -- ^ @since 0.2.0.0
-    , Foldable
-      -- ^ @since 0.2.0.0
-    , Traversable
-      -- ^ @since 0.2.0.0
-    )
-
--- | Streams backed by boxed vectors.
---
--- @since 0.3.0.0
-type VChimera = Chimera V.Vector
-
--- | Streams backed by unboxed vectors.
---
--- @since 0.3.0.0
-type UChimera = Chimera U.Vector
-
--- | 'pure' creates a constant stream.
---
--- @since 0.2.0.0
-instance Applicative (Chimera V.Vector) where
-  pure a =
-    Chimera $
-      A.arrayFromListN (bits + 1) $
-        G.singleton a : map (\k -> G.replicate (1 `shiftL` k) a) [0 .. bits - 1]
-  (<*>) = zipWithSubvectors (<*>)
-#if __GLASGOW_HASKELL__ > 801
-  liftA2 f = zipWithSubvectors (liftA2 f)
-#endif
-
--- | @since 0.3.1.0
-instance Monad (Chimera V.Vector) where
-  m >>= f = tabulate $ \w -> index (f (index m w)) w
-
--- | @since 0.3.1.0
-instance MonadFix (Chimera V.Vector) where
-  mfix = tabulate . mfix . fmap index
-
--- | @since 0.3.1.0
-instance MonadZip (Chimera V.Vector) where
-  mzip = zipWithSubvectors mzip
-  mzipWith = zipWithSubvectors . mzipWith
-
-#ifdef MIN_VERSION_mtl
--- | @since 0.3.1.0
-instance MonadReader Word (Chimera V.Vector) where
-  ask = tabulate id
-  local = flip $ (tabulate .) . (.) . index
-#endif
-
-#ifdef MIN_VERSION_distributive
--- | @since 0.3.1.0
-instance Distributive (Chimera V.Vector) where
-  distribute = tabulate . flip (fmap . flip index)
-  collect f = tabulate . flip ((<$>) . (. f) . flip index)
-
-#ifdef MIN_VERSION_adjunctions
--- | @since 0.3.1.0
-instance Rep.Representable (Chimera V.Vector) where
-  type Rep (Chimera V.Vector) = Word
-  tabulate = tabulate
-  index = index
-#endif
-#endif
-
-bits :: Int
-bits = finiteBitSize (0 :: Word)
-
--- | Create a stream of values of a given function.
--- Once created it can be accessed via 'index' or 'toList'.
---
--- >>> ch = tabulate (^ 2) :: UChimera Word
--- >>> index ch 9
--- 81
--- >>> take 10 (toList ch)
--- [0,1,4,9,16,25,36,49,64,81]
---
--- @since 0.2.0.0
-tabulate :: G.Vector v a => (Word -> a) -> Chimera v a
-tabulate f = runIdentity $ tabulateM (pure . f)
-
--- | Similar to 'V.generateM', but for raw arrays.
-generateArrayM :: Monad m => Int -> (Int -> m a) -> m (A.Array a)
-generateArrayM n f = A.arrayFromListN n <$> traverse f [0 .. n - 1]
-
--- | Monadic version of 'tabulate'.
---
--- @since 0.2.0.0
-tabulateM
-  :: (Monad m, G.Vector v a)
-  => (Word -> m a)
-  -> m (Chimera v a)
-tabulateM f = Chimera <$> generateArrayM (bits + 1) tabulateSubVector
-  where
-    tabulateSubVector 0 = G.singleton <$> f 0
-    tabulateSubVector i = G.generateM ii (\j -> f (int2word (ii + j)))
-      where
-        ii = 1 `unsafeShiftL` (i - 1)
-{-# SPECIALIZE tabulateM :: G.Vector v a => (Word -> Identity a) -> Identity (Chimera v a) #-}
-
--- | For a given @f@ create a stream of values of a recursive function 'fix' @f@.
--- Once created it can be accessed via 'index' or 'toList'.
---
--- For example, imagine that we want to tabulate
--- <https://en.wikipedia.org/wiki/Catalan_number Catalan numbers>:
---
--- >>> catalan n = if n == 0 then 1 else sum [ catalan i * catalan (n - 1 - i) | i <- [0 .. n - 1] ]
---
--- Can we find @catalanF@ such that @catalan@ = 'fix' @catalanF@?
--- Just replace all recursive calls to @catalan@ with @f@:
---
--- >>> catalanF f n = if n == 0 then 1 else sum [ f i * f (n - 1 - i) | i <- [0 .. n - 1] ]
---
--- Now we are ready to use 'tabulateFix':
---
--- >>> ch = tabulateFix catalanF :: VChimera Integer
--- >>> index ch 9
--- 4862
--- >>> take 10 (toList ch)
--- [1,1,2,5,14,42,132,429,1430,4862]
---
--- __Note__: Only recursive function calls with decreasing arguments are memoized.
--- If full memoization is desired, use 'tabulateFix'' instead.
---
--- @since 0.2.0.0
-tabulateFix :: G.Vector v a => ((Word -> a) -> Word -> a) -> Chimera v a
-tabulateFix uf = runIdentity $ tabulateFixM ((pure .) . uf . (runIdentity .))
-
--- | Fully memoizing version of 'tabulateFix'.
--- This function will tabulate every recursive call,
--- but might allocate a lot of memory in doing so.
--- For example, the following piece of code calculates the
--- highest number reached by the
--- <https://en.wikipedia.org/wiki/Collatz_conjecture#Statement_of_the_problem Collatz sequence>
--- of a given number, but also allocates tens of gigabytes of memory,
--- because the Collatz sequence will spike to very high numbers.
---
--- >>> collatzF :: (Word -> Word) -> (Word -> Word)
--- >>> collatzF _ 0 = 0
--- >>> collatzF f n = if n <= 2 then 4 else n `max` f (if even n then n `quot` 2 else 3 * n + 1)
--- >>>
--- >>> maximumBy (comparing $ index $ tabulateFix' collatzF) [0..1000000]
--- ...
---
--- Using 'memoizeFix' instead fixes the problem:
---
--- >>> maximumBy (comparing $ memoizeFix collatzF) [0..1000000]
--- 56991483520
---
--- @since 0.3.2.0
-tabulateFix' :: G.Vector v a => ((Word -> a) -> Word -> a) -> Chimera v a
-tabulateFix' uf = runIdentity $ tabulateFixM' ((pure .) . uf . (runIdentity .))
-
--- | Monadic version of 'tabulateFix'.
--- There are no particular guarantees about the order of recursive calls:
--- they may be executed more than once or executed in different order.
--- That said, monadic effects must be idempotent and commutative.
---
--- @since 0.2.0.0
-tabulateFixM
-  :: (Monad m, G.Vector v a)
-  => ((Word -> m a) -> Word -> m a)
-  -> m (Chimera v a)
-tabulateFixM = tabulateFixM_ Downwards
-{-# SPECIALIZE tabulateFixM :: G.Vector v a => ((Word -> Identity a) -> Word -> Identity a) -> Identity (Chimera v a) #-}
-
--- | Monadic version of 'tabulateFix''.
---
--- @since 0.3.3.0
-tabulateFixM'
-  :: forall m v a
-   . (Monad m, G.Vector v a)
-  => ((Word -> m a) -> Word -> m a)
-  -> m (Chimera v a)
-tabulateFixM' = tabulateFixM_ Full
-{-# SPECIALIZE tabulateFixM' :: G.Vector v a => ((Word -> Identity a) -> Word -> Identity a) -> Identity (Chimera v a) #-}
-
--- | Memoization strategy, only used by @tabulateFixM_@.
-data Strategy = Full | Downwards
-
--- | Internal implementation for 'tabulateFixM' and 'tabulateFixM''.
-tabulateFixM_
-  :: forall m v a
-   . (Monad m, G.Vector v a)
-  => Strategy
-  -> ((Word -> m a) -> Word -> m a)
-  -> m (Chimera v a)
-tabulateFixM_ strat f = result
-  where
-    result :: m (Chimera v a)
-    result = Chimera <$> generateArrayM (bits + 1) tabulateSubVector
-
-    tabulateSubVector :: Int -> m (v a)
-    tabulateSubVector 0 =
-      G.singleton <$> case strat of
-        Downwards -> fix f 0
-        Full -> f (\k -> flip index k <$> result) 0
-    tabulateSubVector i = subResult
-      where
-        subResult = G.generateM ii (\j -> f fixF (int2word (ii + j)))
-        subResultBoxed = V.generateM ii (\j -> f fixF (int2word (ii + j)))
-        ii = 1 `unsafeShiftL` (i - 1)
-
-        fixF :: Word -> m a
-        fixF k
-          | k < int2word ii =
-              flip index k <$> result
-          | k <= int2word ii `shiftL` 1 - 1 =
-              (`V.unsafeIndex` (word2int k - ii)) <$> subResultBoxed
-          | otherwise =
-              case strat of
-                Downwards -> f fixF k
-                Full -> flip index k <$> result
-
--- | 'iterate' @f@ @x@ returns an infinite stream, generated by
--- repeated applications of @f@ to @x@.
---
--- It holds that 'index' ('iterate' @f@ @x@) 0 is equal to @x@.
---
--- >>> ch = iterate (+ 1) 0 :: UChimera Int
--- >>> take 10 (toList ch)
--- [0,1,2,3,4,5,6,7,8,9]
---
--- @since 0.3.0.0
-iterate :: G.Vector v a => (a -> a) -> a -> Chimera v a
-iterate f = runIdentity . iterateM (pure . f)
-
--- | Similar to 'G.iterateNM'.
-iterateListNM :: forall a m. Monad m => Int -> (a -> m a) -> a -> m [a]
-iterateListNM n f = if n <= 0 then const (pure []) else go (n - 1)
-  where
-    go :: Int -> a -> m [a]
-    go 0 s = pure [s]
-    go k s = do
-      fs <- f s
-      (s :) <$> go (k - 1) fs
-
--- | Monadic version of 'iterate'.
---
--- @since 0.3.0.0
-iterateM :: (Monad m, G.Vector v a) => (a -> m a) -> a -> m (Chimera v a)
-iterateM f seed = do
-  nextSeed <- f seed
-  let z = G.singleton seed
-  zs <- iterateListNM bits go (G.singleton nextSeed)
-  pure $ Chimera $ fromListN (bits + 1) (z : zs)
-  where
-    go vec = do
-      nextSeed <- f (G.unsafeLast vec)
-      G.iterateNM (G.length vec `shiftL` 1) f nextSeed
-{-# SPECIALIZE iterateM :: G.Vector v a => (a -> Identity a) -> a -> Identity (Chimera v a) #-}
-
--- | 'unfoldr' @f@ @x@ returns an infinite stream, generated by
--- repeated applications of @f@ to @x@, similar to `Data.List.unfoldr`.
---
--- >>> ch = unfoldr (\acc -> (acc * acc, acc + 1)) 0 :: UChimera Int
--- >>> take 10 (toList ch)
--- [0,1,4,9,16,25,36,49,64,81]
---
--- @since 0.3.3.0
-unfoldr :: G.Vector v b => (a -> (b, a)) -> a -> Chimera v b
-unfoldr f = runIdentity . unfoldrM (pure . f)
-
--- | This is not quite satisfactory, see https://github.com/haskell/vector/issues/447
-unfoldrExactVecNM :: forall m a b v. (Monad m, G.Vector v b) => Int -> (a -> m (b, a)) -> a -> m (v b, a)
-unfoldrExactVecNM n f s = flip LazyState.evalStateT s $ do
-  vec <- G.replicateM n f'
-  seed <- LazyState.get
-  pure (vec, seed)
-  where
-    f' :: LazyState.StateT a m b
-    f' = do
-      seed <- LazyState.get
-      (value, newSeed) <- lift (f seed)
-      LazyState.put newSeed
-      pure value
-
--- | Monadic version of 'unfoldr'.
---
--- @since 0.3.3.0
-unfoldrM :: (Monad m, G.Vector v b) => (a -> m (b, a)) -> a -> m (Chimera v b)
-unfoldrM f seed = do
-  let go n s =
-        if n >= bits
-          then pure []
-          else do
-            (vec, s') <- unfoldrExactVecNM (1 `shiftL` n) f s
-            rest <- go (n + 1) s'
-            pure $ vec : rest
-  (z, seed') <- unfoldrExactVecNM 1 f seed
-  zs <- go 0 seed'
-  pure $ Chimera $ fromListN (bits + 1) (z : zs)
-{-# SPECIALIZE unfoldrM :: G.Vector v b => (a -> Identity (b, a)) -> a -> Identity (Chimera v b) #-}
-
--- | 'iterateWithIndex' @f@ @x@ returns an infinite stream, generated by
--- applications of @f@ to a current index and previous value, starting from @x@.
---
--- It holds that 'index' ('iterateWithIndex' @f@ @x@) 0 is equal to @x@.
---
--- >>> ch = iterateWithIndex (+) 100 :: UChimera Word
--- >>> take 10 (toList ch)
--- [100,101,103,106,110,115,121,128,136,145]
---
--- @since 0.3.4.0
-iterateWithIndex :: G.Vector v a => (Word -> a -> a) -> a -> Chimera v a
-iterateWithIndex f = runIdentity . iterateWithIndexM ((pure .) . f)
-
-iterateWithIndexExactVecNM :: forall m a v. (Monad m, G.Vector v a) => Int -> (Word -> a -> m a) -> a -> m (v a)
-iterateWithIndexExactVecNM n f s = G.unfoldrExactNM n go (int2word n, s)
-  where
-    go :: (Word, a) -> m (a, (Word, a))
-    go (i, x) = do
-      x' <- f i x
-      pure (x', (i + 1, x'))
-
--- | Monadic version of 'iterateWithIndex'.
---
--- @since 0.3.4.0
-iterateWithIndexM :: (Monad m, G.Vector v a) => (Word -> a -> m a) -> a -> m (Chimera v a)
-iterateWithIndexM f seed = do
-  nextSeed <- f 1 seed
-  let z = G.singleton seed
-  zs <- iterateListNM bits go (G.singleton nextSeed)
-  pure $ Chimera $ fromListN (bits + 1) (z : zs)
-  where
-    go vec =
-      iterateWithIndexExactVecNM (G.length vec `shiftL` 1) f (G.unsafeLast vec)
-{-# SPECIALIZE iterateWithIndexM :: G.Vector v a => (Word -> a -> Identity a) -> a -> Identity (Chimera v a) #-}
-
-interleaveVec :: G.Vector v a => v a -> v a -> v a
-interleaveVec as bs =
-  G.generate
-    (G.length as `shiftL` 1)
-    (\n -> (if even n then as else bs) G.! (n `shiftR` 1))
-
--- | Intertleave two streams, sourcing even elements from the first one
--- and odd elements from the second one.
---
--- >>> ch = interleave (tabulate id) (tabulate (+ 100)) :: UChimera Word
--- >>> take 10 (toList ch)
--- [0,100,1,101,2,102,3,103,4,104]
---
--- @since 0.3.3.0
-interleave :: G.Vector v a => Chimera v a -> Chimera v a -> Chimera v a
-interleave (Chimera as) (Chimera bs) = Chimera $ A.arrayFromListN (bits + 1) vecs
-  where
-    vecs =
-      A.indexArray as 0
-        : A.indexArray bs 0
-        : map (\i -> interleaveVec (A.indexArray as i) (A.indexArray bs i)) [1 .. bits - 1]
-
--- | Index a stream in a constant time.
---
--- >>> ch = tabulate (^ 2) :: UChimera Word
--- >>> index ch 9
--- 81
---
--- @since 0.2.0.0
-index :: G.Vector v a => Chimera v a -> Word -> a
-index (Chimera vs) i =
-  (vs `A.indexArray` (bits - lz))
-    `G.unsafeIndex` word2int (i .&. complement ((1 `shiftL` (bits - 1)) `unsafeShiftR` lz))
-  where
-    lz :: Int
-    !lz = countLeadingZeros i
-{-# INLINE index #-}
-
--- | Convert a stream to an infinite list.
---
--- >>> ch = tabulate (^ 2) :: UChimera Word
--- >>> take 10 (toList ch)
--- [0,1,4,9,16,25,36,49,64,81]
---
--- @since 0.3.0.0
-toList :: G.Vector v a => Chimera v a -> [a]
-toList (Chimera vs) = foldMap G.toList vs
-
--- | Convert a stream to a proper infinite list.
---
--- @since 0.3.4.0
-toInfinite :: G.Vector v a => Chimera v a -> Infinite a
-toInfinite = foldr (:<)
-
--- | Right-associative fold, necessarily lazy in the accumulator.
--- Any unconditional attempt to force the accumulator even to WHNF
--- will hang the computation. E. g., the following definition isn't productive:
---
--- > import Data.List.NonEmpty (NonEmpty(..))
--- > toNonEmpty = foldr (\a (x :| xs) -> a :| x : xs) :: VChimera a -> NonEmpty a
---
--- One should use lazy patterns, e. g.,
---
--- > toNonEmpty = foldr (\a ~(x :| xs) -> a :| x : xs)
-foldr :: G.Vector v a => (a -> b -> b) -> Chimera v a -> b
-foldr f (Chimera vs) = F.foldr (flip $ G.foldr f) undefined vs
-
-measureOff :: Int -> [a] -> Either Int ([a], [a])
-measureOff n
-  | n <= 0 = Right . ([],)
-  | otherwise = go n
-  where
-    go m [] = Left m
-    go 1 (x : xs) = Right ([x], xs)
-    go m (x : xs) = case go (m - 1) xs of
-      l@Left {} -> l
-      Right (xs', xs'') -> Right (x : xs', xs'')
-
-measureOffVector :: G.Vector v a => Int -> v a -> Either Int (v a, v a)
-measureOffVector n xs
-  | n <= l = Right (G.splitAt n xs)
-  | otherwise = Left (n - l)
-  where
-    l = G.length xs
-
--- | Create a stream of values from a given prefix, followed by default value
--- afterwards.
---
--- @since 0.3.3.0
-fromListWithDef
-  :: G.Vector v a
-  => a
-  -- ^ Default value
-  -> [a]
-  -- ^ Prefix
-  -> Chimera v a
-fromListWithDef a = Chimera . fromListN (bits + 1) . go0
-  where
-    go0 = \case
-      [] -> G.singleton a : map (\k -> G.replicate (1 `shiftL` k) a) [0 .. bits - 1]
-      x : xs -> G.singleton x : go 0 xs
-
-    go k xs = case measureOff kk xs of
-      Left l ->
-        G.fromListN kk (xs ++ replicate l a)
-          : map (\n -> G.replicate (1 `shiftL` n) a) [k + 1 .. bits - 1]
-      Right (ys, zs) -> G.fromListN kk ys : go (k + 1) zs
-      where
-        kk = 1 `shiftL` k
-
--- | Create a stream of values from a given infinite list.
---
--- @since 0.3.4.0
-fromInfinite
-  :: G.Vector v a
-  => Infinite a
-  -> Chimera v a
-fromInfinite = Chimera . fromListN (bits + 1) . go0
-  where
-    go0 (x :< xs) = G.singleton x : go 0 xs
-
-    go k xs = G.fromListN kk ys : go (k + 1) zs
-      where
-        kk = 1 `shiftL` k
-        (ys, zs) = Inf.splitAt kk xs
-
--- | Create a stream of values from a given prefix, followed by default value
--- afterwards.
---
--- @since 0.3.3.0
-fromVectorWithDef
-  :: G.Vector v a
-  => a
-  -- ^ Default value
-  -> v a
-  -- ^ Prefix
-  -> Chimera v a
-fromVectorWithDef a = Chimera . fromListN (bits + 1) . go0
-  where
-    go0 xs = case G.uncons xs of
-      Nothing -> G.singleton a : map (\k -> G.replicate (1 `shiftL` k) a) [0 .. bits - 1]
-      Just (y, ys) -> G.singleton y : go 0 ys
-
-    go k xs = case measureOffVector kk xs of
-      Left l ->
-        (xs G.++ G.replicate l a)
-          : map (\n -> G.replicate (1 `shiftL` n) a) [k + 1 .. bits - 1]
-      Right (ys, zs) -> ys : go (k + 1) zs
-      where
-        kk = 1 `shiftL` k
-
--- | Return an infinite repetition of a given vector.
--- Throw an error on an empty vector.
---
--- >>> ch = cycle (Data.Vector.fromList [4, 2]) :: VChimera Int
--- >>> take 10 (toList ch)
--- [4,2,4,2,4,2,4,2,4,2]
---
--- @since 0.3.0.0
-cycle :: G.Vector v a => v a -> Chimera v a
-cycle vec = case l of
-  0 -> error "Data.Chimera.cycle: empty list"
-  _ -> tabulate (G.unsafeIndex vec . word2int . (`rem` l))
-  where
-    l = int2word $ G.length vec
-
--- | Memoize a function:
--- repeating calls to 'memoize' @f@ @n@
--- would compute @f@ @n@ only once
--- and cache the result in 'VChimera'.
--- This is just a shortcut for 'index' '.' 'tabulate'.
--- When @a@ is 'U.Unbox', it is faster to use
--- 'index' ('tabulate' @f@ :: 'UChimera' @a@).
---
--- prop> memoize f n = f n
---
--- @since 0.3.0.0
-memoize :: (Word -> a) -> (Word -> a)
-memoize = index @V.Vector . tabulate
-
--- | For a given @f@ memoize a recursive function 'fix' @f@,
--- caching results in 'VChimera'.
--- This is just a shortcut for 'index' '.' 'tabulateFix'.
--- When @a@ is 'U.Unbox', it is faster to use
--- 'index' ('tabulateFix' @f@ :: 'UChimera' @a@).
---
--- prop> memoizeFix f n = fix f n
---
--- For example, imagine that we want to memoize
--- <https://en.wikipedia.org/wiki/Fibonacci_number Fibonacci numbers>:
---
--- >>> fibo n = if n < 2 then toInteger n else fibo (n - 1) + fibo (n - 2)
---
--- Can we find @fiboF@ such that @fibo@ = 'fix' @fiboF@?
--- Just replace all recursive calls to @fibo@ with @f@:
---
--- >>> fiboF f n = if n < 2 then toInteger n else f (n - 1) + f (n - 2)
---
--- Now we are ready to use 'memoizeFix':
---
--- >>> memoizeFix fiboF 10
--- 55
--- >>> memoizeFix fiboF 100
--- 354224848179261915075
---
--- This function can be used even when arguments
--- of recursive calls are not strictly decreasing,
--- but they might not get memoized. If this is not desired
--- use 'tabulateFix'' instead.
--- For example, here is a routine to measure the length of
--- <https://oeis.org/A006577 Collatz sequence>:
---
--- >>> collatzF f n = if n <= 1 then 0 else 1 + f (if even n then n `quot` 2 else 3 * n + 1)
--- >>> memoizeFix collatzF 27
--- 111
---
--- @since 0.3.0.0
-memoizeFix :: ((Word -> a) -> Word -> a) -> (Word -> a)
-memoizeFix = index @V.Vector . tabulateFix
-
--- | Map subvectors of a stream, using a given length-preserving function.
---
--- @since 0.3.0.0
-mapSubvectors
-  :: (G.Vector u a, G.Vector v b)
-  => (u a -> v b)
-  -> Chimera u a
-  -> Chimera v b
-mapSubvectors f = runIdentity . traverseSubvectors (pure . f)
-
--- | Traverse subvectors of a stream, using a given length-preserving function.
---
--- Be careful, because similar to 'tabulateM', only lazy monadic effects can
--- be executed in a finite time: lazy state monad is fine, but strict one is
--- not.
---
--- @since 0.3.3.0
-traverseSubvectors
-  :: (G.Vector u a, G.Vector v b, Applicative m)
-  => (u a -> m (v b))
-  -> Chimera u a
-  -> m (Chimera v b)
-traverseSubvectors f (Chimera bs) = Chimera <$> traverse safeF bs
-  where
-    -- Computing vector length is cheap, so let's check that @f@ preserves length.
-    safeF x =
-      ( \fx ->
-          if G.length x == G.length fx
-            then fx
-            else error "traverseSubvectors: the function is not length-preserving"
-      )
-        <$> f x
-{-# SPECIALIZE traverseSubvectors :: (G.Vector u a, G.Vector v b) => (u a -> Identity (v b)) -> Chimera u a -> Identity (Chimera v b) #-}
-
--- | Zip subvectors from two streams, using a given length-preserving function.
---
--- @since 0.3.3.0
-zipWithSubvectors
-  :: (G.Vector u a, G.Vector v b, G.Vector w c)
-  => (u a -> v b -> w c)
-  -> Chimera u a
-  -> Chimera v b
-  -> Chimera w c
-zipWithSubvectors f = (runIdentity .) . zipWithMSubvectors ((pure .) . f)
-
--- | Zip subvectors from two streams, using a given monadic length-preserving function.
--- Caveats for 'tabulateM' and 'traverseSubvectors' apply.
---
--- @since 0.3.3.0
-zipWithMSubvectors
-  :: (G.Vector u a, G.Vector v b, G.Vector w c, Applicative m)
-  => (u a -> v b -> m (w c))
-  -> Chimera u a
-  -> Chimera v b
-  -> m (Chimera w c)
-zipWithMSubvectors f (Chimera bs1) (Chimera bs2) = Chimera <$> sequenceA (mzipWith safeF bs1 bs2)
-  where
-    -- Computing vector length is cheap, so let's check that @f@ preserves length.
-    safeF x y =
-      ( \fx ->
-          if G.length x == G.length fx
-            then fx
-            else error "traverseSubvectors: the function is not length-preserving"
-      )
-        <$> f x y
-{-# SPECIALIZE zipWithMSubvectors :: (G.Vector u a, G.Vector v b, G.Vector w c) => (u a -> v b -> Identity (w c)) -> Chimera u a -> Chimera v b -> Identity (Chimera w c) #-}
-
--- | Take a slice of 'Chimera', represented as a list on consecutive subvectors.
---
--- @since 0.3.3.0
-sliceSubvectors
-  :: G.Vector v a
-  => Int
-  -- ^ How many initial elements to drop?
-  -> Int
-  -- ^ How many subsequent elements to take?
-  -> Chimera v a
-  -> [v a]
-sliceSubvectors off len = doTake len . doDrop off . F.toList . unChimera
-  where
-    doTake !_ [] = []
-    doTake n (x : xs)
-      | n <= 0 = []
-      | n >= l = x : doTake (n - l) xs
-      | otherwise = [G.take n x]
-      where
-        l = G.length x
-
-    doDrop !_ [] = []
-    doDrop n (x : xs)
-      | n <= 0 = x : xs
-      | l <= n = doDrop (n - l) xs
-      | otherwise = G.drop n x : xs
-      where
-        l = G.length x
diff --git a/src/Data/Chimera/ContinuousMapping.hs b/src/Data/Chimera/ContinuousMapping.hs
--- a/src/Data/Chimera/ContinuousMapping.hs
+++ b/src/Data/Chimera/ContinuousMapping.hs
@@ -1,3 +1,8 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE TypeApplications #-}
+
 -- |
 -- Module:      Data.Chimera.ContinuousMapping
 -- Copyright:   (c) 2017 Andrew Lelechenko
@@ -63,26 +68,39 @@
 -- Namely, cast the state to 'Word' @->@ 'Bool', ready for memoization:
 --
 -- > cast :: (Int -> Int -> Bool) -> (Word -> Bool)
--- > cast f = \n -> let (x, y) = fromZCurve n in
--- >  f (wordToInt x) (wordToInt y)
+-- > cast f = \n -> let (x, y) = fromZCurve n in f (fromHalf x) (fromHalf y)
+-- >   where
+-- >     fromHalf :: HalfWord -> Int
+-- >     fromHalf = wordToInt . fromIntegral @HalfWord @Word
 --
 -- and then back:
 --
 -- > uncast :: (Word -> Bool) -> (Int -> Int -> Bool)
--- > uncast g = \x y -> g (toZCurve (intToWord x) (intToWord y))
+-- > uncast g = \x y -> g (toZCurve (toHalf x) (toHalf y))
+-- >   where
+-- >     toHalf :: Int -> HalfWord
+-- >     toHalf = fromIntegral @Word @HalfWord . intToWord
 module Data.Chimera.ContinuousMapping (
   intToWord,
   wordToInt,
+  HalfWord,
   toZCurve,
   fromZCurve,
+  throughZCurveFix,
+  ThirdWord,
   toZCurve3,
   fromZCurve3,
+  throughZCurveFix3,
 ) where
 
+import Data.Bifunctor
 import Data.Bits
 import Data.Chimera.FromIntegral
+import Data.Coerce
 import Data.Word
 
+#include "MachDeps.h"
+
 -- | Total map, which satisfies
 --
 -- prop> abs (intToWord x - intToWord y) <= 2 * abs (x - y)
@@ -111,16 +129,27 @@
 wordToInt w = word2int $ (if w .&. 1 == 0 then id else complement) (w `shiftR` 1)
 {-# INLINE wordToInt #-}
 
+-- | 32 bits on 64-bit architecture, 16 bits on 32-bit architecture.
+--
+-- To create a value of type 'HalfWord' use 'fromIntegral'.
+--
+-- @since 0.4.0.0
+#if WORD_SIZE_IN_BITS == 64
+newtype HalfWord = HalfWord Word32
+  deriving newtype (Eq, Ord, Show, Read, Bits, FiniteBits, Bounded, Enum, Num, Integral, Real)
+#else
+newtype HalfWord = HalfWord Word16
+  deriving newtype (Eq, Ord, Show, Read, Bits, FiniteBits, Bounded, Enum, Num, Integral, Real)
+#endif
+
 -- | Total map from plain to line, continuous almost everywhere.
 -- See <https://en.wikipedia.org/wiki/Z-order_curve Z-order curve>.
 --
--- Only lower halfs of bits of arguments are used (32 bits on 64-bit architecture).
---
 -- >>> [ toZCurve x y | x <- [0..3], y <- [0..3] ]
 -- [0,2,8,10,1,3,9,11,4,6,12,14,5,7,13,15]
 --
 -- @since 0.2.0.0
-toZCurve :: Word -> Word -> Word
+toZCurve :: HalfWord -> HalfWord -> Word
 toZCurve x y = part1by1 y `shiftL` 1 .|. part1by1 x
 
 -- | Inverse for 'toZCurve'.
@@ -130,20 +159,102 @@
 -- [(0,0),(1,0),(0,1),(1,1),(2,0),(3,0),(2,1),(3,1),(0,2),(1,2),(0,3),(1,3),(2,2),(3,2),(2,3),(3,3)]
 --
 -- @since 0.2.0.0
-fromZCurve :: Word -> (Word, Word)
+fromZCurve :: Word -> (HalfWord, HalfWord)
 fromZCurve z = (compact1by1 z, compact1by1 (z `shiftR` 1))
 
+-- | Convert a function of two 'HalfWord's to a function of one 'Word'.
+contramapFromZCurve
+  :: (HalfWord -> HalfWord -> a)
+  -> (Word -> a)
+contramapFromZCurve f = uncurry f . fromZCurve
+
+-- | Convert a function of one 'Word' to a function of two 'HalfWord's.
+contramapToZCurve
+  :: (Word -> a)
+  -> (HalfWord -> HalfWord -> a)
+contramapToZCurve f = (f .) . toZCurve
+
+-- | For an input function @f@ return function @g@ such that
+-- 'Data.Function.fix' @f@ = ('Data.Function.fix' @g@ '.') '.' 'toZCurve'.
+--
+-- @since 0.4.0.0
+throughZCurveFix
+  :: ((HalfWord -> HalfWord -> a) -> (HalfWord -> HalfWord -> a))
+  -> (Word -> a)
+  -> (Word -> a)
+throughZCurveFix f = contramapFromZCurve . f . contramapToZCurve
+
+-- | 21 bits on 64-bit architecture, 10 bits on 32-bit architecture.
+--
+-- To create a value of type 'ThirdWord' use 'fromIntegral'.
+--
+-- @since 0.4.0.0
+newtype ThirdWord = ThirdWord Word32
+  deriving newtype (Eq, Ord, Show)
+
+mkThirdWord :: Word32 -> ThirdWord
+mkThirdWord n = t
+  where
+    t = ThirdWord (n .&. ((1 `shiftL` finiteBitSize t) - 1))
+
+instance Read ThirdWord where
+  readsPrec = (fmap (first mkThirdWord) .) . readsPrec
+
+instance Bits ThirdWord where
+  (.&.) = coerce ((.&.) @Word32)
+  (.|.) = coerce ((.|.) @Word32)
+  xor = coerce (xor @Word32)
+  complement (ThirdWord n) = mkThirdWord (complement n)
+  shift (ThirdWord n) k = mkThirdWord (shift n k)
+  bitSize = finiteBitSize
+  bitSizeMaybe = Just . finiteBitSize
+  isSigned = coerce (isSigned @Word32)
+  testBit = coerce (testBit @Word32)
+  bit = mkThirdWord . bit
+  popCount = coerce (popCount @Word32)
+
+  rotate t k'
+    | k == 0 = t
+    | otherwise = (t `shiftL` k) .|. (t `shiftR` (fbs - k))
+    where
+      fbs = finiteBitSize t
+      k = k' `mod` fbs
+
+instance FiniteBits ThirdWord where
+  finiteBitSize = const $ finiteBitSize (0 :: Word) `quot` 3
+
+instance Bounded ThirdWord where
+  minBound = mkThirdWord minBound
+  maxBound = mkThirdWord maxBound
+
+instance Enum ThirdWord where
+  toEnum = mkThirdWord . toEnum
+  fromEnum = coerce (fromEnum @Word32)
+
+instance Num ThirdWord where
+  ThirdWord x + ThirdWord y = mkThirdWord (x + y)
+  ThirdWord x * ThirdWord y = mkThirdWord (x * y)
+  negate (ThirdWord x) = mkThirdWord (negate x)
+  abs = coerce (abs @Word32)
+  signum = coerce (signum @Word32)
+  fromInteger = mkThirdWord . fromInteger
+
+instance Real ThirdWord where
+  toRational = coerce (toRational @Word32)
+
+instance Integral ThirdWord where
+  quotRem = coerce (quotRem @Word32)
+  toInteger = coerce (toInteger @Word32)
+
 -- | Total map from space to line, continuous almost everywhere.
 -- See <https://en.wikipedia.org/wiki/Z-order_curve Z-order curve>.
 --
--- Only lower thirds of bits of arguments are used (21 bits on 64-bit architecture).
---
 -- >>> [ toZCurve3 x y z | x <- [0..3], y <- [0..3], z <- [0..3] ]
 -- [0,4,32,36,2,6,34,38,16,20,48,52,18,22,50,54,1,5,33,37,3,7,35,39,17,21,49,53,19,23,51,55,
 -- 8,12,40,44,10,14,42,46,24,28,56,60,26,30,58,62,9,13,41,45,11,15,43,47,25,29,57,61,27,31,59,63]
 --
 -- @since 0.2.0.0
-toZCurve3 :: Word -> Word -> Word -> Word
+toZCurve3 :: ThirdWord -> ThirdWord -> ThirdWord -> Word
 toZCurve3 x y z = part1by2 z `shiftL` 2 .|. part1by2 y `shiftL` 1 .|. part1by2 x
 
 -- | Inverse for 'toZCurve3'.
@@ -156,11 +267,35 @@
 --  (0,2,2),(1,2,2),(0,3,2),(1,3,2),(0,2,3),(1,2,3),(0,3,3),(1,3,3),(2,2,2),(3,2,2),(2,3,2),(3,3,2),(2,2,3),(3,2,3),(2,3,3),(3,3,3)]
 --
 -- @since 0.2.0.0
-fromZCurve3 :: Word -> (Word, Word, Word)
+fromZCurve3 :: Word -> (ThirdWord, ThirdWord, ThirdWord)
 fromZCurve3 z = (compact1by2 z, compact1by2 (z `shiftR` 1), compact1by2 (z `shiftR` 2))
 
+-- | Convert a function of two 'ThirdWord's to a function of one 'Word'.
+contramapFromZCurve3
+  :: (ThirdWord -> ThirdWord -> ThirdWord -> a)
+  -> (Word -> a)
+contramapFromZCurve3 f = uncurry3 f . fromZCurve3
+  where
+    uncurry3 func (a, b, c) = func a b c
+
+-- | Convert a function of one 'Word' to a function of two 'ThirdWord's.
+contramapToZCurve3
+  :: (Word -> a)
+  -> (ThirdWord -> ThirdWord -> ThirdWord -> a)
+contramapToZCurve3 f = ((f .) .) . toZCurve3
+
+-- | For an input function @f@ return function @g@ such that
+-- 'Data.Function.fix' @f@ = (('Data.Function.fix' @g@ '.') '.') '.' 'toZCurve3'.
+--
+-- @since 0.4.0.0
+throughZCurveFix3
+  :: ((ThirdWord -> ThirdWord -> ThirdWord -> a) -> (ThirdWord -> ThirdWord -> ThirdWord -> a))
+  -> (Word -> a)
+  -> (Word -> a)
+throughZCurveFix3 f = contramapFromZCurve3 . f . contramapToZCurve3
+
 -- Inspired by https://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/
-part1by1 :: Word -> Word
+part1by1 :: HalfWord -> Word
 part1by1 x = fromIntegral (x5 :: Word64)
   where
     x0 = fromIntegral x .&. 0x00000000ffffffff
@@ -171,7 +306,7 @@
     x5 = (x4 `xor` (x4 `shiftL` 1)) .&. 0x5555555555555555
 
 -- Inspired by https://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/
-part1by2 :: Word -> Word
+part1by2 :: ThirdWord -> Word
 part1by2 x = fromIntegral (x5 :: Word64)
   where
     x0 = fromIntegral x .&. 0x00000000ffffffff
@@ -182,7 +317,7 @@
     x5 = (x4 `xor` (x4 `shiftL` 2)) .&. 0x1249249249249249
 
 -- Inspired by https://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/
-compact1by1 :: Word -> Word
+compact1by1 :: Word -> HalfWord
 compact1by1 x = fromIntegral (x5 :: Word64)
   where
     x0 = fromIntegral x .&. 0x5555555555555555
@@ -193,7 +328,7 @@
     x5 = (x4 `xor` (x4 `shiftR` 16)) .&. 0x00000000ffffffff
 
 -- Inspired by https://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/
-compact1by2 :: Word -> Word
+compact1by2 :: Word -> ThirdWord
 compact1by2 x = fromIntegral (x5 :: Word64)
   where
     x0 = fromIntegral x .&. 0x1249249249249249
diff --git a/src/Data/Chimera/Internal.hs b/src/Data/Chimera/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Chimera/Internal.hs
@@ -0,0 +1,785 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE DeriveTraversable #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE LambdaCase #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE PolyKinds #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TupleSections #-}
+{-# LANGUAGE TypeApplications #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE ViewPatterns #-}
+
+-- |
+-- Module:      Data.Chimera.Internal
+-- Copyright:   (c) 2018-2019 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+module Data.Chimera.Internal (
+  -- * Chimera
+  Chimera,
+  VChimera,
+  UChimera,
+
+  -- * Construction
+  tabulate,
+  tabulateFix,
+  tabulateFix',
+  iterate,
+  iterateWithIndex,
+  unfoldr,
+  cycle,
+  fromListWithDef,
+  fromVectorWithDef,
+  fromInfinite,
+
+  -- * Manipulation
+  interleave,
+  prependVector,
+
+  -- * Elimination
+  index,
+  foldr,
+  toList,
+  toInfinite,
+
+  -- * Monadic construction
+  tabulateM,
+  tabulateFixM,
+  tabulateFixM',
+  iterateM,
+  iterateWithIndexM,
+  unfoldrM,
+
+  -- * Subvectors
+  mapSubvectors,
+  imapSubvectors,
+  traverseSubvectors,
+  zipWithSubvectors,
+  zipWithMSubvectors,
+  sliceSubvectors,
+) where
+
+import Control.Applicative
+import Control.Monad.Fix
+import Control.Monad.Trans.Class
+import qualified Control.Monad.Trans.State.Lazy as LazyState
+import Control.Monad.Zip
+import Data.Bits
+import Data.Coerce
+import qualified Data.Foldable as F
+import Data.Functor.Identity
+import Data.List.Infinite (Infinite (..))
+import qualified Data.List.Infinite as Inf
+import qualified Data.Primitive.Array as A
+import Data.Typeable
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+import GHC.Exts (fromListN)
+import Prelude hiding (Applicative (..), and, cycle, div, drop, foldr, fromIntegral, iterate, not, or, (*), (^))
+
+#ifdef MIN_VERSION_mtl
+import Control.Monad.Reader (MonadReader, ask, local)
+#endif
+#ifdef MIN_VERSION_distributive
+import Data.Distributive
+#ifdef MIN_VERSION_adjunctions
+import qualified Data.Functor.Rep as Rep
+#endif
+#endif
+
+import Data.Chimera.FromIntegral
+
+-- | Lazy infinite streams with elements from @a@,
+-- backed by a 'G.Vector' @v@ (boxed, unboxed, storable, etc.).
+-- Use 'tabulate', 'tabulateFix', etc. to create a stream
+-- and 'index' to access its arbitrary elements
+-- in constant time.
+--
+-- @since 0.2.0.0
+newtype Chimera v a = Chimera {unChimera :: A.Array (v a)}
+  deriving
+    ( Functor
+      -- ^ @since 0.2.0.0
+    )
+
+-- | Streams backed by boxed vectors.
+--
+-- @since 0.3.0.0
+type VChimera = Chimera V.Vector
+
+-- | Streams backed by unboxed vectors.
+--
+-- @since 0.3.0.0
+type UChimera = Chimera U.Vector
+
+-- | 'pure' creates a constant stream.
+--
+-- @since 0.2.0.0
+instance Applicative (Chimera V.Vector) where
+  pure a =
+    Chimera $
+      A.arrayFromListN (bits + 1) $
+        G.singleton a : map (\k -> G.replicate (1 `shiftL` k) a) [0 .. bits - 1]
+  (<*>) = zipWithSubvectors (<*>)
+  liftA2 f = zipWithSubvectors (liftA2 f)
+
+-- | @since 0.3.1.0
+instance Monad (Chimera V.Vector) where
+  m >>= f = tabulate $ \w -> index (f (index m w)) w
+
+-- | @since 0.3.1.0
+instance MonadFix (Chimera V.Vector) where
+  mfix = tabulate . mfix . fmap index
+
+-- | @since 0.3.1.0
+instance MonadZip (Chimera V.Vector) where
+  mzip = zipWithSubvectors mzip
+  mzipWith = zipWithSubvectors . mzipWith
+
+#ifdef MIN_VERSION_mtl
+-- | @since 0.3.1.0
+instance MonadReader Word (Chimera V.Vector) where
+  ask = tabulate id
+  local = flip $ (tabulate .) . (.) . index
+#endif
+
+#ifdef MIN_VERSION_distributive
+-- | @since 0.3.1.0
+instance Distributive (Chimera V.Vector) where
+  distribute = tabulate . flip (fmap . flip index)
+  collect f = tabulate . flip ((<$>) . (. f) . flip index)
+
+#ifdef MIN_VERSION_adjunctions
+-- | @since 0.3.1.0
+instance Rep.Representable (Chimera V.Vector) where
+  type Rep (Chimera V.Vector) = Word
+  tabulate = tabulate
+  index = index
+#endif
+#endif
+
+bits :: Int
+bits = finiteBitSize (0 :: Word)
+
+-- | Create a stream of values of a given function.
+-- Once created it can be accessed via 'index' or 'toList'.
+--
+-- >>> ch = tabulate (^ 2) :: UChimera Word
+-- >>> index ch 9
+-- 81
+-- >>> take 10 (toList ch)
+-- [0,1,4,9,16,25,36,49,64,81]
+--
+-- Note that @a@ could be a function type itself,
+-- so one can tabulate a function of multiple arguments
+-- as a nested 'Chimera' of 'Chimera's.
+--
+-- @since 0.2.0.0
+tabulate :: G.Vector v a => (Word -> a) -> Chimera v a
+tabulate f = runIdentity $ tabulateM (coerce f)
+{-# INLINEABLE tabulate #-}
+
+-- | Similar to 'V.generateM', but for raw arrays.
+generateArrayM :: Monad m => Int -> (Int -> m a) -> m (A.Array a)
+generateArrayM n f = A.arrayFromListN n <$> traverse f [0 .. n - 1]
+
+-- | Monadic version of 'tabulate'.
+--
+-- @since 0.2.0.0
+tabulateM
+  :: (Monad m, G.Vector v a)
+  => (Word -> m a)
+  -> m (Chimera v a)
+tabulateM f = Chimera <$> generateArrayM (bits + 1) tabulateSubVector
+  where
+    tabulateSubVector 0 = G.singleton <$> f 0
+    tabulateSubVector i = G.generateM ii (\j -> f (int2word (ii + j)))
+      where
+        ii = 1 `unsafeShiftL` (i - 1)
+{-# INLINEABLE tabulateM #-}
+{-# SPECIALIZE tabulateM :: G.Vector v a => (Word -> Identity a) -> Identity (Chimera v a) #-}
+
+-- | For a given @f@ create a stream of values of a recursive function 'Data.Function.fix' @f@.
+-- Once created it can be accessed via 'index' or 'toList'.
+--
+-- For example, imagine that we want to tabulate
+-- <https://en.wikipedia.org/wiki/Catalan_number Catalan numbers>:
+--
+-- >>> catalan n = if n == 0 then 1 else sum [ catalan i * catalan (n - 1 - i) | i <- [0 .. n - 1] ]
+--
+-- Can we find @catalanF@ such that @catalan@ = 'Data.Function.fix' @catalanF@?
+-- Just replace all recursive calls to @catalan@ with @f@:
+--
+-- >>> catalanF f n = if n == 0 then 1 else sum [ f i * f (n - 1 - i) | i <- [0 .. n - 1] ]
+--
+-- Now we are ready to use 'tabulateFix':
+--
+-- >>> ch = tabulateFix catalanF :: VChimera Integer
+-- >>> index ch 9
+-- 4862
+-- >>> take 10 (toList ch)
+-- [1,1,2,5,14,42,132,429,1430,4862]
+--
+-- __Note__: Only recursive function calls with decreasing arguments are memoized.
+-- If full memoization is desired, use 'tabulateFix'' instead.
+--
+-- Using unboxed \/ storable \/ primitive vectors with 'tabulateFix' is not always a win:
+-- the internal memoizing routine necessarily uses boxed vectors to achieve
+-- a certain degree of laziness, so converting to 'UChimera' is extra work.
+-- This could pay off in a long run by reducing memory residence though.
+--
+-- @since 0.2.0.0
+tabulateFix :: (G.Vector v a, Typeable v) => ((Word -> a) -> Word -> a) -> Chimera v a
+tabulateFix uf = runIdentity $ tabulateFixM (coerce uf)
+{-# INLINEABLE tabulateFix #-}
+
+-- | Fully memoizing version of 'tabulateFix'.
+-- This function will tabulate every recursive call,
+-- but might allocate a lot of memory in doing so.
+-- For example, the following piece of code calculates the
+-- highest number reached by the
+-- <https://en.wikipedia.org/wiki/Collatz_conjecture#Statement_of_the_problem Collatz sequence>
+-- of a given number, but also allocates tens of gigabytes of memory,
+-- because the Collatz sequence will spike to very high numbers.
+--
+-- >>> collatzF :: (Word -> Word) -> (Word -> Word)
+-- >>> collatzF _ 0 = 0
+-- >>> collatzF f n = if n <= 2 then 4 else n `max` f (if even n then n `quot` 2 else 3 * n + 1)
+-- >>>
+-- >>> maximumBy (comparing $ index $ tabulateFix' collatzF) [0..1000000]
+-- ...
+--
+-- Using 'Data.Chimera.memoizeFix' instead fixes the problem:
+--
+-- >>> maximumBy (comparing $ memoizeFix collatzF) [0..1000000]
+-- 56991483520
+--
+-- Since 'tabulateFix'' memoizes all recursive calls, even with increasing argument,
+-- you most likely do not want to use it with anything else than boxed vectors ('VChimera').
+--
+-- @since 0.3.2.0
+tabulateFix' :: (G.Vector v a, Typeable v) => ((Word -> a) -> Word -> a) -> Chimera v a
+tabulateFix' uf = runIdentity $ tabulateFixM' (coerce uf)
+{-# INLINEABLE tabulateFix' #-}
+
+-- | Monadic version of 'tabulateFix'.
+-- There are no particular guarantees about the order of recursive calls:
+-- they may be executed more than once or executed in different order.
+-- That said, monadic effects must be idempotent and commutative.
+--
+-- @since 0.2.0.0
+tabulateFixM
+  :: (Monad m, G.Vector v a, Typeable v)
+  => ((Word -> m a) -> Word -> m a)
+  -> m (Chimera v a)
+tabulateFixM = tabulateFixM_ Downwards
+{-# INLINEABLE tabulateFixM #-}
+{-# SPECIALIZE tabulateFixM :: (G.Vector v a, Typeable v) => ((Word -> Identity a) -> Word -> Identity a) -> Identity (Chimera v a) #-}
+
+-- | Monadic version of 'tabulateFix''.
+--
+-- @since 0.3.3.0
+tabulateFixM'
+  :: forall m v a
+   . (Monad m, G.Vector v a, Typeable v)
+  => ((Word -> m a) -> Word -> m a)
+  -> m (Chimera v a)
+tabulateFixM' = tabulateFixM_ Full
+{-# INLINEABLE tabulateFixM' #-}
+{-# SPECIALIZE tabulateFixM' :: (G.Vector v a, Typeable v) => ((Word -> Identity a) -> Word -> Identity a) -> Identity (Chimera v a) #-}
+
+-- | Memoization strategy, only used by @tabulateFixM_@.
+data Strategy = Full | Downwards
+
+-- | Internal implementation for 'tabulateFixM' and 'tabulateFixM''.
+tabulateFixM_
+  :: forall m v a
+   . (Monad m, G.Vector v a, Typeable v)
+  => Strategy
+  -> ((Word -> m a) -> Word -> m a)
+  -> m (Chimera v a)
+tabulateFixM_ strat f = result
+  where
+    result :: m (Chimera v a)
+    result = Chimera <$> generateArrayM (bits + 1) tabulateSubVector
+
+    tabulateSubVector :: Int -> m (v a)
+    tabulateSubVector 0 =
+      G.singleton <$> case strat of
+        Downwards -> fix f 0
+        Full -> f (\k -> flip index k <$> result) 0
+    tabulateSubVector i = subResult
+      where
+        subResult = fromBoxedVector <$> subResultBoxed
+        subResultBoxed = V.generateM ii (\j -> f fixF (int2word (ii + j)))
+        ii = 1 `unsafeShiftL` (i - 1)
+
+        fixF :: Word -> m a
+        fixF k
+          | k < int2word ii =
+              flip index k <$> result
+          | k <= int2word ii `shiftL` 1 - 1 =
+              (`V.unsafeIndex` (word2int k - ii)) <$> subResultBoxed
+          | otherwise =
+              case strat of
+                Downwards -> f fixF k
+                Full -> flip index k <$> result
+-- It's crucial to inline into tabulateFixM and tabulateFixM'.
+{-# INLINE tabulateFixM_ #-}
+
+fromBoxedVector :: forall v a. (G.Vector v a, Typeable v) => V.Vector a -> v a
+fromBoxedVector = case eqT @V.Vector @v of
+  Just Refl -> id
+  Nothing -> G.convert
+
+-- | 'iterate' @f@ @x@ returns an infinite stream, generated by
+-- repeated applications of @f@ to @x@.
+--
+-- It holds that 'index' ('iterate' @f@ @x@) 0 is equal to @x@.
+--
+-- >>> ch = iterate (+ 1) 0 :: UChimera Int
+-- >>> take 10 (toList ch)
+-- [0,1,2,3,4,5,6,7,8,9]
+--
+-- @since 0.3.0.0
+iterate :: G.Vector v a => (a -> a) -> a -> Chimera v a
+iterate f = runIdentity . iterateM (coerce f)
+
+-- | Similar to 'G.iterateNM'.
+iterateListNM :: forall a m. Monad m => Int -> (a -> m a) -> a -> m [a]
+iterateListNM n f = if n <= 0 then const (pure []) else go (n - 1)
+  where
+    go :: Int -> a -> m [a]
+    go 0 s = pure [s]
+    go k s = do
+      fs <- f s
+      (s :) <$> go (k - 1) fs
+
+-- | Monadic version of 'iterate'.
+--
+-- @since 0.3.0.0
+iterateM :: (Monad m, G.Vector v a) => (a -> m a) -> a -> m (Chimera v a)
+iterateM f seed = do
+  nextSeed <- f seed
+  let z = G.singleton seed
+  zs <- iterateListNM bits go (G.singleton nextSeed)
+  pure $ Chimera $ fromListN (bits + 1) (z : zs)
+  where
+    go vec = do
+      nextSeed <- f (G.unsafeLast vec)
+      G.iterateNM (G.length vec `shiftL` 1) f nextSeed
+{-# SPECIALIZE iterateM :: G.Vector v a => (a -> Identity a) -> a -> Identity (Chimera v a) #-}
+
+-- | 'unfoldr' @f@ @x@ returns an infinite stream, generated by
+-- repeated applications of @f@ to @x@, similar to `Data.List.unfoldr`.
+--
+-- >>> ch = unfoldr (\acc -> (acc * acc, acc + 1)) 0 :: UChimera Int
+-- >>> take 10 (toList ch)
+-- [0,1,4,9,16,25,36,49,64,81]
+--
+-- @since 0.3.3.0
+unfoldr :: G.Vector v b => (a -> (b, a)) -> a -> Chimera v b
+unfoldr f = runIdentity . unfoldrM (coerce f)
+
+-- | This is not quite satisfactory, see https://github.com/haskell/vector/issues/447
+unfoldrExactVecNM :: forall m a b v. (Monad m, G.Vector v b) => Int -> (a -> m (b, a)) -> a -> m (v b, a)
+unfoldrExactVecNM n f s = flip LazyState.evalStateT s $ do
+  vec <- G.replicateM n f'
+  seed <- LazyState.get
+  pure (vec, seed)
+  where
+    f' :: LazyState.StateT a m b
+    f' = do
+      seed <- LazyState.get
+      (value, newSeed) <- lift (f seed)
+      LazyState.put newSeed
+      pure value
+
+-- | Monadic version of 'unfoldr'.
+--
+-- @since 0.3.3.0
+unfoldrM :: (Monad m, G.Vector v b) => (a -> m (b, a)) -> a -> m (Chimera v b)
+unfoldrM f seed = do
+  let go n s =
+        if n >= bits
+          then pure []
+          else do
+            (vec, s') <- unfoldrExactVecNM (1 `shiftL` n) f s
+            rest <- go (n + 1) s'
+            pure $ vec : rest
+  (z, seed') <- unfoldrExactVecNM 1 f seed
+  zs <- go 0 seed'
+  pure $ Chimera $ fromListN (bits + 1) (z : zs)
+{-# SPECIALIZE unfoldrM :: G.Vector v b => (a -> Identity (b, a)) -> a -> Identity (Chimera v b) #-}
+
+-- | 'iterateWithIndex' @f@ @x@ returns an infinite stream, generated by
+-- applications of @f@ to a current index and previous value, starting from @x@.
+--
+-- It holds that 'index' ('iterateWithIndex' @f@ @x@) 0 is equal to @x@.
+--
+-- >>> ch = iterateWithIndex (+) 100 :: UChimera Word
+-- >>> take 10 (toList ch)
+-- [100,101,103,106,110,115,121,128,136,145]
+--
+-- @since 0.3.4.0
+iterateWithIndex :: G.Vector v a => (Word -> a -> a) -> a -> Chimera v a
+iterateWithIndex f = runIdentity . iterateWithIndexM (coerce f)
+
+iterateWithIndexExactVecNM :: forall m a v. (Monad m, G.Vector v a) => Int -> (Word -> a -> m a) -> a -> m (v a)
+iterateWithIndexExactVecNM n f s = G.unfoldrExactNM n go (int2word n, s)
+  where
+    go :: (Word, a) -> m (a, (Word, a))
+    go (i, x) = do
+      x' <- f i x
+      pure (x', (i + 1, x'))
+
+-- | Monadic version of 'iterateWithIndex'.
+--
+-- @since 0.3.4.0
+iterateWithIndexM :: (Monad m, G.Vector v a) => (Word -> a -> m a) -> a -> m (Chimera v a)
+iterateWithIndexM f seed = do
+  nextSeed <- f 1 seed
+  let z = G.singleton seed
+  zs <- iterateListNM bits go (G.singleton nextSeed)
+  pure $ Chimera $ fromListN (bits + 1) (z : zs)
+  where
+    go vec =
+      iterateWithIndexExactVecNM (G.length vec `shiftL` 1) f (G.unsafeLast vec)
+{-# SPECIALIZE iterateWithIndexM :: G.Vector v a => (Word -> a -> Identity a) -> a -> Identity (Chimera v a) #-}
+
+interleaveVec :: G.Vector v a => v a -> v a -> v a
+interleaveVec as bs =
+  G.generate
+    (G.length as `shiftL` 1)
+    (\n -> (if even n then as else bs) G.! (n `shiftR` 1))
+
+-- | Intertleave two streams, sourcing even elements from the first one
+-- and odd elements from the second one.
+--
+-- >>> ch = interleave (tabulate id) (tabulate (+ 100)) :: UChimera Word
+-- >>> take 10 (toList ch)
+-- [0,100,1,101,2,102,3,103,4,104]
+--
+-- @since 0.3.3.0
+interleave :: G.Vector v a => Chimera v a -> Chimera v a -> Chimera v a
+interleave (Chimera as) (Chimera bs) = Chimera $ A.arrayFromListN (bits + 1) vecs
+  where
+    vecs =
+      A.indexArray as 0
+        : A.indexArray bs 0
+        : map (\i -> interleaveVec (A.indexArray as i) (A.indexArray bs i)) [1 .. bits - 1]
+
+-- | Index a stream in a constant time.
+--
+-- >>> ch = tabulate (^ 2) :: UChimera Word
+-- >>> index ch 9
+-- 81
+--
+-- @since 0.2.0.0
+index :: G.Vector v a => Chimera v a -> Word -> a
+index (Chimera vs) i =
+  (vs `A.indexArray` (bits - lz))
+    `G.unsafeIndex` word2int (i .&. complement ((1 `shiftL` (bits - 1)) `unsafeShiftR` lz))
+  where
+    lz :: Int
+    !lz = countLeadingZeros i
+{-# INLINE index #-}
+
+-- | Convert a stream to an infinite list.
+--
+-- >>> ch = tabulate (^ 2) :: UChimera Word
+-- >>> take 10 (toList ch)
+-- [0,1,4,9,16,25,36,49,64,81]
+--
+-- @since 0.3.0.0
+toList :: G.Vector v a => Chimera v a -> [a]
+toList (Chimera vs) = foldMap G.toList vs
+
+-- | Convert a stream to a proper infinite list.
+--
+-- @since 0.3.4.0
+toInfinite :: G.Vector v a => Chimera v a -> Infinite a
+toInfinite = foldr (:<)
+
+-- | Right-associative fold, necessarily lazy in the accumulator.
+-- Any unconditional attempt to force the accumulator even to WHNF
+-- will hang the computation. E. g., the following definition isn't productive:
+--
+-- > import Data.List.NonEmpty (NonEmpty(..))
+-- > toNonEmpty = foldr (\a (x :| xs) -> a :| x : xs) :: VChimera a -> NonEmpty a
+--
+-- One should use lazy patterns, e. g.,
+--
+-- > toNonEmpty = foldr (\a ~(x :| xs) -> a :| x : xs)
+foldr :: G.Vector v a => (a -> b -> b) -> Chimera v a -> b
+foldr f (Chimera vs) = F.foldr (flip $ G.foldr f) undefined vs
+
+measureOff :: Int -> [a] -> Either Int ([a], [a])
+measureOff n
+  | n <= 0 = Right . ([],)
+  | otherwise = go n
+  where
+    go m [] = Left m
+    go 1 (x : xs) = Right ([x], xs)
+    go m (x : xs) = case go (m - 1) xs of
+      l@Left {} -> l
+      Right (xs', xs'') -> Right (x : xs', xs'')
+
+measureOffVector :: G.Vector v a => Int -> v a -> Either Int (v a, v a)
+measureOffVector n xs
+  | n <= l = Right (G.splitAt n xs)
+  | otherwise = Left (n - l)
+  where
+    l = G.length xs
+
+-- | Create a stream of values from a given prefix, followed by default value
+-- afterwards.
+--
+-- @since 0.3.3.0
+fromListWithDef
+  :: G.Vector v a
+  => a
+  -- ^ Default value
+  -> [a]
+  -- ^ Prefix
+  -> Chimera v a
+fromListWithDef a = Chimera . fromListN (bits + 1) . go0
+  where
+    go0 = \case
+      [] -> G.singleton a : map (\k -> G.replicate (1 `shiftL` k) a) [0 .. bits - 1]
+      x : xs -> G.singleton x : go 0 xs
+
+    go k xs =
+      if k == bits
+        then []
+        else v : go (k + 1) zs
+      where
+        kk = 1 `shiftL` k
+        (v, zs) =
+          case measureOff kk xs of
+            Left l ->
+              ( if l == kk
+                  then G.replicate kk a
+                  else G.fromListN kk (xs ++ replicate l a)
+              , []
+              )
+            Right (ys, zs') -> (G.fromListN kk ys, zs')
+
+-- | Create a stream of values from a given infinite list.
+--
+-- @since 0.3.4.0
+fromInfinite
+  :: G.Vector v a
+  => Infinite a
+  -> Chimera v a
+fromInfinite = Chimera . fromListN (bits + 1) . go0
+  where
+    go0 (x :< xs) = G.singleton x : go 0 xs
+
+    go k xs =
+      if k == bits
+        then []
+        else G.fromListN kk ys : go (k + 1) zs
+      where
+        kk = 1 `shiftL` k
+        (ys, zs) = Inf.splitAt kk xs
+
+-- | Create a stream of values from a given prefix, followed by default value
+-- afterwards.
+--
+-- @since 0.3.3.0
+fromVectorWithDef
+  :: G.Vector v a
+  => a
+  -- ^ Default value
+  -> v a
+  -- ^ Prefix
+  -> Chimera v a
+fromVectorWithDef a = Chimera . fromListN (bits + 1) . go0
+  where
+    go0 xs = case G.uncons xs of
+      Nothing -> G.singleton a : map (\k -> G.replicate (1 `shiftL` k) a) [0 .. bits - 1]
+      Just (y, ys) -> G.singleton y : go 0 ys
+
+    go k xs = case measureOffVector kk xs of
+      Left l ->
+        (xs G.++ G.replicate l a)
+          : map (\n -> G.replicate (1 `shiftL` n) a) [k + 1 .. bits - 1]
+      Right (ys, zs) -> ys : go (k + 1) zs
+      where
+        kk = 1 `shiftL` k
+
+-- | Prepend a given vector to a stream of values.
+--
+-- @since 0.4.0.0
+prependVector
+  :: forall v a
+   . G.Vector v a
+  => v a
+  -> Chimera v a
+  -> Chimera v a
+prependVector (G.uncons -> Nothing) ch = ch
+prependVector (G.uncons -> Just (pref0, pref)) (Chimera as) =
+  Chimera $
+    fromListN (bits + 1) $
+      fmap sliceAndConcat $
+        [LazySlice 0 1 $ G.singleton pref0] : go 0 1 0 inputs
+  where
+    inputs :: [(Word, v a)]
+    inputs =
+      (int2word $ G.length pref, pref)
+        : zip (1 : map (1 `unsafeShiftL`) [0 .. bits - 1]) (F.toList as)
+
+    go :: Int -> Word -> Word -> [(Word, t)] -> [[LazySlice t]]
+    go _ _ _ [] = []
+    go n need off orig@((lt, t) : rest)
+      | n >= bits = []
+      | otherwise = case compare (off + need) lt of
+          LT -> [LazySlice off need t] : go (n + 1) (1 `shiftL` (n + 1)) (off + need) orig
+          EQ -> [LazySlice off need t] : go (n + 1) (1 `shiftL` (n + 1)) 0 rest
+          GT -> case go n (off + need - lt) 0 rest of
+            [] -> error "prependVector: the stream should not get exhausted prematurely"
+            hd : tl -> (LazySlice off (lt - off) t : hd) : tl
+
+data LazySlice a = LazySlice !Word !Word a
+
+sliceAndConcat :: G.Vector v a => [LazySlice (v a)] -> v a
+sliceAndConcat =
+  G.concat
+    . map (\(LazySlice from len vec) -> G.slice (word2int from) (word2int len) vec)
+
+-- | Return an infinite repetition of a given vector.
+-- Throw an error on an empty vector.
+--
+-- >>> ch = cycle (Data.Vector.fromList [4, 2]) :: VChimera Int
+-- >>> take 10 (toList ch)
+-- [4,2,4,2,4,2,4,2,4,2]
+--
+-- @since 0.3.0.0
+cycle :: G.Vector v a => v a -> Chimera v a
+cycle vec = case l of
+  0 -> error "Data.Chimera.cycle: empty list"
+  _ -> tabulate (G.unsafeIndex vec . word2int . (`rem` l))
+  where
+    l = int2word $ G.length vec
+
+-- | Map subvectors of a stream, using a given length-preserving function.
+--
+-- @since 0.3.0.0
+mapSubvectors
+  :: (G.Vector u a, G.Vector v b)
+  => (u a -> v b)
+  -> Chimera u a
+  -> Chimera v b
+mapSubvectors f = runIdentity . traverseSubvectors (coerce f)
+
+-- | Map subvectors of a stream, using a given length-preserving function.
+-- The first argument of the function is the index of the first element of subvector
+-- in the 'Chimera'.
+--
+-- @since 0.4.0.0
+imapSubvectors
+  :: (G.Vector u a, G.Vector v b)
+  => (Word -> u a -> v b)
+  -> Chimera u a
+  -> Chimera v b
+imapSubvectors f (Chimera bs) = Chimera $ mzipWith safeF (fromListN (bits + 1) [0 .. bits]) bs
+  where
+    -- Computing vector length is cheap, so let's check that @f@ preserves length.
+    safeF i x =
+      if xLen == G.length fx
+        then fx
+        else error "imapSubvectors: the function is not length-preserving"
+      where
+        xLen = G.length x
+        fx = f (if i == 0 then 0 else 1 `unsafeShiftL` (i - 1)) x
+
+-- | Traverse subvectors of a stream, using a given length-preserving function.
+--
+-- Be careful, because similar to 'tabulateM', only lazy monadic effects can
+-- be executed in a finite time: lazy state monad is fine, but strict one is
+-- not.
+--
+-- @since 0.3.3.0
+traverseSubvectors
+  :: (G.Vector u a, G.Vector v b, Applicative m)
+  => (u a -> m (v b))
+  -> Chimera u a
+  -> m (Chimera v b)
+traverseSubvectors f (Chimera bs) = Chimera <$> traverse safeF bs
+  where
+    -- Computing vector length is cheap, so let's check that @f@ preserves length.
+    safeF x =
+      ( \fx ->
+          if G.length x == G.length fx
+            then fx
+            else error "traverseSubvectors: the function is not length-preserving"
+      )
+        <$> f x
+{-# SPECIALIZE traverseSubvectors :: (G.Vector u a, G.Vector v b) => (u a -> Identity (v b)) -> Chimera u a -> Identity (Chimera v b) #-}
+
+-- | Zip subvectors from two streams, using a given length-preserving function.
+--
+-- @since 0.3.3.0
+zipWithSubvectors
+  :: (G.Vector u a, G.Vector v b, G.Vector w c)
+  => (u a -> v b -> w c)
+  -> Chimera u a
+  -> Chimera v b
+  -> Chimera w c
+zipWithSubvectors f = (runIdentity .) . zipWithMSubvectors (coerce f)
+
+-- | Zip subvectors from two streams, using a given monadic length-preserving function.
+-- Caveats for 'tabulateM' and 'traverseSubvectors' apply.
+--
+-- @since 0.3.3.0
+zipWithMSubvectors
+  :: (G.Vector u a, G.Vector v b, G.Vector w c, Applicative m)
+  => (u a -> v b -> m (w c))
+  -> Chimera u a
+  -> Chimera v b
+  -> m (Chimera w c)
+zipWithMSubvectors f (Chimera bs1) (Chimera bs2) = Chimera <$> sequenceA (mzipWith safeF bs1 bs2)
+  where
+    -- Computing vector length is cheap, so let's check that @f@ preserves length.
+    safeF x y =
+      ( \fx ->
+          if G.length x == G.length fx
+            then fx
+            else error "traverseSubvectors: the function is not length-preserving"
+      )
+        <$> f x y
+{-# SPECIALIZE zipWithMSubvectors :: (G.Vector u a, G.Vector v b, G.Vector w c) => (u a -> v b -> Identity (w c)) -> Chimera u a -> Chimera v b -> Identity (Chimera w c) #-}
+
+-- | Take a slice of 'Chimera', represented as a list on consecutive subvectors.
+--
+-- @since 0.3.3.0
+sliceSubvectors
+  :: G.Vector v a
+  => Int
+  -- ^ How many initial elements to drop?
+  -> Int
+  -- ^ How many subsequent elements to take?
+  -> Chimera v a
+  -> [v a]
+sliceSubvectors off len = doTake len . doDrop off . F.toList . unChimera
+  where
+    doTake !_ [] = []
+    doTake n (x : xs)
+      | n <= 0 = []
+      | n >= l = x : doTake (n - l) xs
+      | otherwise = [G.take n x]
+      where
+        l = G.length x
+
+    doDrop !_ [] = []
+    doDrop n (x : xs)
+      | n <= 0 = x : xs
+      | l <= n = doDrop (n - l) xs
+      | otherwise = G.drop n x : xs
+      where
+        l = G.length x
diff --git a/src/Data/Chimera/Memoize.hs b/src/Data/Chimera/Memoize.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Chimera/Memoize.hs
@@ -0,0 +1,91 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
+{-# LANGUAGE TypeFamilies #-}
+
+-- |
+-- Module:      Data.Chimera.Memoize
+-- Copyright:   (c) 2018-2019 Andrew Lelechenko
+-- Licence:     BSD3
+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
+--
+-- High-level functions for memoization.
+module Data.Chimera.Memoize (
+  memoize,
+  memoizeFix,
+) where
+
+import qualified Data.Vector as V
+import Prelude hiding (Applicative (..), and, cycle, div, drop, foldr, fromIntegral, iterate, not, or, (*), (^))
+
+import Data.Chimera.Internal
+
+-- | Memoize a function:
+-- repeating calls to 'memoize' @f@ @n@
+-- would compute @f@ @n@ only once
+-- and cache the result in 'VChimera'.
+-- This is just a shortcut for 'index' '.' 'tabulate'.
+--
+-- prop> memoize f n = f n
+--
+-- Note that @a@ could be a function type itself. This allows, for instance,
+-- to define
+--
+-- > memoize2 :: (Word -> Word -> a) -> Word -> Word -> a
+-- > memoize2 = memoize . (memoize .)
+-- >
+-- > memoize3 :: (Word -> Word -> Word -> a) -> Word -> Word -> Word -> a
+-- > memoize3 = memoize . (memoize2 .)
+--
+-- @since 0.3.0.0
+memoize :: (Word -> a) -> (Word -> a)
+memoize = index @V.Vector . tabulate
+
+-- | For a given @f@ memoize a recursive function 'Data.Function.fix' @f@,
+-- caching results in 'VChimera'.
+-- This is just a shortcut for 'index' '.' 'tabulateFix'.
+--
+-- prop> memoizeFix f n = fix f n
+--
+-- For example, imagine that we want to memoize
+-- <https://en.wikipedia.org/wiki/Fibonacci_number Fibonacci numbers>:
+--
+-- >>> fibo n = if n < 2 then toInteger n else fibo (n - 1) + fibo (n - 2)
+--
+-- Can we find @fiboF@ such that @fibo@ = 'Data.Function.fix' @fiboF@?
+-- Just replace all recursive calls to @fibo@ with @f@:
+--
+-- >>> fiboF f n = if n < 2 then toInteger n else f (n - 1) + f (n - 2)
+--
+-- Now we are ready to use 'memoizeFix':
+--
+-- >>> memoizeFix fiboF 10
+-- 55
+-- >>> memoizeFix fiboF 100
+-- 354224848179261915075
+--
+-- This function can be used even when arguments
+-- of recursive calls are not strictly decreasing,
+-- but they might not get memoized.
+-- For example, here is a routine to measure the length of
+-- <https://oeis.org/A006577 Collatz sequence>:
+--
+-- >>> collatzF f n = if n <= 1 then 0 else 1 + f (if even n then n `quot` 2 else 3 * n + 1)
+-- >>> memoizeFix collatzF 27
+-- 111
+--
+-- If you want to memoize all recursive calls, even with increasing arguments,
+-- you can employ another function of the same signature:
+-- 'Data.Function.fix' '.' ('memoize' '.'). It is less efficient though.
+--
+-- To memoize recursive functions of multiple arguments, one can use
+--
+-- > memoizeFix2 :: ((Word -> Word -> a) -> Word -> Word -> a) -> Word -> Word -> a
+-- > memoizeFix2 = let memoize2 = memoize . (memoize .) in Data.Function.fix . (memoize2 .)
+--
+-- @since 0.3.0.0
+memoizeFix :: ((Word -> a) -> Word -> a) -> (Word -> a)
+memoizeFix = index @V.Vector . tabulateFix
diff --git a/test/Test.hs b/test/Test.hs
--- a/test/Test.hs
+++ b/test/Test.hs
@@ -12,7 +12,9 @@
 import Data.Bits
 import Data.Foldable
 import Data.Function (fix)
+import qualified Data.List.Infinite as I
 import qualified Data.List as L
+import qualified Data.List.NonEmpty as NE
 import qualified Data.Vector.Generic as G
 
 import Data.Chimera.ContinuousMapping
@@ -152,11 +154,27 @@
   , QC.testProperty "toList" $
     \x xs -> xs === take (length xs) (Ch.toList (Ch.fromListWithDef x xs :: UChimera Bool))
 
-  , QC.testProperty "fromListWithDef" $
+  , testGroup "fromListWithDef"
+    [ QC.testProperty "finite list" $
+      \x xs ix ->
+        let jx = ix `mod` 65536 in
+          (if fromIntegral jx < length xs then xs !! fromIntegral jx else x) ===
+            Ch.index (Ch.fromListWithDef x xs :: UChimera Bool) jx
+
+    , QC.testProperty "infinite list" $
+      \x xs ix ->
+        let jx = ix `mod` 65536 in
+          let xs' = QC.getInfiniteList xs in
+            (xs' !! fromIntegral jx) ===
+              Ch.index (Ch.fromListWithDef x xs' :: UChimera Bool) jx
+    ]
+
+  , QC.testProperty "fromInfinite" $
     \x xs ix ->
       let jx = ix `mod` 65536 in
-        (if fromIntegral jx < length xs then xs !! fromIntegral jx else x) ===
-          Ch.index (Ch.fromListWithDef x xs :: UChimera Bool) jx
+        let ys = I.cycle (x NE.:| xs) in
+          (ys I.!! jx) ===
+            Ch.index (Ch.fromInfinite ys :: UChimera Bool) jx
 
   , QC.testProperty "fromVectorWithDef" $
     \x xs ix ->
@@ -165,12 +183,26 @@
           (if fromIntegral jx < length xs then vs G.! fromIntegral jx else x) ===
             Ch.index (Ch.fromVectorWithDef x vs :: UChimera Bool) jx
 
-  , QC.testProperty "mapWithKey" $
+  , QC.testProperty "prependVector" $
+    \(Blind bs) xs ix ->
+      let jx = ix `mod` 65536 in
+        let vs = G.fromList xs in
+          (if fromIntegral jx < length xs then vs G.! fromIntegral jx else Ch.index bs (min 65555 $ jx - fromIntegral (length xs))) ===
+            Ch.index (Ch.prependVector vs bs :: UChimera Bool) jx
+
+  , QC.testProperty "mapSubvectors" $
     \(Blind bs) (Fun _ (g :: Word -> Word)) ix ->
       let jx = ix `mod` 65536 in
         g (Ch.index bs jx) === Ch.index (Ch.mapSubvectors (G.map g) bs :: UChimera Word) jx
 
-  , QC.testProperty "zipWithKey" $
+  , QC.testProperty "imapSubvectors" $
+    \(Blind bs) (Fun _ (g :: (Word, Int) -> Char)) ix ->
+      let jx = ix `mod` 65536 in
+        curry g jx (Ch.index bs jx) ===
+          Ch.index (Ch.imapSubvectors (\off ->
+            G.imap (curry g . (+ off) . fromIntegral)) bs :: UChimera Char) jx
+
+  , QC.testProperty "zipWithSubvectors" $
     \(Blind bs1) (Blind bs2) (Fun _ (g :: (Word, Word) -> Word)) ix ->
       let jx = ix `mod` 65536 in
         g (Ch.index bs1 jx, Ch.index bs2 jx) === Ch.index (Ch.zipWithSubvectors (G.zipWith (curry g)) bs1 bs2 :: UChimera Word) jx
@@ -199,3 +231,9 @@
 
 iterateWithIndex :: (Word -> a -> a) -> a -> [a]
 iterateWithIndex f seed = L.unfoldr (\(ix, a) -> let a' = f (ix + 1) a in Just (a, (ix + 1, a'))) (0, seed)
+
+instance Arbitrary HalfWord where
+  arbitrary = fromIntegral <$> (arbitrary :: Gen Word)
+
+instance Arbitrary ThirdWord where
+  arbitrary = fromIntegral <$> (arbitrary :: Gen Word)