diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,19 @@
 # Revision history for acl-hs
 
+## 1.5.1.0 -- June 2025
+
+- Added `transformers` as a new dependency.
+- Added more functions:
+  - `AtCoder.Extra.Lct`: `read`, `lcaMaybe`, `lengthBetween` and `prodTree`.
+  - `AtCoder.Extra.Graph`: `adj1`.
+  - `AtCoder.Internal.Csr`: `adj1`.
+  - `AtCoder.Extra.Vector`:
+    - `iconcatMap`, `concatMapM`, `iconcatMapM`
+    - `compress`, `mapAccumL`, `chunks`
+    - `scanlM` variants,
+    - `maxRangeSum`, `minRangeSum`
+    - `slideMinIndices`, `slideMaxIndices`
+
 ## 1.5.0.0 -- May 2025
 
 - Changed the parameter orders of `bfs01` and `dijkstra` in `AtCoder.Extra.Graph`.
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -1,5 +1,8 @@
 # ac-library-hs
 
+[![DeepWiki](https://img.shields.io/badge/DeepWiki-toyboot4e%2Fac--library--hs-blue.svg?logo=data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAACwAAAAyCAYAAAAnWDnqAAAAAXNSR0IArs4c6QAAA05JREFUaEPtmUtyEzEQhtWTQyQLHNak2AB7ZnyXZMEjXMGeK/AIi+QuHrMnbChYY7MIh8g01fJoopFb0uhhEqqcbWTp06/uv1saEDv4O3n3dV60RfP947Mm9/SQc0ICFQgzfc4CYZoTPAswgSJCCUJUnAAoRHOAUOcATwbmVLWdGoH//PB8mnKqScAhsD0kYP3j/Yt5LPQe2KvcXmGvRHcDnpxfL2zOYJ1mFwrryWTz0advv1Ut4CJgf5uhDuDj5eUcAUoahrdY/56ebRWeraTjMt/00Sh3UDtjgHtQNHwcRGOC98BJEAEymycmYcWwOprTgcB6VZ5JK5TAJ+fXGLBm3FDAmn6oPPjR4rKCAoJCal2eAiQp2x0vxTPB3ALO2CRkwmDy5WohzBDwSEFKRwPbknEggCPB/imwrycgxX2NzoMCHhPkDwqYMr9tRcP5qNrMZHkVnOjRMWwLCcr8ohBVb1OMjxLwGCvjTikrsBOiA6fNyCrm8V1rP93iVPpwaE+gO0SsWmPiXB+jikdf6SizrT5qKasx5j8ABbHpFTx+vFXp9EnYQmLx02h1QTTrl6eDqxLnGjporxl3NL3agEvXdT0WmEost648sQOYAeJS9Q7bfUVoMGnjo4AZdUMQku50McDcMWcBPvr0SzbTAFDfvJqwLzgxwATnCgnp4wDl6Aa+Ax283gghmj+vj7feE2KBBRMW3FzOpLOADl0Isb5587h/U4gGvkt5v60Z1VLG8BhYjbzRwyQZemwAd6cCR5/XFWLYZRIMpX39AR0tjaGGiGzLVyhse5C9RKC6ai42ppWPKiBagOvaYk8lO7DajerabOZP46Lby5wKjw1HCRx7p9sVMOWGzb/vA1hwiWc6jm3MvQDTogQkiqIhJV0nBQBTU+3okKCFDy9WwferkHjtxib7t3xIUQtHxnIwtx4mpg26/HfwVNVDb4oI9RHmx5WGelRVlrtiw43zboCLaxv46AZeB3IlTkwouebTr1y2NjSpHz68WNFjHvupy3q8TFn3Hos2IAk4Ju5dCo8B3wP7VPr/FGaKiG+T+v+TQqIrOqMTL1VdWV1DdmcbO8KXBz6esmYWYKPwDL5b5FA1a0hwapHiom0r/cKaoqr+27/XcrS5UwSMbQAAAABJRU5ErkJggg==)](https://deepwiki.com/toyboot4e/ac-library-hs)
+<!-- DeepWiki badge generated by https://deepwiki.ryoppippi.com/ -->
+
 Haskell port of [ac-library](https://github.com/atcoder/ac-library), a library for competitive programming on [AtCoder](https://atcoder.jp/).
 
 The `main` branch is the version avaiable on AtCoder (+ documentation fixes). See `dev` branch for upcoming updates.
diff --git a/ac-library-hs.cabal b/ac-library-hs.cabal
--- a/ac-library-hs.cabal
+++ b/ac-library-hs.cabal
@@ -4,11 +4,11 @@
 -- PVP summary:  +-+------- breaking API changes
 --               | | +----- non-breaking API additions
 --               | | | +--- code changes with no API change
-version:         1.5.0.0
+version:         1.5.1.0
 synopsis:        Data structures and algorithms
 description:
   Haskell port of [ac-library](https://github.com/atcoder/ac-library), a library for competitive
-  programming on [AtCoder](https://atcoder.jp/).
+  programming on [AtCoder](https://atcoder.jp/). ([GitHub](https://github.com/toyboot4e/ac-library-hs), [DeepWiki](https://deepwiki.com/toyboot4e/ac-library-hs))
 
   - Functions primarily use half-open interval \([l, r)\).
   - The `Extra` module contains additional utilities beyond the original C++ library.
@@ -39,7 +39,8 @@
     , bitvec             <1.2
     , bytestring         <0.14
     , primitive          >=0.6.4.0 && <0.10
-    , random             >=1.2.0 && < 1.3
+    , transformers       >= 0.2.0.0
+    , random             >=1.2.0   && <1.3
     , vector             >=0.13.0  && <0.14
     , vector-algorithms  <0.10
     , wide-word          <0.2
@@ -170,6 +171,8 @@
     Tests.Extra.Semigroup.Permutation
     Tests.Extra.Seq
     Tests.Extra.Seq.Map
+    Tests.Extra.Tree
+    Tests.Extra.Tree.Lct
     Tests.Extra.WaveletMatrix
     Tests.Extra.WaveletMatrix.BitVector
     Tests.Extra.WaveletMatrix.Raw
@@ -230,17 +233,22 @@
     Bench.RepeatWithIndex
     Bench.RepeatWithoutIndex
     Bench.SwapDupe
+    Bench.Vector.ConcatMapM
+    Bench.Vector.IConcatMapM
+    Bench.Vector.MapAccumL
     BenchLib.AddMod
     BenchLib.Matrix
     BenchLib.ModInt.ModIntNats
     BenchLib.ModInt.Modulus
+    BenchLib.Montgomery64.Inline
+    BenchLib.Montgomery64.Noinline
     BenchLib.MulMod.Barrett64
     BenchLib.MulMod.BarrettWideWord
     BenchLib.MulMod.Montgomery
-    BenchLib.Montgomery64.Inline
-    BenchLib.Montgomery64.Noinline
     BenchLib.PowMod
     BenchLib.SwapDupe
+    BenchLib.Vector.ConcatMapM
+    BenchLib.Vector.MapAccumL
 
   build-depends:
     , ac-library-hs
diff --git a/benchmarks/Bench/Vector/ConcatMapM.hs b/benchmarks/Bench/Vector/ConcatMapM.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/Bench/Vector/ConcatMapM.hs
@@ -0,0 +1,129 @@
+module Bench.Vector.ConcatMapM (benches) where
+
+import BenchLib.Vector.ConcatMapM qualified as ConcatMapM
+import Control.Monad.Primitive (PrimMonad)
+import Control.Monad.ST (runST)
+import Control.Monad.Trans.State.Strict (StateT (..), evalState, evalStateT)
+import Control.Monad.State.Class (MonadState, modify')
+import Criterion
+import Data.Vector.Unboxed qualified as VU
+import Data.Vector.Unboxed.Mutable qualified as VUM
+import System.Random
+
+benches :: Benchmark
+benches =
+  bgroup
+    "concatMapM"
+    [ -- FIXME: too fast, something is wrong.
+      bench "list" $ nf (concatMap fList) (VU.toList vec),
+      bench "prim    IO" $ nfIO (prim1 vec),
+      bench "monad 1 IO" $ nfIO (io1 vec),
+      bench "monad 2 IO" $ nfIO (io2 vec),
+      bench "prim    ST" $ nf (\xs -> runST (prim1 xs)) vec,
+      bench "monad 1 ST" $ nf st1 vec,
+      bench "monad 2 ST" $ nf st2 vec,
+      bench "prim    StateT" $ nf (\vec -> runST $ (`evalStateT` (0 :: Int)) (ConcatMapM.primConcatMapM1 g vec)) vec,
+      bench "monad 1 State" $ nf ((`evalState` (0 :: Int)) . ConcatMapM.concatMapM1 g) vec,
+      bench "monad 2 State" $ nf ((`evalState` (0 :: Int)) . ConcatMapM.concatMapM2 g) vec,
+      bench "prime   StateT + IO" $ nfIO (evalStateT (primStateT1 vec) (0 :: Int)),
+      bench "monad 1 StateT + IO" $ nfIO (evalStateT (stateT1 vec) (0 :: Int)),
+      bench "monad 2 StateT + IO" $ nfIO (evalStateT (stateT2 vec) (0 :: Int))
+    ]
+  where
+    n = 10 ^ 3 :: Int
+    vec :: VU.Vector Int
+    vec = VU.unfoldrExactN n (uniformR (0, n - 1)) (mkStdGen (1 + 123456789))
+
+    f :: Int -> VU.Vector Int
+    f x = VU.fromList [x, x, x, x, x]
+
+    fList :: Int -> [Int]
+    fList x = [x, x, x, x, x]
+
+    g :: (MonadState Int m) => Int -> m (VU.Vector Int)
+    g x = do
+      modify' (+ 1)
+      pure $ f x
+
+    prim1 :: (PrimMonad m) => VU.Vector Int -> m (VU.Vector Int)
+    prim1 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.primConcatMapM1
+        ( \x -> do
+            VUM.modify ref (+ 1) 0
+            pure $ f x
+        )
+        vec
+
+    io1 :: VU.Vector Int -> IO (VU.Vector Int)
+    io1 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.concatMapM1
+        ( \x -> do
+            VUM.modify ref (+ 1) 0
+            pure $ f x
+        )
+        vec
+
+    io2 :: VU.Vector Int -> IO (VU.Vector Int)
+    io2 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.concatMapM2
+        ( \x -> do
+            VUM.modify ref (+ 1) 0
+            pure $ f x
+        )
+        vec
+
+    st1 :: VU.Vector Int -> VU.Vector Int
+    st1 vec = runST $ do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.concatMapM1
+        ( \x -> do
+            VUM.modify ref (+ 1) 0
+            pure $ f x
+        )
+        vec
+
+    st2 :: VU.Vector Int -> VU.Vector Int
+    st2 vec = runST $ do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.concatMapM2
+        ( \x -> do
+            VUM.modify ref (+ 1) 0
+            pure $ f x
+        )
+        vec
+
+    stateT1 :: VU.Vector Int -> StateT Int IO (VU.Vector Int)
+    stateT1 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.concatMapM1
+        ( \x -> do
+            VUM.modify ref (+ 1) 0
+            modify' (+ 1)
+            pure $ f x
+        )
+        vec
+
+    stateT2 :: VU.Vector Int -> StateT Int IO (VU.Vector Int)
+    stateT2 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.concatMapM2
+        ( \x -> do
+            VUM.modify ref (+ 1) 0
+            modify' (+ 1)
+            pure $ f x
+        )
+        vec
+
+    primStateT1 :: VU.Vector Int -> StateT Int IO (VU.Vector Int)
+    primStateT1 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.primConcatMapM1
+        ( \x -> do
+            VUM.modify ref (+ 1) 0
+            modify' (+ 1)
+            pure $ f x
+        )
+        vec
diff --git a/benchmarks/Bench/Vector/IConcatMapM.hs b/benchmarks/Bench/Vector/IConcatMapM.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/Bench/Vector/IConcatMapM.hs
@@ -0,0 +1,127 @@
+module Bench.Vector.IConcatMapM (benches) where
+
+import BenchLib.Vector.ConcatMapM qualified as ConcatMapM
+import Control.Monad.Primitive (PrimMonad)
+import Control.Monad.ST (runST)
+import Control.Monad.Trans.State.Strict (State (..), StateT (..), evalState, evalStateT)
+import Control.Monad.State.Class (MonadState, modify')
+import Criterion
+import Data.Vector.Fusion.Bundle.Monadic qualified as BundleM
+import Data.Vector.Generic qualified as VG
+import Data.Vector.Generic.Mutable qualified as VGM
+import Data.Vector.Unboxed qualified as VU
+import Data.Vector.Unboxed.Mutable qualified as VUM
+import System.Random
+
+benches :: Benchmark
+benches =
+  bgroup
+    "iconcatMapM"
+    [ bench "prim    IO" $ nfIO (prim1 vec),
+      bench "monad 1 IO" $ nfIO (io1 vec),
+      bench "monad 2 IO" $ nfIO (io2 vec),
+      bench "prim    ST" $ whnf (\xs -> runST (prim1 xs)) vec,
+      bench "monad 1 ST" $ whnf st1 vec,
+      bench "monad 2 ST" $ whnf st2 vec,
+      bench "prim    StateT" $ whnf (\vec -> runST $ (`evalStateT` (0 :: Int)) (ConcatMapM.primIConcatMapM1 g vec)) vec,
+      bench "monad 1 State" $ whnf ((`evalState` (0 :: Int)) . ConcatMapM.iconcatMapM1 g) vec,
+      bench "monad 2 State" $ whnf ((`evalState` (0 :: Int)) . ConcatMapM.iconcatMapM2 g) vec,
+      bench "prime   StateT + IO" $ nfIO (evalStateT (primStateT1 vec) (0 :: Int)),
+      bench "monad 1 StateT + IO" $ nfIO (evalStateT (stateT1 vec) (0 :: Int)),
+      bench "monad 2 StateT + IO" $ nfIO (evalStateT (stateT2 vec) (0 :: Int))
+    ]
+  where
+    n = 10 ^ 3 :: Int
+    vec :: VU.Vector Int
+    vec = VU.unfoldrExactN n (uniformR (0, n - 1)) (mkStdGen (1 + 123456789))
+
+    f :: Int -> Int -> VU.Vector Int
+    f i x = VU.fromList [i + x, i + x, i + x, i + x, i + x]
+
+    g :: (MonadState Int m) => Int -> Int -> m (VU.Vector Int)
+    g i x = do
+      modify' (+ i)
+      pure $ f i x
+
+    prim1 :: (PrimMonad m) => VU.Vector Int -> m (VU.Vector Int)
+    prim1 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.primIConcatMapM1
+        ( \i x -> do
+            VUM.modify ref (+ i) 0
+            pure $ f i x
+        )
+        vec
+
+    io1 :: VU.Vector Int -> IO (VU.Vector Int)
+    io1 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.iconcatMapM1
+        ( \i x -> do
+            VUM.modify ref (+ i) 0
+            pure $ f i x
+        )
+        vec
+
+    io2 :: VU.Vector Int -> IO (VU.Vector Int)
+    io2 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.iconcatMapM2
+        ( \i x -> do
+            VUM.modify ref (+ i) 0
+            pure $ f i x
+        )
+        vec
+
+    st1 :: VU.Vector Int -> VU.Vector Int
+    st1 vec = runST $ do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.iconcatMapM1
+        ( \i x -> do
+            VUM.modify ref (+ i) 0
+            pure $ f i x
+        )
+        vec
+
+    st2 :: VU.Vector Int -> VU.Vector Int
+    st2 vec = runST $ do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.iconcatMapM2
+        ( \i x -> do
+            VUM.modify ref (+ i) 0
+            pure $ f i x
+        )
+        vec
+
+    stateT1 :: VU.Vector Int -> StateT Int IO (VU.Vector Int)
+    stateT1 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.iconcatMapM1
+        ( \i x -> do
+            VUM.modify ref (+ i) 0
+            modify' (+ i)
+            pure $ f i x
+        )
+        vec
+
+    stateT2 :: VU.Vector Int -> StateT Int IO (VU.Vector Int)
+    stateT2 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.iconcatMapM2
+        ( \i x -> do
+            VUM.modify ref (+ i) 0
+            modify' (+ i)
+            pure $ f i x
+        )
+        vec
+
+    primStateT1 :: VU.Vector Int -> StateT Int IO (VU.Vector Int)
+    primStateT1 vec = do
+      ref <- VUM.replicate 1 (0 :: Int)
+      ConcatMapM.primIConcatMapM1
+        ( \i x -> do
+            VUM.modify ref (+ i) 0
+            modify' (+ i)
+            pure $ f i x
+        )
+        vec
diff --git a/benchmarks/Bench/Vector/MapAccumL.hs b/benchmarks/Bench/Vector/MapAccumL.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/Bench/Vector/MapAccumL.hs
@@ -0,0 +1,26 @@
+module Bench.Vector.MapAccumL (benches) where
+
+import BenchLib.Vector.MapAccumL qualified as MapAccumL
+import Criterion
+import Data.List qualified as L
+import Data.Vector.Unboxed qualified as VU
+import System.Random
+
+benches :: Benchmark
+benches =
+  bgroup
+    "mapAccumL"
+    [ -- whnf did not work (somehow) for mapAccumL1, so I'm using nf
+      -- NOTE: list is not fair to compare though
+      bench "list" $ nf (L.mapAccumL f (0 :: Int)) (VU.toList vec),
+      bench "mapM + State" $ nf (MapAccumL.mapAccumL1 f (0 :: Int)) vec,
+      bench "ifoldM'" $ nf (MapAccumL.mapAccumL2 f (0 :: Int)) vec,
+      bench "bundle" $ nf (MapAccumL.mapAccumL3 f (0 :: Int)) vec,
+      bench "mapM + State + PrimMonad" $ nf (MapAccumL.mapAccumL4 f (0 :: Int)) vec
+    ]
+  where
+    n = 10 ^ 3 :: Int
+    vec :: VU.Vector Int
+    vec = VU.unfoldrExactN n (uniformR (0, n - 1)) (mkStdGen (1 + 123456789))
+    f :: Int -> Int -> (Int, Int)
+    f s x = (s + 10, s * x)
diff --git a/benchmarks/BenchLib/Vector/ConcatMapM.hs b/benchmarks/BenchLib/Vector/ConcatMapM.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/BenchLib/Vector/ConcatMapM.hs
@@ -0,0 +1,101 @@
+{-# LANGUAGE RecordWildCards #-}
+
+-- | Fast modular multiplication by barrett reduction.
+module BenchLib.Vector.ConcatMapM
+  ( primConcatMapM1,
+    concatMapM1,
+    concatMapM2,
+    primIConcatMapM1,
+    iconcatMapM1,
+    iconcatMapM2,
+  )
+where
+
+import Control.Monad.Primitive (PrimMonad)
+import Data.Vector.Fusion.Bundle.Monadic qualified as BundleM
+import Data.Vector.Generic qualified as VG
+import Data.Vector.Generic.Mutable qualified as VGM
+
+-- I expected these functions do not perform well with State monad, but
+-- they're working as fast as IO and ST?
+
+-- | https://github.com/haskell/vector/issues/416
+--
+-- With this wrapper
+unstreamPrimM :: (PrimMonad m, VG.Vector v a) => BundleM.Bundle m u a -> m (v a)
+{-# INLINE [1] unstreamPrimM #-}
+unstreamPrimM s = VGM.munstream s >>= VG.unsafeFreeze
+
+-- | PrimMonad only. This addes more constraints, but performs well with complex monads.
+{-# INLINE primConcatMapM1 #-}
+primConcatMapM1 ::
+  (PrimMonad m, VG.Vector v a, VG.Vector v b) =>
+  (a -> m (v b)) ->
+  v a ->
+  m (v b)
+primConcatMapM1 f =
+  unstreamPrimM
+    . BundleM.concatVectors
+    . BundleM.mapM f
+    . BundleM.fromVector
+
+{-# INLINE concatMapM1 #-}
+concatMapM1 ::
+  (Monad m, VG.Vector v a, VG.Vector v b) =>
+  (a -> m (v b)) ->
+  v a ->
+  m (v b)
+concatMapM1 f =
+  VG.unstreamM
+    . BundleM.concatVectors
+    . BundleM.mapM f
+    . BundleM.fromVector
+
+{-# INLINE concatMapM2 #-}
+concatMapM2 ::
+  (Monad m, VG.Vector v a, VG.Vector v b) =>
+  (a -> m (v b)) ->
+  v a ->
+  m (v b)
+concatMapM2 f xs =
+  VG.unstreamM
+    . BundleM.concatVectors
+    $ BundleM.generateM (VG.length xs) (f . VG.unsafeIndex xs)
+
+-- | PrimMonad only. This addes more constraints, but performs well with complex monads.
+{-# INLINE primIConcatMapM1 #-}
+primIConcatMapM1 ::
+  (PrimMonad m, VG.Vector v a, VG.Vector v b) =>
+  (Int -> a -> m (v b)) ->
+  v a ->
+  m (v b)
+primIConcatMapM1 f =
+  unstreamPrimM
+    . BundleM.concatVectors
+    . BundleM.mapM (uncurry f)
+    . BundleM.indexed
+    . BundleM.fromVector
+
+{-# INLINE iconcatMapM1 #-}
+iconcatMapM1 ::
+  (Monad m, VG.Vector v a, VG.Vector v b) =>
+  (Int -> a -> m (v b)) ->
+  v a ->
+  m (v b)
+iconcatMapM1 f =
+  VG.unstreamM
+    . BundleM.concatVectors
+    . BundleM.mapM (uncurry f)
+    . BundleM.indexed
+    . BundleM.fromVector
+
+{-# INLINE iconcatMapM2 #-}
+iconcatMapM2 ::
+  (Monad m, VG.Vector v a, VG.Vector v b) =>
+  (Int -> a -> m (v b)) ->
+  v a ->
+  m (v b)
+iconcatMapM2 f xs =
+  VG.unstreamM
+    . BundleM.concatVectors
+    $ BundleM.generateM (VG.length xs) (\i -> f i (VG.unsafeIndex xs i))
diff --git a/benchmarks/BenchLib/Vector/MapAccumL.hs b/benchmarks/BenchLib/Vector/MapAccumL.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/BenchLib/Vector/MapAccumL.hs
@@ -0,0 +1,111 @@
+-- | Fast modular multiplication by barrett reduction.
+module BenchLib.Vector.MapAccumL
+  ( mapAccumL1,
+    mapAccumL2,
+    mapAccumL3,
+    mapAccumL4,
+  )
+where
+
+import Control.Monad.Primitive (PrimMonad)
+import Control.Monad.ST (ST, runST)
+import Control.Monad.Trans.State.Strict (StateT, runState, runStateT, state)
+import Data.Vector.Fusion.Bundle qualified as Bundle
+import Data.Vector.Fusion.Bundle.Monadic qualified as BundleM
+import Data.Vector.Fusion.Bundle.Size qualified as Bundle
+import Data.Vector.Fusion.Stream.Monadic qualified as MS
+import Data.Vector.Generic qualified as VG
+import Data.Vector.Generic.Mutable qualified as VGM
+
+{-# INLINE mapAccumL1 #-}
+mapAccumL1 ::
+  (VG.Vector v a, VG.Vector v b) =>
+  (s -> a -> (s, b)) ->
+  s ->
+  v a ->
+  (s, v b)
+mapAccumL1 f s0 xs = (\(!x, !s) -> (s, x)) $ runState (VG.mapM g xs) s0
+  where
+    g a =
+      state
+        ( \s ->
+            let (!s', !b) = f s a
+             in (b, s')
+        )
+
+{-# INLINE mapAccumL2 #-}
+mapAccumL2 ::
+  (VG.Vector v a, VG.Vector v b) =>
+  (s -> a -> (s, b)) ->
+  s ->
+  v a ->
+  (s, v b)
+mapAccumL2 f s0 xs = runST $ do
+  vec <- VGM.unsafeNew (VG.length xs)
+  !res <-
+    VG.ifoldM'
+      ( \ !s i x -> do
+          let (!s', !x') = f s x
+          VGM.write vec i x'
+          pure s'
+      )
+      s0
+      xs
+  vec' <- VG.unsafeFreeze vec
+  pure (res, vec')
+
+{-# INLINE mapAccumL3 #-}
+mapAccumL3 ::
+  forall v s a b.
+  (VG.Vector v a, VG.Vector v b) =>
+  (s -> a -> (s, b)) ->
+  s ->
+  v a ->
+  (s, v b)
+mapAccumL3 f s0 xs = (\(!x, !s) -> (s, x)) $ runST $ (`runStateT` s0) $ do
+  unstreamPrimM
+    . BundleM.mapM g
+    $ BundleM.fromVector xs
+  where
+    g :: forall st. a -> StateT s (ST st) b
+    g a =
+      state
+        ( \s ->
+            let (!s', !b) = f s a
+             in (b, s')
+        )
+
+-- | https://github.com/haskell/vector/issues/416
+unstreamPrimM :: (PrimMonad m, VG.Vector v a) => BundleM.Bundle m u a -> m (v a)
+{-# INLINE [1] unstreamPrimM #-}
+unstreamPrimM s = VGM.munstream s >>= VG.unsafeFreeze
+
+-- @cojna/iota. Note that it does not return state value.
+{-# INLINE mapAccumL4 #-}
+mapAccumL4 ::
+  (VG.Vector v a, VG.Vector v b) =>
+  (s -> a -> (s, b)) ->
+  s ->
+  v a ->
+  v b
+mapAccumL4 f x =
+  VG.unstream
+    . Bundle.inplace
+      (streamAccumM (\s a -> pure (f s a)) x)
+      Bundle.toMax
+    . VG.stream
+
+-- @cojna/iota
+streamAccumM :: (Monad m) => (s -> a -> m (s, b)) -> s -> MS.Stream m a -> MS.Stream m b
+streamAccumM f s0 (MS.Stream step x0) = MS.Stream step' (s0, x0)
+  where
+    step' (!s, x) = do
+      r <- step x
+      case r of
+        MS.Yield a x' -> do
+          (s', b) <- f s a
+          return $ MS.Yield b (s', x')
+        MS.Skip x' -> return $ MS.Skip (s, x')
+        MS.Done -> return MS.Done
+    {-# INLINE [0] step' #-}
+{-# INLINE [1] streamAccumM #-}
diff --git a/benchmarks/Main.hs b/benchmarks/Main.hs
--- a/benchmarks/Main.hs
+++ b/benchmarks/Main.hs
@@ -1,5 +1,3 @@
-{-# LANGUAGE CPP #-}
-
 module Main where
 
 import Bench.AddMod qualified
@@ -11,12 +9,15 @@
 import Bench.RepeatWithIndex qualified
 import Bench.RepeatWithoutIndex qualified
 import Bench.SwapDupe qualified
+import Bench.Vector.ConcatMapM qualified
+import Bench.Vector.IConcatMapM qualified
+import Bench.Vector.MapAccumL qualified
 import Criterion.Main
 
 main :: IO ()
 main =
   defaultMain
-    -- TODO: generate criterion graph by benchmark group?
+    -- TODO: Filter benchmark cases with command line arguments. Currently I'm filtering benchmarks by hand.
     [ Bench.MulMod.benches,
       Bench.ModInt.benches,
       Bench.AddMod.benches,
@@ -25,5 +26,8 @@
       Bench.Montgomery64.benches,
       Bench.RepeatWithIndex.benches,
       Bench.RepeatWithoutIndex.benches,
-      Bench.SwapDupe.benches
+      Bench.SwapDupe.benches,
+      Bench.Vector.ConcatMapM.benches,
+      Bench.Vector.IConcatMapM.benches,
+      Bench.Vector.MapAccumL.benches
     ]
diff --git a/src/AtCoder/Extra/Bisect.hs b/src/AtCoder/Extra/Bisect.hs
--- a/src/AtCoder/Extra/Bisect.hs
+++ b/src/AtCoder/Extra/Bisect.hs
@@ -10,7 +10,7 @@
 -- @
 --
 -- ==== __Example__
--- Perform index compression:
+-- Perform index compression with `lowerBound`:
 --
 -- >>> import AtCoder.Extra.Bisect
 -- >>> import Data.Vector.Algorithms.Intro qualified as VAI
@@ -22,13 +22,13 @@
 --
 -- @since 1.3.0.0
 module AtCoder.Extra.Bisect
-  ( -- * C++-like binary search
+  ( -- * C++-like binary searches
     lowerBound,
     lowerBoundIn,
     upperBound,
     upperBoundIn,
 
-    -- * Generic bisection method
+    -- * Generic bisection methods
     maxRight,
     maxRightM,
     minLeft,
@@ -45,7 +45,7 @@
 --
 -- @
 -- Y Y Y Y Y N N N N N      Y: x_i < x_ref
--- --------- *---------> x  N: x_i >= x_ref
+-- --------- *---------> x  N: not Y
 --           R              R: the right boundary point returned
 -- @
 --
@@ -106,8 +106,8 @@
 -- | \(O(\log n)\) Returns the maximum \(r\) where \(x_i \le x_{ref}\) holds for \(i \in [0, r)\).
 --
 -- @
--- Y Y Y Y Y N N N N N      Y: x_i <= x_ref,
--- --------- *---------> x  N: x_i > x_ref,
+-- Y Y Y Y Y N N N N N      Y: x_i <= x_ref
+-- --------- *---------> x  N: not Y
 --           R              R: the right boundary point returned
 -- @
 --
@@ -175,8 +175,8 @@
 -- right boundary point.
 --
 -- @
--- Y Y Y Y Y N N N N N      Y: p(i) returns `true`,
--- --------- *---------> x  N: p(i) returns `false`,
+-- Y Y Y Y Y N N N N N      Y: p(i) returns `true`
+-- --------- *---------> x  N: not Y
 --           R              R: the right boundary point returned
 -- @
 --
@@ -223,8 +223,8 @@
 -- left boundary point.
 --
 -- @
--- N N N N N Y Y Y Y Y      Y: p(i) returns `true`,
--- --------* ----------> x  N: p(i) returns `false`,
+-- N N N N N Y Y Y Y Y      Y: p(i) returns `true`
+-- --------* ----------> x  N: not Y
 --         L                L: the left boundary point returned
 -- @
 --
diff --git a/src/AtCoder/Extra/Graph.hs b/src/AtCoder/Extra/Graph.hs
--- a/src/AtCoder/Extra/Graph.hs
+++ b/src/AtCoder/Extra/Graph.hs
@@ -1686,7 +1686,7 @@
   where
     idx !from !to = nVerts * from + to
 
--- | \(O(n)\) Given a predecessor array, retrieves a path from the root to a vertex.
+-- | \(O(n)\) Given a predecessor array, reconstructs a path from the root to a vertex.
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
@@ -1698,7 +1698,7 @@
 constructPathFromRoot :: (HasCallStack) => VU.Vector Int -> Int -> VU.Vector Int
 constructPathFromRoot parents = VU.reverse . constructPathToRoot parents
 
--- | \(O(n)\) Given a predecessor array, retrieves a path from a vertex to the root.
+-- | \(O(n)\) Given a predecessor array, reconstructs a path from a vertex to the root.
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
@@ -1713,7 +1713,7 @@
     f (-1) = Nothing
     f v = Just (v, parents VG.! v)
 
--- | \(O(n)\) Given a NxN predecessor matrix (created with `trackingFloydWarshall`), retrieves a
+-- | \(O(n)\) Given a NxN predecessor matrix (created with `trackingFloydWarshall`), reconstructs a
 -- path from the root to a sink vertex.
 --
 -- ==== Constraints
@@ -1735,7 +1735,7 @@
   VU.Vector Int
 constructPathFromRootMat parents source = VU.reverse . constructPathToRootMat parents source
 
--- | \(O(n)\) Given a NxN predecessor matrix(created with `trackingFloydWarshall`), retrieves a
+-- | \(O(n)\) Given a NxN predecessor matrix(created with `trackingFloydWarshall`), reconstructs a
 -- path from a vertex to the root.
 --
 -- ==== Constraints
@@ -1763,8 +1763,8 @@
     -- https://zenn.dev/mod_poppo/articles/atcoder-beginner-contest-284-d#%E8%A7%A3%E6%B3%953%EF%BC%9Asqrt%E3%81%A8round%E3%82%92%E4%BD%BF%E3%81%86
     n :: Int = round . sqrt $ (fromIntegral (VU.length parents) :: Double)
 
--- | \(O(n)\) Given a NxN predecessor matrix (created with `newTrackingFloydWarshall`), retrieves a
--- path from the root to a sink vertex.
+-- | \(O(n)\) Given a NxN predecessor matrix (created with `newTrackingFloydWarshall`), reconstructs
+-- a path from the root to a sink vertex.
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
@@ -1785,8 +1785,8 @@
   m (VU.Vector Int)
 constructPathFromRootMatM parents source = (VU.reverse <$>) . constructPathToRootMatM parents source
 
--- | \(O(n)\) Given a NxN predecessor matrix (created with `newTrackingFloydWarshall`), retrieves a
--- path from a vertex to the root.
+-- | \(O(n)\) Given a NxN predecessor matrix (created with `newTrackingFloydWarshall`),
+-- reconstructs a path from a vertex to the root.
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
diff --git a/src/AtCoder/Extra/Math/Montgomery64.hs b/src/AtCoder/Extra/Math/Montgomery64.hs
--- a/src/AtCoder/Extra/Math/Montgomery64.hs
+++ b/src/AtCoder/Extra/Math/Montgomery64.hs
@@ -84,7 +84,7 @@
   let !m128 :: Word128 = fromIntegral m
       !n2 = word128Lo64 $ (-m128) `mod` m128
       !r = getR m 0
-      !_ = ACIA.runtimeAssert (r * m == -1) "AtCoder.Extra.Montgomery64.fromVal: internal implementation error"
+      !_ = ACIA.runtimeAssert (r * m == maxBound) "AtCoder.Extra.Montgomery64.fromVal: internal implementation error"
    in Montgomery64 m r n2
   where
     !_ = ACIA.runtimeAssert (odd m && m <= bit 62) $ "AtCoder.Extra.Montgomery64.fromVal: not given odd modulus value that is less than or equal to 2^62: " ++ show m
diff --git a/src/AtCoder/Extra/SegTree2d/Dense.hs b/src/AtCoder/Extra/SegTree2d/Dense.hs
--- a/src/AtCoder/Extra/SegTree2d/Dense.hs
+++ b/src/AtCoder/Extra/SegTree2d/Dense.hs
@@ -180,7 +180,7 @@
 -- @since 1.2.3.0
 {-# INLINE write #-}
 write :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => DenseSegTree2d (PrimState m) a -> Int -> Int -> a -> m ()
-write seg@DenseSegTree2d {..} x y a = stToPrim $ do
+write seg x y a = stToPrim $ do
   modifyM seg (pure . const a) x y
 
 -- | \(O(\log  h \log w)\) Given a user function \(f\), modifies the monoid value at \((x, y)\) with
diff --git a/src/AtCoder/Extra/Semigroup/Permutation.hs b/src/AtCoder/Extra/Semigroup/Permutation.hs
--- a/src/AtCoder/Extra/Semigroup/Permutation.hs
+++ b/src/AtCoder/Extra/Semigroup/Permutation.hs
@@ -61,6 +61,9 @@
 
 -- | \(O(1)\) Creates a `Permutation`, performing boundary check on input vector.
 --
+-- ==== Constraints
+-- - \(-1 \le x \lt n\)
+--
 -- @since 1.1.0.0
 {-# INLINE new #-}
 new :: (HasCallStack) => VU.Vector Int -> Permutation
@@ -91,7 +94,7 @@
 zero :: Int -> Permutation
 zero n = Permutation $ VU.replicate n (-1)
 
--- | \(O(1)\) Maps an index.
+-- | \(O(1)\) Maps an index with the permutation.
 --
 -- @since 1.1.0.0
 {-# INLINE act #-}
@@ -112,6 +115,6 @@
   {-# INLINE (<>) #-}
   Permutation r2 <> Permutation r1 = Permutation $ VU.map f r1
     where
-      !_ = ACIA.runtimeAssert (VU.length r2 == VU.length r1) "AtCoder.Extra.Semigroup.Permutation.(<>): legth mismatch"
+      !_ = ACIA.runtimeAssert (VU.length r2 == VU.length r1) "AtCoder.Extra.Semigroup.Permutation.(<>): length mismatch"
       f (-1) = -1
       f i = VG.unsafeIndex r2 i
diff --git a/src/AtCoder/Extra/Tree.hs b/src/AtCoder/Extra/Tree.hs
--- a/src/AtCoder/Extra/Tree.hs
+++ b/src/AtCoder/Extra/Tree.hs
@@ -101,9 +101,10 @@
       !w = bfs2 VG.! to
    in (Gr.constructPathFromRoot parents to, w)
 
--- | \(O(m \log m)\) Kruskal's algorithm. Returns edge indices for building a minimum spanning tree.
+-- | \(O(m \log m)\) Kruskal's algorithm. Returns \((w, \mathrm{use}, \mathrm{mst})\): (weight of the minimum
+-- spanning tree, edge use flags, minimum spanning tree).
 --
--- NOTE: The edges should not be duplicated: only one of \((u, v, w)\) or \((v, u, w)\) is required
+-- NOTE: The edges need not be duplicated: only one of \((u, v, w)\) or \((v, u, w)\) is required
 -- for each edge.
 --
 -- ==== __Example__
@@ -127,10 +128,12 @@
 mst :: (Num w, Ord w, VU.Unbox w) => Int -> VU.Vector (Int, Int, w) -> (w, VU.Vector Bit, Gr.Csr w)
 mst = mstBy (comparing id)
 
--- | \(O(m \log m)\) Kruskal's algorithm. Returns edge indices for building a minimum/maximum
--- spanning tree.
+-- | \(O(m \log m)\) Kruskal's algorithm for creating either minimum or maximum spanning tree. Returns
+-- \((w, \mathrm{use}, \mathrm{mst})\): (weight of the minimum spanning tree, edge use flags, minimum
+-- spanning tree).
 --
--- NOTE: The edges should not be duplicated: only one of \((u, v, w)\) or \((v, u, w)\) is required
+--
+-- NOTE: The edges need not be duplicated: only one of \((u, v, w)\) or \((v, u, w)\) is required
 -- for each edge.
 --
 -- ==== __Example__
diff --git a/src/AtCoder/Extra/Tree/Hld.hs b/src/AtCoder/Extra/Tree/Hld.hs
--- a/src/AtCoder/Extra/Tree/Hld.hs
+++ b/src/AtCoder/Extra/Tree/Hld.hs
@@ -84,6 +84,9 @@
 -- >>> Hld.ancestor hld 5 3 -- go up three parents from `5`
 -- 0
 --
+-- >>> Hld.lca hld 2 5 -- lowest common ancestor of `2` and `5`:
+-- 1
+--
 -- >>> Hld.jump hld 5 2 3   -- go to the third vertex from `5` to `2`:
 -- Just 2
 --
@@ -93,7 +96,8 @@
 -- >>> Hld.path hld 5 3     -- get the path between `5` and `3`:
 -- [5,4,1,2,3]
 --
--- Our `Hld` is rooted at @0@ vertex and subtree queries are available:
+-- So the subtree queries are run with recpect to root vertex `0`, since our HLD is created with
+-- `new`:
 --
 -- >>> Hld.isInSubtree hld 2 3 -- `3` is in the subtree of `2`
 -- True
@@ -218,7 +222,7 @@
 --
 -- ==== Constraints
 -- - \(n \ge 1\)
--- - The input graph must be a tree; the edges must be undirected (both \((u, v, w)\) and
+-- - The input graph must be a tree. Note that the edges must be undirected (both \((u, v, w)\) and
 --   \((v, u, w)\) edges are required).
 --
 -- @since 1.1.0.0
@@ -230,7 +234,7 @@
 --
 -- ==== Constraints
 -- - \(n \ge 1\)
--- - The input graph must be a tree; the edges must be undirected (both \((u, v, w)\) and
+-- - The input graph must be a tree. Note that the edges must be undirected (both \((u, v, w)\) and
 --   \((v, u, w)\) edges are required).
 --
 -- @since 1.1.0.0
diff --git a/src/AtCoder/Extra/Tree/Lct.hs b/src/AtCoder/Extra/Tree/Lct.hs
--- a/src/AtCoder/Extra/Tree/Lct.hs
+++ b/src/AtCoder/Extra/Tree/Lct.hs
@@ -2,11 +2,16 @@
 {-# LANGUAGE RecordWildCards #-}
 {-# LANGUAGE TypeFamilies #-}
 
--- | Link/cut tree: forest with monoid values.
+-- | Link/cut tree: dynamic forest with monoid values on vertices. If you need to have monoid values
+-- on edges, treat the original edges as new vertices.
 --
+-- - __Most operations are unsafe__; user must ensure connectivities of \(u\) and \(v\) before
+--   running each query.
+-- - This specific implementation is not capable of applying monoid action to a subtree.
+--
 -- ==== __Example__
 --
--- Create a link/cut tree of @Sum Int@ with inverse operator `negate`:
+-- Create a link/cut tree of monoid @Sum Int@ with inverse operator `negate`:
 --
 -- >>> import AtCoder.Extra.Tree.Lct qualified as Lct
 -- >>> import Data.Semigroup (Sum (..))
@@ -15,11 +20,14 @@
 -- >>> --    +--3
 -- >>> lct <- Lct.buildInv negate (VU.generate 4 Sum) $ VU.fromList [(0, 1), (1, 2), (1, 3)]
 --
+-- ===== `prodPath`, `prodSubtree`, `prodTree`
+--
 -- Monoid products can be calculated for paths or subtrees:
 --
 -- >>> Lct.prodPath lct 0 2
 -- Sum {getSum = 3}
 --
+-- >>> -- If we create the LCT with `buildInv`, we can use `prodSubtree`:
 -- >>> Lct.prodSubtree lct 1 {- parent -} 2
 -- Sum {getSum = 4}
 --
@@ -28,6 +36,8 @@
 -- >>> Lct.root lct 3
 -- 2
 --
+-- ===== `lca`, `jump`, `lengthBetween`
+--
 -- Set (`evert`) the root of the underlying tree to \(0\) and get the `lca` of vertices \(2\) and
 -- \(3\):
 --
@@ -37,12 +47,27 @@
 --
 -- Similar to @Hld@, `Lct` allows various tree queries:
 --
--- >>> Lct.parent lct 3
--- Just 1
---
--- >>> Lct.jump lct 2 3 2
+-- >>> Lct.jump lct {- path -} 2 3 {- k -} 2
 -- 3
 --
+-- >>> Lct.jumpMaybe lct {- path -} 2 3 {- k -} 1000
+-- Nothing
+--
+-- >>> Lct.lengthBetween lct {- path -} 2 3
+-- 2
+--
+-- ===== `parent`
+--
+-- >>> Lct.evert lct 0  -- set root `0`
+-- >>> Lct.parent lct 0 -- under root `0`, parent of `0` is `Nothing`:
+-- Nothing
+--
+-- >>> Lct.evert lct 0  -- set root `0`
+-- >>> Lct.parent lct 1 -- under root `0`, parent of `1` is `0`:
+-- Just 0
+--
+-- ===== `link` / `cut`
+--
 -- Edges can be dynamically added (`link`) or removed (`cut`):
 --
 -- >>> -- 0  1  2
@@ -71,13 +96,14 @@
     build,
     buildInv,
 
-    -- * Modifications
-
-    -- ** Write
+    -- * Monoid value access
+    read,
     write,
     modify,
     modifyM,
 
+    -- * Tree operations
+
     -- ** Link/cut
     link,
     cut,
@@ -88,16 +114,19 @@
     expose_,
 
     -- * Tree queries
-
-    -- ** Root, parent, jump, LCA
     root,
+    same,
     parent,
     jump,
+    jumpMaybe,
+    lengthBetween,
     lca,
+    lcaMaybe,
 
-    -- ** Products
+    -- ** Monoid products
     prodPath,
     prodSubtree,
+    prodTree,
   )
 where
 
@@ -107,10 +136,12 @@
 import Control.Monad.ST (ST)
 import Data.Bit
 import Data.Bits
+import Data.Maybe
 import Data.Vector.Generic.Mutable qualified as VGM
 import Data.Vector.Unboxed qualified as VU
 import Data.Vector.Unboxed.Mutable qualified as VUM
 import GHC.Stack (HasCallStack)
+import Prelude hiding (read)
 
 -- import GHC.Stack (HasCallStack)
 
@@ -151,7 +182,7 @@
     --
     -- @since 1.1.1.0
     sLct :: !(VUM.MVector s Int),
-    -- | Decomposed node data storage: reverse flag.
+    -- | Decomposed node data storage: reverse flags.
     --
     -- @since 1.1.1.0
     revLct :: !(VUM.MVector s Bit),
@@ -163,12 +194,12 @@
     --
     -- @since 1.1.1.0
     prodLct :: !(VUM.MVector s a),
-    -- | Decomposed node data storage: dual monod product (right fold). This is required for
+    -- | Decomposed node data storage: dual monoid products (right foldings). This is required for
     -- non-commutative monoids only.
     --
     -- @since 1.1.1.0
     dualProdLct :: !(VUM.MVector s a),
-    -- | Decomposed node data storage: path-parent monoid product. This works for subtree product
+    -- | Decomposed node data storage: path-parent monoid products. This works for subtree product
     -- queries over commutative monoids only.
     --
     -- @since 1.1.1.0
@@ -193,8 +224,8 @@
 new :: (PrimMonad m, Monoid a, VU.Unbox a) => Int -> m (Lct (PrimState m) a)
 new = newInv id
 
--- | \(O(n + m \log n)\) Creates a link/cut tree with an inverse operator, initial monoid values and
--- no edges. This setup enables subtree queries (`prodSubtree`).
+-- | \(O(n)\) Creates a link/cut tree with an inverse operator, initial monoid values and no edges.
+-- This setup enables subtree queries (`prodSubtree`).
 --
 -- @since 1.1.1.0
 {-# INLINE newInv #-}
@@ -204,6 +235,9 @@
 -- | \(O(n + m \log n)\) Creates a link/cut tree of initial monoid values and initial edges. This
 -- setup disables subtree queries (`prodSubtree`).
 --
+-- ==== Constraints
+-- - \(0 \le u, v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE build #-}
 build ::
@@ -219,6 +253,9 @@
 -- | \(O(n + m \log n)\) Creates a link/cut tree with an inverse operator, initial monoid values and
 -- initial edges. This setup enables subtree queries (`prodSubtree`).
 --
+-- ==== Constraints
+-- - \(0 \le u, v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE buildInv #-}
 buildInv ::
@@ -234,46 +271,68 @@
 buildInv invOpLct xs es = stToPrim $ buildST invOpLct xs es
 
 -- -------------------------------------------------------------------------------------------------
--- Write
+-- Monoid value access
 -- -------------------------------------------------------------------------------------------------
 
--- | Amortized \(O(\log n)\). Writes the monoid value of a vertex.
+-- | \(O(1)\). Reads the monoid value on a vertex \(v\).
 --
+-- ==== Constraints
+-- - \(0 \le v \lt n\)
+--
+-- @since 1.5.1.0
+{-# INLINE read #-}
+read :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> m a
+read lct v = stToPrim $ do
+  VGM.read (vLct lct) v
+  where
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.read" v (nLct lct)
+
+-- | Amortized \(O(\log n)\). Writes to the monoid value of a vertex \(v\).
+--
+-- ==== Constraints
+-- - \(0 \le v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE write #-}
 write :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> a -> m ()
 write lct v x = stToPrim $ do
-  -- make @v@ the new root of the underlying tree:
+  -- make @v@ a new root of the underlying tree:
   evertST lct v
   VGM.unsafeWrite (vLct lct) v x
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.write" v (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.write" v (nLct lct)
 
 -- | Amortized \(O(\log n)\). Given a user function \(f\), modifies the monoid value of a vertex
 -- \(v\).
 --
+-- ==== Constraints
+-- - \(0 \le v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE modify #-}
 modify :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> (a -> a) -> Vertex -> m ()
 modify lct f v = stToPrim $ do
-  -- make @v@ the new root of the underlying tree:
+  -- make @v@ a new root of the underlying tree:
   evertST lct v
   VGM.unsafeModify (vLct lct) f v
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.modify" v (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.modify" v (nLct lct)
 
 -- | Amortized \(O(\log n)\). Given a user function \(f\), modifies the monoid value of a vertex
 -- \(v\).
 --
+-- ==== Constraints
+-- - \(0 \le v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE modifyM #-}
 modifyM :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> (a -> m a) -> Vertex -> m ()
 modifyM lct f v = do
-  -- make @v@ the new root of the underlying tree:
+  -- make @v@ a new root of the underlying tree:
   stToPrim $ evertST lct v
   VGM.unsafeModifyM (vLct lct) f v
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.modifyM" v (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.modifyM" v (nLct lct)
 
 -- -------------------------------------------------------------------------------------------------
 -- Link/cut operations
@@ -282,99 +341,179 @@
 -- | Amortized \(O(\log n)\). Creates an edge between \(c\) and \(p\). In the represented tree, the
 -- \(p\) will be the parent of \(c\).
 --
+-- ==== Constraints
+-- - \(0 \le c, p \lt n\)
+-- - \(u \ne v\)
+-- - \(c\) and \(p\) are in different trees, otherwise the behavior is undefined.
+--
 -- @since 1.1.1.0
 {-# INLINE link #-}
 link :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> Vertex -> m ()
 link lct c p = stToPrim $ linkST lct c p
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.link" c (nLct lct)
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.link" p (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.link" c (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.link" p (nLct lct)
 
 -- | Amortized \(O(\log n)\). Deletes an edge between \(u\) and \(v\).
 --
+-- ==== Constraints
+-- - \(0 \le u, v \lt n\)
+-- - \(u \ne v\)
+-- - There's an edge between \(u\) and \(v\), otherwise the behavior is undefined.
+--
 -- @since 1.1.1.0
 {-# INLINE cut #-}
 cut :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> Vertex -> m ()
 cut lct u v = stToPrim $ cutST lct u v
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.cut" u (nLct lct)
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.cut" v (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.cut" u (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.cut" v (nLct lct)
 
 -- | Amortized \(O(\log n)\). Makes \(v\) a new root of the underlying tree.
 --
+-- ==== Constraints
+-- - \(0 \le v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE evert #-}
 evert :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> m ()
 evert lct v = stToPrim $ evertST lct v
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.evert" v (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.evert" v (nLct lct)
 
 -- | Amortized \(O(\log n)\). Makes \(v\) and the root to be in the same preferred path (auxiliary
--- tree). After the opeartion, \(v\) will be the new root and all the children will be detached from
+-- tree). After the operation, \(v\) will be the new root and all the children will be detached from
 -- the preferred path.
 --
+-- ==== Constraints
+-- - \(0 \le v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE expose #-}
 expose :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> m Vertex
 expose lct v = stToPrim $ exposeST lct v
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.expose_" v (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.expose" v (nLct lct)
 
 -- | Amortized \(O(\log n)\). `expose` with the return value discarded.
 --
+-- ==== Constraints
+-- - \(0 \le v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE expose_ #-}
 expose_ :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> m ()
 expose_ lct v0 = stToPrim $ do
   _ <- exposeST lct v0
   pure ()
-  where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.expose_" v0 (nLct lct)
 
 -- -------------------------------------------------------------------------------------------------
 -- Jump, LCA
 -- -------------------------------------------------------------------------------------------------
 
--- | \(O(\log n)\) Returns the root of the underlying tree. Two vertices in the same connected
--- component have the same root vertex.
+-- | \(O(\log n)\) Returns the root of the underlying tree.
 --
+-- ==== Constraints
+-- - \(0 \le v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE root #-}
-root :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Int -> m Vertex
+root :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> m Vertex
 root lct c0 = stToPrim $ rootST lct c0
 
+-- | \(O(\log n)\) Returns whether the vertices \(u\) and \(v\) are in the same connected component
+-- (have the same `root`).
+--
+-- ==== Constraints
+-- - \(0 \le u, v \lt n\)
+--
+-- @since 1.5.1.0
+{-# INLINE same #-}
+same :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> Vertex -> m Bool
+same lct u v = stToPrim $ do
+  r1 <- rootST lct u
+  r2 <- rootST lct v
+  pure $ r1 == r2
+
 -- | \(O(\log n)\) Returns the parent vertex in the underlying tree.
 --
+-- ==== Constraints
+-- - \(0 \le v \lt n\)
+--
 -- @since 1.1.1.0
 {-# INLINE parent #-}
-parent :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Int -> m (Maybe Vertex)
-parent lct x = stToPrim $ parentST lct x
+parent :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> m (Maybe Vertex)
+parent lct v = stToPrim $ parentST lct v
 
--- | \(O(\log n)\) Given a path between \(u\) and \(v\), returns the \(k\)-th vertex of the path.
+-- | \(O(\log n)\) Given a path between \(u\) and \(v\), returns the \(k\)-th vertex from \(u\) in
+-- the path.
 --
 -- ==== Constraints
--- - The \(k\)-th vertex must exist.
+-- - \(0 \le u, v \lt n\)
+-- - \(0 \le k \lt \mathrm{|path|}\)
+-- - \(u\) and \(v\) must be in the same connected component, otherwise the vehavior is undefined.
 --
 -- @since 1.1.1.0
 {-# INLINE jump #-}
 jump :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> Vertex -> Int -> m Vertex
-jump lct u v k = stToPrim $ jumpST lct u v k
+jump lct u v k = stToPrim $ do
+  fromMaybe (error "AtCoder.Extra.Tree.Lct.jump: invalid jump") <$> jumpMaybeST lct u v k
 
+-- | \(O(\log n)\) Given a path between \(u\) and \(v\), returns the \(k\)-th vertex from \(u\) in
+-- the path.
+--
+-- ==== Constraints
+-- - \(0 \le u, v \lt n\)
+-- - \(u\) and \(v\) must be in the same connected component, otherwise the vehavior is undefined.
+--
+-- @since 1.5.1.0
+{-# INLINE jumpMaybe #-}
+jumpMaybe :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> Vertex -> Int -> m (Maybe Vertex)
+jumpMaybe lct u v k = stToPrim $ jumpMaybeST lct u v k
+
+-- | \(O(\log n)\) Returns the length of path between \(u\) and \(v\).
+--
+-- ==== Constraints
+-- - \(0 \le u, v \lt n\)
+-- - \(u\) and \(v\) must be in the same connected component.
+--
+-- @since 1.5.1.0
+{-# INLINE lengthBetween #-}
+lengthBetween :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> Vertex -> m Vertex
+lengthBetween lct u v = stToPrim $ do
+  lengthBetweenST lct u v
+
 -- | \(O(\log n)\) Returns the LCA of \(u\) and \(v\). Because the root of the underlying tree changes
 -- in almost every operation, one might want to use `evert` beforehand.
 --
 -- ==== Constraints
+-- - \(0 \le u, v \lt n\)
 -- - \(u\) and \(v\) must be in the same connected component.
 --
 -- @since 1.1.1.0
 {-# INLINE lca #-}
-lca :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Int -> Int -> m Vertex
+lca :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> Vertex -> m Vertex
 lca lct u v = stToPrim $ do
+  fromMaybe (error ("AtCoder.Extra.Tree.Lct.lca: given two vertices in different connected components " ++ show (u, v)))
+    <$> lcaMaybe lct u v
+
+-- | \(O(\log n)\) Returns the LCA of \(u\) and \(v\). Because the root of the underlying tree changes
+-- in almost every operation, one might want to use `evert` beforehand.
+--
+-- ==== Constraints
+-- - \(0 \le u, v \lt n\)
+--
+-- @since 1.5.1.0
+{-# INLINE lcaMaybe #-}
+lcaMaybe :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> Vertex -> m (Maybe Vertex)
+lcaMaybe lct u v = stToPrim $ do
   ru <- rootST lct u
   rv <- rootST lct v
-  let !_ = ACIA.runtimeAssert (ru == rv) $ "AtCoder.Extra.Lct.lca: given two vertices in different connected components " ++ show (u, v)
-  _ <- exposeST lct u
-  exposeST lct v
+  if ru == rv
+    then do
+      _ <- exposeST lct u
+      Just <$> exposeST lct v
+    else pure Nothing
 
 -- -------------------------------------------------------------------------------------------------
 -- Monoid product
@@ -382,6 +521,11 @@
 
 -- | Amortized \(O(\log n)\). Folds a path between \(u\) and \(v\) (inclusive).
 --
+-- ==== Constraints
+-- - \(0 \le u, v \lt n\)
+-- - \(u\) and \(v\) must be in the same connected component, otherwise the return value is
+--   nonsense.
+--
 -- @since 1.1.1.0
 {-# INLINE prodPath #-}
 prodPath :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => Lct (PrimState m) a -> Vertex -> Vertex -> m a
@@ -393,14 +537,15 @@
   -- now that @v@ is at the root of the auxiliary tree, its aggregation value is the path folding:
   VGM.unsafeRead prodLct v
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.prodPath" u (nLct lct)
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.prodPath" v (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.prodPath" u (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.prodPath" v (nLct lct)
 
 -- | Amortized \(O(\log n)\). Fold the subtree under \(v\), considering \(p\) as the root-side
 -- vertex. Or, if \(p\) equals \(v\), \(v\) will be the new root.
 --
 -- ==== Constraints
--- - The inverse operator has to be set on construction (`newInv` or `buildInv`).
+-- - The inverse operator must be set on construction (`newInv` or `buildInv`).
+-- - \(0 \le u, v \lt n\)
 --
 -- @since 1.1.1.0
 {-# INLINE prodSubtree #-}
@@ -408,14 +553,32 @@
   (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) =>
   -- | Link/cut tree
   Lct (PrimState m) a ->
-  -- | Vertex
+  -- | Vertex \(v\)
   Vertex ->
-  -- | Root or parent
+  -- | Parent \(p\) (not need be adjacent to \(v\)), or same as \(v\), making it a new root.
   Vertex ->
   -- | Subtree's monoid product
   m a
 prodSubtree lct v rootOrParent = stToPrim $ prodSubtreeST lct v rootOrParent
 
+-- | Amortized \(O(\log n)\). Fold a tree that contains \(v\).
+--
+-- ==== Constraints
+-- - The inverse operator must be set on construction (`newInv` or `buildInv`).
+-- - \(0 \le v \lt n\)
+--
+-- @since 1.5.1.0
+{-# INLINE prodTree #-}
+prodTree ::
+  (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) =>
+  -- | Link/cut tree
+  Lct (PrimState m) a ->
+  -- | Vertex \(v\)
+  Vertex ->
+  -- | Subtree's monoid product
+  m a
+prodTree lct v = stToPrim $ prodSubtreeST lct v v
+
 -- -------------------------------------------------------------------------------------------------
 -- Internal
 -- -------------------------------------------------------------------------------------------------
@@ -590,7 +753,7 @@
 -- Node operations
 -- -------------------------------------------------------------------------------------------------
 
--- | \(O(1)\) Propgates the lazily propagated values on a node.
+-- | \(O(1)\) Propagates the lazily propagated values on a node.
 {-# INLINEABLE pushNodeST #-}
 pushNodeST :: (VU.Unbox a) => Lct s a -> Vertex -> ST s ()
 pushNodeST lct@Lct {lLct, rLct, revLct} v = do
@@ -614,7 +777,7 @@
 {-# INLINEABLE swapLrNodeST #-}
 swapLrNodeST :: (VU.Unbox a) => Lct s a -> Vertex -> ST s ()
 swapLrNodeST Lct {lLct, rLct, prodLct, dualProdLct} i = do
-  -- swap chidlren
+  -- swap children
   VGM.unsafeModifyM lLct (VGM.unsafeExchange rLct i) i
   -- swap prodLct[i] and dualProdLct[i]
   VGM.unsafeModifyM prodLct (VGM.unsafeExchange dualProdLct i) i
@@ -686,6 +849,7 @@
   _ <- exposeST lct p
   pushNodeST lct p
 
+  -- no assertion and very unsafe..
   -- dbgM $ do
   --   cp <- VGM.unsafeRead pLct c
   --   let !_ = ACIA.runtimeAssert (nullLct cp) $ "cp must be null: " ++ show (c, cp)
@@ -698,6 +862,7 @@
   VGM.unsafeWrite rLct p c
   updateNodeST lct p
 
+-- no assertion and very unsafe..
 {-# INLINEABLE cutST #-}
 cutST :: (Monoid a, VU.Unbox a) => Lct s a -> Vertex -> Vertex -> ST s ()
 cutST lct@Lct {pLct, lLct} u v = do
@@ -773,7 +938,7 @@
   -- do
   --   -- FIXME: remove
   --   pRes <- VGM.unsafeRead pLct res
-  --   unless (nullLct pRes) $ error $ "xxx must be null!!! " ++ show (res, pRes)
+  --   unless (nullLct pRes) $ error $ "must be null!!! " ++ show (res, pRes)
 
   splayST lct v0
 
@@ -802,7 +967,7 @@
   splayST lct c'
   pure c'
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.rootST" c0 (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.rootST" c0 (nLct lct)
 
 {-# INLINEABLE parentST #-}
 parentST :: (HasCallStack, Monoid a, VU.Unbox a) => Lct s a -> Int -> ST s (Maybe Vertex)
@@ -823,36 +988,35 @@
                 inner yr
       Just <$> inner xl
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.parentST" x (nLct lct)
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.parentST" x (nLct lct)
 
-{-# INLINEABLE jumpST #-}
-jumpST :: (HasCallStack, Monoid a, VU.Unbox a) => Lct s a -> Vertex -> Vertex -> Int -> ST s Vertex
-jumpST lct@Lct {lLct, rLct, sLct} u0 v0 k0 = do
+{-# INLINEABLE jumpMaybeST #-}
+jumpMaybeST :: (HasCallStack, Monoid a, VU.Unbox a) => Lct s a -> Vertex -> Vertex -> Int -> ST s (Maybe Vertex)
+jumpMaybeST lct@Lct {lLct, rLct, sLct} u0 v0 k0 = do
   -- make @v0@ a new root of the underlying tree
   evertST lct v0
   -- make @u0@ in the same preferred path as the root (@v0)
   _ <- exposeST lct u0
-
-  do
-    size <- VGM.unsafeRead sLct u0
-    let !_ = ACIA.runtimeAssert (0 <= k0 && k0 < size) "invalid jump"
-    pure ()
-
-  let inner k u = do
-        pushNodeST lct u
-        -- TODO: what is happening?
-        ur <- VGM.unsafeRead rLct u
-        urSize <- if nullLct ur then pure 0 else VGM.unsafeRead sLct ur
-        case compare k urSize of
-          LT -> inner k ur
-          EQ -> pure u
-          GT -> do
-            ul <- VGM.unsafeRead lLct u
-            inner (k - (urSize + 1)) ul
+  size <- VGM.unsafeRead sLct u0
+  if 0 <= k0 && k0 < size
+    then do
+      let inner k u = do
+            pushNodeST lct u
+            -- TODO: what is happening?
+            ur <- VGM.unsafeRead rLct u
+            urSize <- if nullLct ur then pure 0 else VGM.unsafeRead sLct ur
+            case compare k urSize of
+              LT -> inner k ur
+              EQ -> pure u
+              GT -> do
+                ul <- VGM.unsafeRead lLct u
+                inner (k - (urSize + 1)) ul
 
-  res <- inner k0 u0
-  splayST lct res
-  pure res
+      res <- inner k0 u0
+      splayST lct res
+      pure $ Just res
+    else do
+      pure Nothing
 
 {-# INLINEABLE prodSubtreeST #-}
 prodSubtreeST ::
@@ -861,7 +1025,7 @@
   Lct s a ->
   -- | Vertex
   Vertex ->
-  -- | Root or parent
+  -- | Parent \(p\) (not need be adjacent to \(v\)), or same as \(v\), making it a new root.
   Vertex ->
   -- | Subtree's monoid product
   ST s a
@@ -872,8 +1036,8 @@
       evertST lct v
       VGM.unsafeRead subtreeProdLct v
     else do
-      -- @rootOrParent@ can be far. retrieve the adjacent vertex:
-      parent_ <- jumpST lct v rootOrParent 1
+      -- retrieve the adjacent parent:
+      parent_ <- fromJust <$> jumpMaybeST lct v rootOrParent 1
       -- detach @v@ from the parent. now that it's the root of the subtree vertices, the aggregation
       -- value is the aggregation of all the subtree vertices.
       cutST lct v parent_
@@ -882,5 +1046,15 @@
       linkST lct v parent_
       pure res
   where
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.prodSubtree" v nLct
-    !_ = ACIA.checkIndex "AtCoder.Extra.Lct.prodSubtree" rootOrParent nLct
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.prodSubtreeST" v nLct
+    !_ = ACIA.checkIndex "AtCoder.Extra.Tree.Lct.prodSubtreeST" rootOrParent nLct
+
+{-# INLINEABLE lengthBetweenST #-}
+lengthBetweenST :: (HasCallStack, Monoid a, VU.Unbox a) => Lct s a -> Vertex -> Vertex -> ST s Int
+lengthBetweenST lct@Lct {sLct} u0 v0 = do
+  -- make @v0@ a new root of the underlying tree
+  evertST lct v0
+  -- make @u0@ in the same preferred path as the root (@v0)
+  _ <- exposeST lct u0
+  size <- VGM.unsafeRead sLct u0
+  pure $ size - 1
diff --git a/src/AtCoder/Extra/Vector.hs b/src/AtCoder/Extra/Vector.hs
--- a/src/AtCoder/Extra/Vector.hs
+++ b/src/AtCoder/Extra/Vector.hs
@@ -1,35 +1,410 @@
--- | Miscellaneous vector methods.
+-- | Miscellaneous vector functions. These functions are __not__ the fastest implementations, but
+-- fills in some lacking features.
 --
 -- @since 1.2.2.0
 module AtCoder.Extra.Vector
-  ( argsort,
+  ( -- * Sort functions
+    argsort,
+
+    -- * Index compression
+    compress,
+
+    -- * Vector Utilities
+    iconcatMap,
+    concatMapM,
+    iconcatMapM,
+    mapAccumL,
+    chunks,
+
+    -- ** Monadic scanl
+    prescanlM,
+    prescanlM',
+    postscanlM,
+    postscanlM',
+    scanlM,
+    scanlM',
+    scanl1M,
+    scanl1M',
+
+    -- * Queries
+    maxRangeSum,
+    minRangeSum,
+    slideMinIndices,
+    slideMaxIndices,
   )
 where
 
+-- TODO: maybe add lexicographic permutations, combinations, and subsequences.
+
+import AtCoder.Internal.Assert qualified as ACIA
+import AtCoder.Internal.Queue qualified as Q
+import AtCoder.Extra.Bisect (lowerBound)
+import Control.Monad (when)
+import Control.Monad.Fix (fix)
+import Control.Monad.Primitive (PrimMonad)
+import Control.Monad.ST (ST, runST)
+import Control.Monad.Trans.State.Strict (StateT, runStateT, state)
+import Data.Ord (Down (..))
+import Data.Vector qualified as V
 import Data.Vector.Algorithms.Intro qualified as VAI
+import Data.Vector.Fusion.Bundle qualified as Bundle
+import Data.Vector.Fusion.Bundle.Monadic qualified as BundleM
+import Data.Vector.Fusion.Stream.Monadic qualified as S
 import Data.Vector.Generic qualified as VG
+import Data.Vector.Generic.Mutable qualified as VGM
 import Data.Vector.Unboxed qualified as VU
+import GHC.Stack (HasCallStack)
 
--- | \(O(n \log n)\) Returns indices of the vector, stably sorted by their value.
+-- | \(O(n \log n)\) Returns indices of the vector elements, stably sorted by their value.
 --
 -- ==== Example
--- >>> import Data.Vector.Algorithms.Intro qualified as VAI
+-- >>> import AtCoder.Extra.Vector qualified as EV
 -- >>> import Data.Vector.Unboxed qualified as VU
--- >>> argsort $ VU.fromList [0, 1, 0, 1, 0]
+-- >>> EV.argsort $ VU.fromList @Int [0, 1, 0, 1, 0]
 -- [0,2,4,1,3]
 --
 -- @since 1.2.3.0
 {-# INLINEABLE argsort #-}
-argsort :: (Ord a, VU.Unbox a) => VU.Vector a -> VU.Vector Int
+-- TODO: use generic vector
+argsort :: (HasCallStack, Ord a, VU.Unbox a) => VU.Vector a -> VU.Vector Int
 argsort xs =
   VU.modify
-    ( VAI.sortBy
-        ( \i j ->
-            ( compare (xs VG.! i) (xs VG.! j) <> compare i j
-            )
-        )
-    )
+    (VAI.sortBy (\i j -> compare (xs VG.! i) (xs VG.! j) <> compare i j))
     $ VU.generate (VU.length xs) id
 
--- TODO: maybe add lexicographic permutations, combinations, and subsequences.
+-- | \(O(n \log n)\) One dimensional index compression: xs -> (nubSortXs, xs')
+--
+-- ==== Example
+-- >>> import AtCoder.Extra.Bisect (lowerBound)
+-- >>> import AtCoder.Extra.Vector qualified as EV
+-- >>> import Data.Vector.Unboxed qualified as VU
+-- >>> let xs = VU.fromList [0, 20, 40, 10, 30]
+-- >>> let (dict, xs') = EV.compress xs
+-- >>> xs'
+-- [0,2,4,1,3]
+-- >>> lowerBound dict 10
+-- 1
+--
+-- @since 1.5.1.0
+{-# INLINE compress #-}
+compress :: VU.Vector Int -> (VU.Vector Int, VU.Vector Int)
+compress xs = (dict, VG.map (lowerBound dict) xs)
+  where
+    !dict = VG.uniq $ VG.modify VAI.sort xs
 
+-- | Maps each element to a vector and concatenate the results.
+--
+-- ==== Example
+-- >>> import AtCoder.Extra.Vector qualified as EV
+-- >>> import Data.Vector.Unboxed qualified as VU
+-- >>> EV.iconcatMap (\i x -> VU.fromList [i + x, i + x]) $ VU.replicate @Int 3 0
+-- [0,0,1,1,2,2]
+--
+-- @since 1.5.1.0
+{-# INLINE iconcatMap #-}
+iconcatMap :: (HasCallStack, VG.Vector v a, VG.Vector v b) => (Int -> a -> v b) -> v a -> v b
+iconcatMap f =
+  VG.unstream
+    . Bundle.concatVectors
+    . Bundle.inplace (S.map (uncurry f) . S.indexed) id
+    . VG.stream
+
+-- | Maps each element to a vector and concatenate the results.
+--
+-- ==== Example
+-- >>> import AtCoder.Extra.Vector qualified as EV
+-- >>> import Data.Vector.Unboxed qualified as VU
+-- >>> EV.iconcatMap (\x -> pure (VU.fromList [x, x])) $ VU.generate @Int 3 id
+-- [0,0,1,1,2,2]
+--
+-- @since 1.5.1.0
+{-# INLINE concatMapM #-}
+concatMapM ::
+  (HasCallStack, Monad m, VG.Vector v a, VG.Vector v b) =>
+  (a -> m (v b)) ->
+  v a ->
+  m (v b)
+concatMapM f =
+  VG.unstreamM
+    . BundleM.concatVectors
+    . BundleM.mapM f
+    . BundleM.fromVector
+
+-- | Maps each element to a vector and concatenate the results.
+--
+-- ==== Example
+-- >>> import AtCoder.Extra.Vector qualified as EV
+-- >>> import Data.Vector.Unboxed qualified as VU
+-- >>> EV.iconcatMapM (\i x -> pure (VU.fromList [i + x, i + x])) $ VU.replicate @Int 3 0
+-- [0,0,1,1,2,2]
+--
+-- @since 1.5.1.0
+{-# INLINE iconcatMapM #-}
+iconcatMapM ::
+  (HasCallStack, Monad m, VG.Vector v a, VG.Vector v b) =>
+  (Int -> a -> m (v b)) ->
+  v a ->
+  m (v b)
+iconcatMapM f =
+  VG.unstreamM
+    . BundleM.concatVectors
+    . BundleM.mapM (uncurry f)
+    . BundleM.indexed
+    . BundleM.fromVector
+
+-- | \(O(n)\) Maps a vector with an accumulator.
+--
+-- ==== Example
+-- >>> import AtCoder.Extra.Vector qualified as EV
+-- >>> import Data.Vector.Unboxed qualified as VU
+-- >>> EV.mapAccumL (\s x -> (s + 1, s * x)) (0 :: Int) $ VU.generate @Int 4 id
+-- (4,[0,1,4,9])
+--
+-- @since 1.5.1.0
+{-# INLINE mapAccumL #-}
+mapAccumL ::
+  forall v s a b.
+  (HasCallStack, VG.Vector v a, VG.Vector v b) =>
+  (s -> a -> (s, b)) ->
+  s ->
+  v a ->
+  (s, v b)
+mapAccumL f s0 xs = (\(!x, !s) -> (s, x)) $ runST $ (`runStateT` s0) $ do
+  unstreamPrimM
+    . BundleM.mapM g
+    $ BundleM.fromVector xs
+  where
+    g :: forall st. a -> StateT s (ST st) b
+    g a =
+      state
+        ( \s ->
+            let (!s', !b) = f s a
+             in (b, s')
+        )
+
+-- | https://github.com/haskell/vector/issues/416
+{-# INLINE [1] unstreamPrimM #-}
+unstreamPrimM :: (PrimMonad m, VG.Vector v a) => BundleM.Bundle m u a -> m (v a)
+unstreamPrimM s = VGM.munstream s >>= VG.unsafeFreeze
+
+-- | @since 1.5.1.0
+{-# INLINE prescanlM #-}
+prescanlM :: (HasCallStack, Monad m, VG.Vector v a, VG.Vector v b) => (a -> b -> m a) -> a -> v b -> m (v a)
+prescanlM f x0 =
+  VG.unstreamM
+    . prescanlMB f x0
+    . VG.stream
+
+-- | @since 1.5.1.0
+{-# INLINE prescanlM' #-}
+prescanlM' :: (HasCallStack, Monad m, VG.Vector v a, VG.Vector v b) => (a -> b -> m a) -> a -> v b -> m (v a)
+prescanlM' f x0 =
+  VG.unstreamM
+    . prescanlMB' f x0
+    . VG.stream
+
+-- | @since 1.5.1.0
+{-# INLINE postscanlM #-}
+postscanlM :: (HasCallStack, Monad m, VG.Vector v a, VG.Vector v b) => (a -> b -> m a) -> a -> v b -> m (v a)
+postscanlM f x0 =
+  VG.unstreamM
+    . postscanlMB f x0
+    . VG.stream
+
+-- | @since 1.5.1.0
+{-# INLINE postscanlM' #-}
+postscanlM' :: (HasCallStack, Monad m, VG.Vector v a, VG.Vector v b) => (a -> b -> m a) -> a -> v b -> m (v a)
+postscanlM' f x0 =
+  VG.unstreamM
+    . postscanlMB' f x0
+    . VG.stream
+
+-- | @since 1.5.1.0
+{-# INLINE scanlM #-}
+scanlM :: (HasCallStack, Monad m, VG.Vector v a, VG.Vector v b) => (a -> b -> m a) -> a -> v b -> m (v a)
+scanlM f x0 =
+  VG.unstreamM
+    . scanlMB f x0
+    . VG.stream
+
+-- | @since 1.5.1.0
+{-# INLINE scanlM' #-}
+scanlM' :: (HasCallStack, Monad m, VG.Vector v a, VG.Vector v b) => (a -> b -> m a) -> a -> v b -> m (v a)
+scanlM' f x0 =
+  VG.unstreamM
+    . scanlMB' f x0
+    . VG.stream
+
+-- | @since 1.5.1.0
+{-# INLINE scanl1M #-}
+scanl1M :: (HasCallStack, Monad m, VG.Vector v a) => (a -> a -> m a) -> v a -> m (v a)
+scanl1M f =
+  VG.unstreamM
+    . scanl1MB f
+    . VG.stream
+
+-- | @since 1.5.1.0
+{-# INLINE scanl1M' #-}
+scanl1M' :: (HasCallStack, Monad m, VG.Vector v a) => (a -> a -> m a) -> v a -> m (v a)
+scanl1M' f =
+  VG.unstreamM
+    . scanl1MB' f
+    . VG.stream
+
+{-# INLINE prescanlMB #-}
+prescanlMB :: (Monad m) => (a -> b -> m a) -> a -> Bundle.Bundle v b -> BundleM.Bundle m v a
+prescanlMB f x0 = BundleM.prescanlM f x0 . Bundle.lift
+
+{-# INLINE prescanlMB' #-}
+prescanlMB' :: (Monad m) => (a -> b -> m a) -> a -> Bundle.Bundle v b -> BundleM.Bundle m v a
+prescanlMB' f x0 = BundleM.prescanlM' f x0 . Bundle.lift
+
+{-# INLINE postscanlMB #-}
+postscanlMB :: (Monad m) => (a -> b -> m a) -> a -> Bundle.Bundle v b -> BundleM.Bundle m v a
+postscanlMB f x0 = BundleM.postscanlM f x0 . Bundle.lift
+
+{-# INLINE postscanlMB' #-}
+postscanlMB' :: (Monad m) => (a -> b -> m a) -> a -> Bundle.Bundle v b -> BundleM.Bundle m v a
+postscanlMB' f x0 = BundleM.postscanlM' f x0 . Bundle.lift
+
+{-# INLINE scanlMB #-}
+scanlMB :: (Monad m) => (a -> b -> m a) -> a -> Bundle.Bundle v b -> BundleM.Bundle m v a
+scanlMB f x0 = BundleM.scanlM f x0 . Bundle.lift
+
+{-# INLINE scanlMB' #-}
+scanlMB' :: (Monad m) => (a -> b -> m a) -> a -> Bundle.Bundle v b -> BundleM.Bundle m v a
+scanlMB' f x0 = BundleM.scanlM' f x0 . Bundle.lift
+
+{-# INLINE scanl1MB #-}
+scanl1MB :: (Monad m) => (a -> a -> m a) -> Bundle.Bundle v a -> BundleM.Bundle m v a
+scanl1MB f = BundleM.scanl1M f . Bundle.lift
+
+{-# INLINE scanl1MB' #-}
+scanl1MB' :: (Monad m) => (a -> a -> m a) -> Bundle.Bundle v a -> BundleM.Bundle m v a
+scanl1MB' f = BundleM.scanl1M' f . Bundle.lift
+
+-- | \(O(n)\) Converts a vector into chunks of vectors with lenth \(k\). The last vector may have
+-- smaller length than \(k\).
+--
+-- >>> import AtCoder.Extra.Vector qualified as EV
+-- >>> import Data.Vector.Unboxed qualified as VU
+-- >>> EV.chunks 3 $ VU.fromList ([1, 2, 3, 4, 5, 6, 7] :: [Int])
+-- [[1,2,3],[4,5,6],[7]]
+--
+-- @since 1.5.1.0
+{-# INLINE chunks #-}
+chunks :: (VG.Vector v a) => Int -> v a -> V.Vector (v a)
+chunks len xs0 = V.unfoldrExactN n step xs0
+  where
+    n = (VG.length xs0 + len - 1) `div` len
+    step xs = (VG.take len xs, VG.drop len xs)
+
+-- | \(O(n)\) Returns maximum range sum.
+--
+-- ==== Example
+-- >>> import AtCoder.Extra.Vector qualified as EV
+-- >>> import Data.Vector.Unboxed qualified as VU
+-- >>> EV.maxRangeSum $ VU.fromList @Int [-3, 1, 6, -2, 7, -5]
+-- 12
+--
+-- @since 1.5.1.0
+{-# INLINE maxRangeSum #-}
+maxRangeSum :: forall v a. (VG.Vector v a, Ord a, Num a) => v a -> a
+maxRangeSum xs = fst $ VG.foldl' f (0 :: a, 0 :: a) csum
+  where
+    csum = VG.postscanl' (+) (0 :: a) xs
+    f (!acc, !minL) x = (max acc (x - minL), min minL x)
+
+-- | \(O(n)\) Returns minimum range sum.
+--
+-- ==== Example
+-- >>> import AtCoder.Extra.Vector qualified as EV
+-- >>> import Data.Vector.Unboxed qualified as VU
+-- >>> EV.minRangeSum $ VU.fromList @Int[-3, 1, 6, -20, 7, -9]
+-- -22
+--
+-- @since 1.5.1.0
+{-# INLINE minRangeSum #-}
+minRangeSum :: forall v a. (VG.Vector v a, Ord a, Num a) => v a -> a
+minRangeSum xs = fst $ VG.foldl' f (0 :: a, 0 :: a) csum
+  where
+    csum = VG.postscanl' (+) (0 :: a) xs
+    f (!acc, !maxL) x = (min acc (x - maxL), max maxL x)
+
+-- | \(O(N)\) Returns indices of minimum values in the windows with the specified length.
+--
+-- ==== Constraints
+-- - \(k \gt 0\)
+--
+-- ==== Example
+--
+-- >>> slideMinIndices 3 (VU.fromList [0 .. 5])
+-- [0,1,2,3]
+-- >>> slideMinIndices 3 (VU.fromList [5, 4 .. 0])
+-- [2,3,4,5]
+--
+-- @since 1.5.1.0
+{-# INLINE slideMinIndices #-}
+slideMinIndices :: (HasCallStack) => Int -> VU.Vector Int -> VU.Vector Int
+slideMinIndices k xs
+  | VU.null xs = VU.empty
+  | k >= VU.length xs = VU.singleton $ VU.minIndex xs
+  | otherwise = slideCmpIndicesOn Down k xs
+  where
+    !_ = ACIA.runtimeAssert (k > 0) "AtCoder.Extra.Vector.slideMinIndices: given non-positive k"
+
+-- | \(O(N)\) Returns indices of maximum values in the windows with the specified length.
+--
+-- ==== Constraints
+-- - \(k \gt 0\)
+--
+-- ==== Example
+--
+-- @
+-- indices: 0 1 2 3 4 5
+-- values:  0 1 2 3 4 5   max value indices:
+--          [---]         2
+--            [---]       3
+--              [---]     4
+--                [---]   5
+-- @
+--
+-- >>> slideMaxIndices 3 (VU.fromList [0 .. 5])
+-- [2,3,4,5]
+-- >>> slideMaxIndices 3 (VU.fromList [5, 4 .. 0])
+-- [0,1,2,3]
+--
+-- @since 1.5.1.0
+{-# INLINE slideMaxIndices #-}
+slideMaxIndices :: (HasCallStack) => Int -> VU.Vector Int -> VU.Vector Int
+slideMaxIndices k xs
+  | VU.null xs = VU.empty
+  | k >= VU.length xs = VU.singleton $ VU.maxIndex xs
+  | otherwise = slideCmpIndicesOn id k xs
+  where
+    !_ = ACIA.runtimeAssert (k > 0) "AtCoder.Extra.Vector.slideMaxIndices: given non-positive k"
+
+-- | \(O(N)\) (1) in <https://qiita.com/kuuso1/items/318d42cd089a49eeb332>
+{-# INLINE slideCmpIndicesOn #-}
+slideCmpIndicesOn :: (VG.Vector v a, Ord b) => (a -> b) -> Int -> v a -> VU.Vector Int
+slideCmpIndicesOn wrap len xs = runST $ do
+  -- dequeue of maximum number indices.
+  !buf <- Q.new (VG.length xs)
+
+  fmap (VU.drop (len - 1)) $ VG.generateM (VG.length xs) $ \i -> do
+    -- remove the front indices that are no longer in the span
+    fix $ \loop -> do
+      b <- maybe False (<= i - len) <$> Q.readMaybeFront buf 0
+      when b $ do
+        Q.popFront_ buf
+        loop
+
+    -- remove the last indices that are less attractive to the new coming value
+    fix $ \loop -> do
+      b <- maybe False ((< wrap (xs VG.! i)) . wrap . (xs VG.!)) <$> Q.readMaybeBack buf 0
+      when b $ do
+        Q.popBack_ buf
+        loop
+
+    Q.pushBack buf i
+    Q.readFront buf 0
diff --git a/src/AtCoder/Internal/Csr.hs b/src/AtCoder/Internal/Csr.hs
--- a/src/AtCoder/Internal/Csr.hs
+++ b/src/AtCoder/Internal/Csr.hs
@@ -43,6 +43,7 @@
     -- * Accessors
     adj,
     adjW,
+    adj1,
     eAdj,
     eAdjW,
   )
@@ -157,6 +158,16 @@
   let il = startCsr VG.! i
       ir = startCsr VG.! (i + 1)
    in VU.zip (VU.slice il (ir - il) adjCsr) (VU.slice il (ir - il) wCsr)
+
+-- | \(O(1)\) Returns the adjacent vertices with \(1\) weights.
+--
+-- @since 1.5.1.0
+{-# INLINE adj1 #-}
+adj1 :: (HasCallStack, VU.Unbox w) => Csr w -> Int -> VU.Vector (Int, Int)
+adj1 Csr {..} i =
+  let il = startCsr VG.! i
+      ir = startCsr VG.! (i + 1)
+   in VU.zip (VU.slice il (ir - il) adjCsr) (VU.replicate (ir - il) 1)
 
 -- | \(O(n)\) Returns a vector of @(csrEdgeIndex, adjacentVertex)@.
 --
diff --git a/test/Main.hs b/test/Main.hs
--- a/test/Main.hs
+++ b/test/Main.hs
@@ -32,6 +32,9 @@
 import Tests.Extra.Semigroup.Permutation qualified
 import Tests.Extra.Seq qualified
 import Tests.Extra.Seq.Map qualified
+import Tests.Extra.Tree qualified
+import Tests.Extra.Tree.Lct qualified
+import Tests.Extra.Vector qualified
 import Tests.Extra.WaveletMatrix qualified
 import Tests.Extra.WaveletMatrix.BitVector qualified
 import Tests.Extra.WaveletMatrix.Raw qualified
@@ -91,6 +94,9 @@
             testGroup "Semigroup.Permutation" Tests.Extra.Semigroup.Permutation.tests,
             testGroup "Seq" Tests.Extra.Seq.tests,
             testGroup "Seq.Map" Tests.Extra.Seq.Map.tests,
+            testGroup "Tree" Tests.Extra.Tree.tests,
+            testGroup "Tree.Lct" Tests.Extra.Tree.Lct.tests,
+            testGroup "Vector" Tests.Extra.Vector.tests,
             testGroup "WaveletMatrix" Tests.Extra.WaveletMatrix.tests,
             testGroup "WaveletMatrix.BitVector" Tests.Extra.WaveletMatrix.BitVector.tests,
             testGroup "WaveletMatrix.Raw" Tests.Extra.WaveletMatrix.Raw.tests,
diff --git a/test/Tests/Extra/KdTree.hs b/test/Tests/Extra/KdTree.hs
--- a/test/Tests/Extra/KdTree.hs
+++ b/test/Tests/Extra/KdTree.hs
@@ -33,7 +33,7 @@
     -- n <- QC.chooseInt (1, 256)
     n <- QC.chooseInt (1, 8)
     -- q <- QC.chooseInt (1, 256)
-    q <- QC.chooseInt (1, 4)
+    q <- QC.chooseInt (1, 32)
     refVec <- (VU.fromList <$>) $ QC.vectorOf n $ do
       x <- QC.chooseInt rng
       y <- QC.chooseInt rng
diff --git a/test/Tests/Extra/LazyKdTree.hs b/test/Tests/Extra/LazyKdTree.hs
--- a/test/Tests/Extra/LazyKdTree.hs
+++ b/test/Tests/Extra/LazyKdTree.hs
@@ -44,7 +44,7 @@
 instance QC.Arbitrary Init where
   arbitrary = do
     n <- QC.chooseInt (1, 256)
-    q <- QC.chooseInt (1, 256)
+    q <- QC.chooseInt (1, 1024)
     xyws <- (VU.fromList <$>) $ QC.vectorOf n $ do
       x <- QC.chooseInt rng
       y <- QC.chooseInt rng
diff --git a/test/Tests/Extra/Tree.hs b/test/Tests/Extra/Tree.hs
new file mode 100644
--- /dev/null
+++ b/test/Tests/Extra/Tree.hs
@@ -0,0 +1,27 @@
+module Tests.Extra.Tree where
+
+import AtCoder.Extra.Graph qualified as Gr
+import AtCoder.Extra.Tree qualified as Tree
+import Data.Bit (Bit (..))
+import Data.Ord (Down (..), comparing)
+import Data.Semigroup (Sum (..))
+import Data.Vector.Unboxed qualified as VU
+import Test.Tasty
+import Test.Tasty.HUnit
+
+unit_one :: TestTree
+unit_one = testCase "one" $ do
+  let n = 1
+  let !t = Gr.build @Int n VU.empty
+  Tree.diameter n (Gr.adjW t) (-1) @?= ((0, 0), 0)
+  Tree.diameterPath 1 (Gr.adjW t) (-1) @?= (VU.singleton 0, 0)
+  Tree.mst n (VU.singleton (0, 0, 100)) @?= (0, VU.singleton (Bit False), t)
+  Tree.mstBy (comparing Down) 1 (VU.singleton (0, 0, 100)) @?= (0, VU.singleton (Bit False), t)
+  Tree.fold (Gr.adjW t) (const (Sum (10 :: Int))) (\(Sum a) (!_, !dw) -> Sum (a + dw)) (<>) 0 @?= 10
+  Tree.scan n (Gr.adjW t) (const (Sum (10 :: Int))) (\(Sum a) (!_, !dw) -> Sum (a + dw)) (<>) 0 @?= VU.singleton 10
+  Tree.foldReroot n (Gr.adjW t) (const (Sum (10 :: Int))) (\(Sum a) (!_, !dw) -> Sum (a + dw)) (<>) @?= VU.singleton 10
+
+tests :: [TestTree]
+tests =
+  [ unit_one
+  ]
diff --git a/test/Tests/Extra/Tree/Lct.hs b/test/Tests/Extra/Tree/Lct.hs
new file mode 100644
--- /dev/null
+++ b/test/Tests/Extra/Tree/Lct.hs
@@ -0,0 +1,315 @@
+{-# LANGUAGE RecordWildCards #-}
+
+-- | Not complete!
+module Tests.Extra.Tree.Lct where
+
+import AtCoder.Extra.Graph qualified as Gr
+import AtCoder.Extra.Tree.Lct (Lct (..))
+import AtCoder.Extra.Tree.Lct qualified as Lct
+import AtCoder.Internal.Buffer qualified as B
+import Control.Monad (when)
+import Control.Monad.Primitive (PrimMonad, PrimState)
+import Data.Foldable (for_)
+import Data.Maybe (isJust, fromJust)
+import Data.Semigroup (Sum (..))
+import Data.Vector qualified as V
+import Data.Vector.Generic qualified as VG
+import Data.Vector.Generic.Mutable qualified as VGM
+import Data.Vector.Unboxed qualified as VU
+import Data.Vector.Unboxed.Mutable qualified as VUM
+import Test.QuickCheck.Monadic as QCM
+import Test.Tasty
+import Test.Tasty.QuickCheck as QC
+
+-- | Maximum number of vertices.
+maxN :: Int
+maxN = 16
+
+-- | Maximum number of queries.
+maxQ :: Int
+maxQ = 1000
+
+-- | Range of monoid value.
+rng :: (Int, Int)
+rng = (-rngI, rngI)
+
+-- | Range of monoid value.
+rngI :: Int
+rngI = 10
+
+data Init s = Init
+  { n :: {-# UNPACK #-} !Int,
+    q :: {-# UNPACK #-} !Int,
+    initialMonoids :: {-# UNPACK #-} !(VU.Vector (Sum Int)),
+    refM :: !(IO Ref),
+    lctM :: !(IO (Lct.Lct s (Sum Int)))
+  }
+
+-- | Referencial implementation
+data Ref = Ref
+  { nR :: {-# UNPACK #-} !Int,
+    -- | Monoid values
+    vecR :: !(VUM.MVector VUM.RealWorld (Sum Int)),
+    -- | NxN
+    edgeR :: !(VUM.MVector VUM.RealWorld Bool)
+  }
+
+newRef :: Int -> VU.Vector (Sum Int) -> IO Ref
+newRef nR xs = do
+  vecR <- VU.thaw xs
+  edgeR <- VUM.replicate (nR * nR) False
+  pure Ref {..}
+
+instance Show (Init s) where
+  show Init {..} = show (n, q, initialMonoids)
+
+instance QC.Arbitrary (Init VUM.RealWorld) where
+  arbitrary = do
+    n <- QC.chooseInt (1, maxN)
+    q <- QC.chooseInt (1, maxQ)
+    initialMonoids <- (VU.fromList <$>) $ QC.vectorOf n (Sum <$> QC.chooseInt rng)
+    let lctM = Lct.buildInv negate initialMonoids VU.empty
+    let refM = let nR = n in newRef n initialMonoids
+    pure Init {..}
+
+data Query
+  = Read Int
+  | Write Int (Sum Int)
+  | Modify Int (Sum Int)
+  | ModifyM Int (Sum Int)
+  | Link Int Int
+  | Cut Int Int
+  | Evert Int
+  | Expose Int
+  | Expose_ Int
+  | -- | Root Int
+    Same Int Int
+  | -- TODO: Test parent!
+
+    -- | Parent Int
+    -- | Jump
+    JumpMaybe Int Int Int
+  | LengthBetween Int Int
+  | -- | Lca Int Int
+    -- | LcaMaybe Int Int
+    ProdPath Int Int
+  | ProdSubtree Int Int
+  -- -- | ProdTree Int Int
+  deriving (Eq, Show)
+
+-- | Arbitrary return type for the `Query` result.
+data Result
+  = None
+  | I !Int
+  | MI !(Maybe Int)
+  | B Bool
+  | SI (Sum Int)
+  deriving (Show, Eq)
+
+queryGen :: Int -> QC.Gen Query
+queryGen n =
+  QC.oneof
+    [ Read <$> u,
+      Write <$> u <*> x,
+      Modify <$> u <*> x,
+      ModifyM <$> u <*> x,
+      Link <$> u <*> v,
+      Cut <$> u <*> v,
+      Evert <$> u,
+      Expose <$> u,
+      Expose_ <$> u,
+      -- \| Root <$> Int
+      Same <$> u <*> v,
+      -- \| Parent <$> Int
+      -- \| Jump
+      JumpMaybe <$> u <*> v <*> k,
+      LengthBetween <$> u <*> v,
+      -- Lca <$> u <*> v,
+      -- LcaMaybe <$> u <*> v,
+      ProdPath <$> u <*> v
+      -- TODO: test prod subtree
+      -- ProdSubtree <$> u <*> v
+    ]
+  where
+    u = QC.chooseInt (0, n - 1)
+    v = QC.chooseInt (0, n - 1)
+    x = Sum <$> QC.chooseInt rng
+    k = QC.chooseInt (-1, n)
+
+toCsrM :: Ref -> IO (Gr.Csr ())
+toCsrM Ref {..} = do
+  es <- VU.unsafeFreeze edgeR
+  let uvs = (`VU.mapMaybe` VU.generate (nR * nR) id) $ \i ->
+        let (!u, !v) = i `divMod` nR
+         in if u /= v && es VG.! (nR * u + v)
+              then Just (u, v)
+              else Nothing
+  pure $ Gr.build' nR uvs
+
+findPath :: Gr.Csr () -> Int -> Int -> Maybe (VU.Vector Int)
+findPath gr v1 v2
+  | ws VG.! v2 == -1 = Nothing
+  | otherwise = Just $ Gr.constructPathFromRoot tree v2
+  where
+    (!ws, !tree) = Gr.trackingBfs (Gr.nCsr gr) grF (-1 :: Int) src
+    grF = Gr.adj1 gr
+    src = VU.singleton (v1, 0)
+
+-- | containers. (referencial implementation)
+handleRef :: Ref -> Query -> IO Result
+handleRef ref@Ref {..} q = do
+  case q of
+    Read i -> do
+      SI <$> VGM.read vecR i
+    Write i x -> do
+      VGM.write vecR i x
+      pure None
+    Modify i x -> do
+      VGM.modify vecR (+ x) i
+      pure None
+    ModifyM i x -> do
+      VGM.modifyM vecR (pure . (+ x)) i
+      pure None
+    Link u v -> do
+      VGM.write edgeR (nR * u + v) True
+      VGM.write edgeR (nR * v + u) True
+      pure None
+    Cut u v -> do
+      VGM.write edgeR (nR * u + v) False
+      VGM.write edgeR (nR * v + u) False
+      pure None
+    Evert _ -> do
+      pure None
+    Expose _ -> do
+      pure None
+    Expose_ _ -> do
+      pure None
+    -- \| Root Int
+    Same u v -> do
+      gr <- toCsrM ref
+      pure . B . isJust $ findPath gr u v
+    -- \| Parent Int
+    -- \| Jump
+    JumpMaybe u v k -> do
+      gr <- toCsrM ref
+      case findPath gr u v of
+        Nothing -> pure $ MI Nothing
+        Just path -> pure . MI $ path VG.!? k
+    LengthBetween u v -> do
+      gr <- toCsrM ref
+      let !path = fromJust $ findPath gr u v
+      pure . I $ VU.length path - 1
+    -- LcaMaybe u v -> do
+    --   gr <- toCsrM ref
+    --   -- FIXME: HLD cannot be used
+    --   let connected = isJust $ findPath gr u v
+    --   if connected
+    --     then pure . MI . Just $ Hld.lca (Hld.new gr) u v
+    --     else pure $ MI Nothing
+    ProdPath u v -> do
+      gr <- toCsrM ref
+      ws <- VU.unsafeFreeze vecR
+      case findPath gr u v of
+        Nothing -> pure $ SI mempty
+        Just path -> pure . SI $ VU.foldMap (ws VG.!) path
+    ProdSubtree u v -> do
+      -- TODO: fix this
+      pure None
+
+-- | ACL
+handleAcl :: (PrimMonad m) => Lct (PrimState m) (Sum Int) -> Query -> m Result
+handleAcl lct q = case q of
+  Read i -> do
+    SI <$> Lct.read lct i
+  Write i x -> do
+    Lct.write lct i x
+    pure None
+  Modify i x -> do
+    Lct.modify lct (+ x) i
+    pure None
+  ModifyM i x -> do
+    Lct.modifyM lct (pure . (+ x)) i
+    pure None
+  Link u v -> do
+    -- TODO: detect if already connect
+    Lct.link lct u v
+    pure None
+  Cut u v -> do
+    Lct.cut lct u v
+    pure None
+  Evert i -> do
+    Lct.evert lct i
+    pure None
+  Expose i -> do
+    _ <- Lct.expose lct i
+    pure None
+  Expose_ i -> do
+    Lct.expose_ lct i
+    pure None
+  -- \| Root Int
+  Same u v -> do
+    B <$> Lct.same lct u v
+  -- \| Parent Int
+  -- \| Jump
+  JumpMaybe u v k -> do
+    MI <$> Lct.jumpMaybe lct u v k
+  LengthBetween u v -> do
+    I <$> Lct.lengthBetween lct u v
+  -- LcaMaybe u v -> do
+  --   MI <$> Lct.lcaMaybe lct u v
+  ProdPath u v -> do
+    SI <$> Lct.prodPath lct u v
+  ProdSubtree u v -> do
+    SI <$> Lct.prodSubtree lct u v
+
+prop_randomTest :: Init VUM.RealWorld -> QC.Property
+prop_randomTest Init {..} = QCM.monadicIO $ do
+  lct <- QCM.run lctM
+  ref <- QCM.run refM
+  qs <- QCM.pick $ QC.vectorOf q (queryGen n)
+  record <- QCM.run $ B.new q
+  for_ qs $ \query -> do
+    b <- QCM.run $ filterRefM ref query
+    when b $ do
+      QCM.run $ B.pushBack record $ VU.DoNotUnboxStrict query
+      record' <- QCM.run $ B.unsafeFreeze record
+      expected <- QCM.run $ handleRef ref query
+      result <- QCM.run $ handleAcl lct query
+      when (expected /= result) $ do
+        -- show all the run queries on failure
+        let s = show . V.map (\(VU.DoNotUnboxStrict x) -> x) $ VU.convert record'
+        QCM.monitor $ QC.counterexample $ "actual queries in use: " ++ s
+      QCM.assertWith (expected == result) $ show (query, expected, result)
+
+-- | Filter invalid queries
+filterRefM :: Ref -> Query -> IO Bool
+filterRefM ref@Ref {..} q = case q of
+  Read _ -> pure True
+  Write _ _ -> pure True
+  Modify _ _ -> pure True
+  ModifyM _ _ -> pure True
+  Link u v
+    | u == v -> pure False
+    | otherwise -> not <$> same u v
+  Cut u v -> do
+    es <- VU.unsafeFreeze edgeR
+    pure $ es VG.! (nR * u + v)
+  Evert _ -> pure True
+  Expose _ -> pure True
+  Expose_ _ -> pure True
+  Same _ _ -> pure True
+  JumpMaybe u v _ -> same u v
+  LengthBetween u v -> same u v
+  -- LcaMaybe u v -> same u v
+  ProdPath u v -> same u v
+  ProdSubtree u v -> same u v
+  where
+    same u v = do
+      gr <- toCsrM ref
+      let connected = isJust $ findPath gr u v
+      pure connected
+
+tests :: [TestTree]
+tests =
+  [ QC.testProperty "randomTest" prop_randomTest
+  ]