diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,32 +1,41 @@
 # Revision history for acl-hs
 
+## 1.3.0.0 -- April 2025
+
+- Added `AtCoder.Extra.Math.isPrimitiveRoot`.
+- Re-created `AtCoder.Extra.Bisect` module.
+- Removed re-exports of `AtCoder.Internal.Math` functions from `AtCoder.Extra.Math`.
+- Removed `Extra.Vector.unsafePermuteInPlace`.
+- Removed some `Handle` function re-exports from `AtCoder.Extra.Seq` module.
+- Changed `AtCoder.Extra.Tree.scan` to non-generic, `Unbox` vector.
+
 ## 1.2.6.0 -- April 2025
 
-- Added `AtCoder.Extra.Math` functions
+- Added `AtCoder.Extra.Math` functions:
   - `isPrime`
   - `primes`
   - `primeFactors`
-- Added `AtCoderExtra.Math.Montgomery64`
-- Added `AtCoderExtra.ModInt64`
+- Added `AtCoderExtra.Math.Montgomery64`.
+- Added `AtCoderExtra.ModInt64`.
 
 ## 1.2.5.0 -- April 2025
 
-- Added `AtCoder.Extra.Mo`
-- Added `AtCoder.Extra.SqrtDecomposition`
+- Added `AtCoder.Extra.Mo`.
+- Added `AtCoder.Extra.SqrtDecomposition`.
 
 ## 1.2.4.0 -- April 2025
 
-- Added `AtCoder.Dsu.mergeMaybe`
-- Added `AtCoder.Extra.Graph` functions
+- Added `AtCoder.Dsu.mergeMaybe`.
+- Added `AtCoder.Extra.Graph` functions:
   - `rev`
   - `connectedComponents`
   - `bipartiteVertexColors`
-  - BFS, Dijkstra, Bellman–ford, Floyd–Warshall
+  - BFS, Dijkstra, Bellman–Ford, Floyd–Warshall
   - path reconstruction functions
-- Added `AtCoder.Extra.Tree` functions
+- Added `AtCoder.Extra.Tree` functions:
   - `diameter`, `diameterPath`
   - `mst`, `mstBy`
-- Added `AtCoder.Internal.Queue.newDeque`
+- Added `AtCoder.Internal.Queue.newDeque`.
 
 ## 1.2.3.0 -- March 2025
 
diff --git a/ac-library-hs.cabal b/ac-library-hs.cabal
--- a/ac-library-hs.cabal
+++ b/ac-library-hs.cabal
@@ -4,7 +4,7 @@
 -- PVP summary:  +-+------- breaking API changes
 --               | | +----- non-breaking API additions
 --               | | | +--- code changes with no API change
-version:         1.2.6.0
+version:         1.3.0.0
 synopsis:        Data structures and algorithms
 description:
   Haskell port of [ac-library](https://github.com/atcoder/ac-library), a library for competitive
@@ -39,7 +39,7 @@
     , bitvec             <1.2
     , bytestring         <0.14
     , primitive          >=0.6.4.0 && <0.10
-    , random             >=1.2.0
+    , random             >=1.2.0 && < 1.3
     , vector             >=0.13.0  && <0.14
     , vector-algorithms  <0.10
     , wide-word          <0.2
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
@@ -13,15 +13,14 @@
 -- Perform index compression:
 --
 -- >>> import AtCoder.Extra.Bisect
--- >>> import Data.Maybe (fromJust)
 -- >>> import Data.Vector.Algorithms.Intro qualified as VAI
 -- >>> import Data.Vector.Unboxed qualified as VU
 -- >>> let xs = VU.fromList ([0, 20, 10, 40, 30] :: [Int])
 -- >>> let dict = VU.uniq $ VU.modify VAI.sort xs
--- >>> VU.map (fromJust . lowerBound dict) xs
+-- >>> VU.map (lowerBound dict) xs
 -- [0,2,1,4,3]
 --
--- @since 1.1.0.0
+-- @since 1.3.0.0
 module AtCoder.Extra.Bisect
   ( -- * C++-like binary search
     lowerBound,
@@ -30,149 +29,122 @@
     upperBoundIn,
 
     -- * Generic bisection method
-    bisectL,
-    bisectLM,
-    bisectR,
-    bisectRM,
+    maxRight,
+    maxRightM,
+    minLeft,
+    minLeftM,
   )
 where
 
 import AtCoder.Internal.Assert qualified as ACIA
-import Data.Functor ((<&>))
 import Data.Functor.Identity
 import Data.Vector.Generic qualified as VG
 import GHC.Stack (HasCallStack)
 
--- | \(O(\log n)\) Bisection method implementation. Works on a half-open interfal \([l, r)\) .
---
--- @since 1.1.0.0
-{-# INLINE bisectLImpl #-}
-bisectLImpl :: (HasCallStack, Monad m) => (Int -> m Bool) -> Int -> Int -> m Int
-bisectLImpl p l0 = inner (l0 - 1)
-  where
-    inner l r
-      | l + 1 == r = pure l
-      | otherwise =
-          p mid >>= \case
-            True -> inner mid r
-            False -> inner l mid
-      where
-        mid = (l + r) `div` 2
-
--- | \(O(\log n)\) Bisection method implementation. Works on a half-open interfal \([l, r)\) .
---
--- @since 1.1.0.0
-{-# INLINE bisectRImpl #-}
-bisectRImpl :: (HasCallStack, Monad m) => (Int -> m Bool) -> Int -> Int -> m Int
-bisectRImpl p l = ((+ 1) <$>) . bisectLImpl p l
-
--- | \(O(\log n)\) Returns the index of the first element \(x\) in the vector such that
--- \(x \ge x_0\), or `Nothing` if no such element exists.
+-- | \(O(\log n)\) Returns the maximum \(r\) where \(x \lt x_i\) holds for \(i \in [0, r)\).
 --
 -- @
 -- Y Y Y Y Y N N N N N      Y: (< x0)
 -- --------- *---------> X  N: (>= x0)
---           R              R: returning point
+--           R              R: the right boundary point returned
 -- @
 --
 -- ==== __Example__
 -- >>> import Data.Vector.Unboxed qualified as VU
 -- >>> let xs = VU.fromList [1, 1, 2, 2, 4, 4]
 -- >>> lowerBound xs 1
--- Just 0
+-- 0
 --
 -- >>> lowerBound xs 2
--- Just 2
+-- 2
 --
 -- >>> lowerBound xs 3
--- Just 4
+-- 4
 --
 -- >>> lowerBound xs 4
--- Just 4
+-- 4
 --
--- Out of range values:
+-- Out of range values also return some index:
 --
 -- >>> lowerBound xs 0
--- Just 0
+-- 0
 --
 -- >>> lowerBound xs 5
--- Nothing
+-- 6
 --
--- @since 1.1.0.0
+-- @since 1.3.0.0
 {-# INLINE lowerBound #-}
-lowerBound :: (HasCallStack, VG.Vector v a, Ord a) => v a -> a -> Maybe Int
+lowerBound :: (HasCallStack, VG.Vector v a, Ord a) => v a -> a -> Int
 lowerBound vec = lowerBoundIn 0 (VG.length vec) vec
 
 -- | \(O(\log n)\) Computes the `lowerBound` for a slice of a vector within the interval \([l, r)\).
 --
 -- - The user predicate evaluates indices in \([l, r)\).
--- - The interval \([l, r)\) is silently clamped to ensure it remains within the bounds \([0, n)\).
 --
+-- ==== Constraints
+-- - \(0 \le l \lt n)
+-- - \(-1 \le r \le n)
+--
 -- ==== __Example__
 -- >>> import Data.Vector.Unboxed qualified as VU
 -- >>> let xs = VU.fromList [10, 10, 20, 20, 40, 40]
 -- >>> --                            *---*---*
 -- >>> lowerBoundIn 2 5 xs 10
--- Just 2
+-- 2
 --
 -- >>> lowerBoundIn 2 5 xs 20
--- Just 2
+-- 2
 --
 -- >>> lowerBoundIn 2 5 xs 30
--- Just 4
+-- 4
 --
 -- >>> lowerBoundIn 2 5 xs 40
--- Just 4
+-- 4
 --
 -- >>> lowerBoundIn 2 5 xs 50
--- Nothing
+-- 5
 --
--- @since 1.1.0.0
+-- @since 1.3.0.0
 {-# INLINE lowerBoundIn #-}
-lowerBoundIn :: (VG.Vector v a, Ord a) => Int -> Int -> v a -> a -> Maybe Int
-lowerBoundIn l_ r_ vec target
-  | ACIA.testInterval l r (VG.length vec) = bisectR l r $ \i -> VG.unsafeIndex vec i < target
-  | otherwise = Nothing
+lowerBoundIn :: (HasCallStack, VG.Vector v a, Ord a) => Int -> Int -> v a -> a -> Int
+lowerBoundIn l r vec target = maxRight l r $ \i -> vec VG.! i < target
   where
-    -- clamp
-    l = max 0 l_
-    r = min (VG.length vec) r_
+    !_ = ACIA.checkIntervalBounded "AtCoder.Extra.Bisect.lowerBoundIn" l r $ VG.length vec
 
--- | \(O(\log n)\) Returns the index of the first element \(x\) in the vector such that
--- \(x \gt x_0\), or `Nothing` if no such element exists.
+-- | \(O(\log n)\) Returns the maximum \(r\) where \(x \le x_i\) holds for \(i \in [0, r)\).
 --
 -- @
 -- Y Y Y Y Y N N N N N      Y: (<= x0)
 -- --------- *---------> X  N: (> x0)
---           R              R: returning point
+--           R              R: the right boundary point returned
 -- @
 --
 -- ==== __Example__
 -- >>> import Data.Vector.Unboxed qualified as VU
 -- >>> let xs = VU.fromList [10, 10, 20, 20, 40, 40]
 -- >>> upperBound xs 10
--- Just 2
+-- 2
 --
 -- >>> upperBound xs 20
--- Just 4
+-- 4
 --
 -- >>> upperBound xs 30
--- Just 4
+-- 4
 --
--- >>> upperBound xs 40
--- Nothing
+-- >>> upperBound xs 39
+-- 4
 --
 -- Out of range values:
 --
 -- >>> upperBound xs 0
--- Just 0
+-- 0
 --
--- >>> upperBound xs 50
--- Nothing
+-- >>> upperBound xs 40
+-- 6
 --
--- @since 1.1.0.0
+-- @since 1.3.0.0
 {-# INLINE upperBound #-}
-upperBound :: (HasCallStack, VG.Vector v a, Ord a) => v a -> a -> Maybe Int
+upperBound :: (HasCallStack, VG.Vector v a, Ord a) => v a -> a -> Int
 upperBound vec = upperBoundIn 0 (VG.length vec) vec
 
 -- | \(O(\log n)\) Computes the `upperBound` for a slice of a vector within the interval \([l, r)\).
@@ -185,115 +157,133 @@
 -- >>> let xs = VU.fromList [10, 10, 20, 20, 40, 40]
 -- >>> --                            *---*---*
 -- >>> upperBoundIn 2 5 xs 0
--- Just 2
+-- 2
 --
 -- >>> upperBoundIn 2 5 xs 10
--- Just 2
+-- 2
 --
 -- >>> upperBoundIn 2 5 xs 20
--- Just 4
+-- 4
 --
 -- >>> upperBoundIn 2 5 xs 30
--- Just 4
+-- 4
 --
 -- >>> upperBoundIn 2 5 xs 40
--- Nothing
+-- 5
 --
 -- >>> upperBoundIn 2 5 xs 50
--- Nothing
+-- 5
 --
--- @since 1.1.0.0
+-- @since 1.3.0.0
 {-# INLINE upperBoundIn #-}
-upperBoundIn :: (VG.Vector v a, Ord a) => Int -> Int -> v a -> a -> Maybe Int
-upperBoundIn l_ r_ vec target
-  | ACIA.testInterval l r (VG.length vec) = bisectR l r $ \i -> VG.unsafeIndex vec i <= target
-  | otherwise = Nothing
+upperBoundIn :: (HasCallStack, VG.Vector v a, Ord a) => Int -> Int -> v a -> a -> Int
+upperBoundIn l r vec target = maxRight l r $ \i -> vec VG.! i <= target
   where
-    -- clamp
-    l = max 0 l_
-    r = min (VG.length vec) r_
+    !_ = ACIA.checkIntervalBounded "AtCoder.Extra.Bisect.upperBoundIn" l r $ VG.length vec
 
 -- | \(O(\log n)\) Applies the bisection method on a half-open interval \([l, r)\) and returns the
--- left boundary point, or `Nothing` if no such point exists.
+-- right boundary point.
 --
 -- @
 -- Y Y Y Y Y N N N N N      Y: user predicate holds
--- --------* ----------> X  N: user predicate does not hold
---         L                L: the left boundary point returned
+-- --------- *---------> X  N: user predicate does not hold
+--           R              R: the right boundary point returned
 -- @
 --
 -- ==== __Example__
 -- >>> import Data.Vector.Unboxed qualified as VU
 -- >>> let xs = VU.fromList [10, 10, 20, 20, 30, 30]
 -- >>> let n = VU.length xs
--- >>> bisectL 0 n ((<= 20) . (xs VU.!))
--- Just 3
+-- >>> maxRight 0 n ((<= 20) . (xs VU.!))
+-- 4
 --
--- >>> bisectL 0 n ((<= 0) . (xs VU.!))
--- Nothing
+-- >>> maxRight 0 n ((<= 0) . (xs VU.!))
+-- 0
 --
--- >>> bisectL 0 n ((<= 100) . (xs VU.!))
--- Just 5
+-- >>> maxRight 0 n ((<= 100) . (xs VU.!))
+-- 6
 --
--- >>> bisectL 0 3 ((<= 20) . (xs VU.!))
--- Just 2
+-- >>> maxRight 0 3 ((<= 20) . (xs VU.!))
+-- 3
 --
--- @since 1.1.0.0
-{-# INLINE bisectL #-}
-bisectL :: (HasCallStack) => Int -> Int -> (Int -> Bool) -> Maybe Int
-bisectL l r p = runIdentity $ bisectLM l r (pure . p)
+-- @since 1.3.0.0
+{-# INLINE maxRight #-}
+maxRight ::
+  (HasCallStack) =>
+  -- | \(l\)
+  Int ->
+  -- | \(r\)
+  Int ->
+  -- | \(p\)
+  (Int -> Bool) ->
+  -- | Maximum \(r'\) where \(p\) holds for \([l, r')\).
+  Int
+maxRight l r p = runIdentity $ maxRightM l r (pure . p)
 
--- | \(O(\log n)\) Monadic variant of `bisectL`.
+-- | \(O(\log n)\) Monadic variant of `maxRight`.
 --
--- @since 1.1.0.0
-{-# INLINE bisectLM #-}
-bisectLM :: (HasCallStack, Monad m) => Int -> Int -> (Int -> m Bool) -> m (Maybe Int)
-bisectLM l r p
-  | l >= r = pure Nothing
-  | otherwise =
-      bisectLImpl p l r <&> \case
-        i | i == (l - 1) -> Nothing
-        i -> Just i
+-- @since 1.3.0.0
+{-# INLINE maxRightM #-}
+maxRightM :: (HasCallStack, Monad m) => Int -> Int -> (Int -> m Bool) -> m Int
+maxRightM l0 r0 p = bisectImpl (l0 - 1) r0 p
+  where
+    !_ = ACIA.checkInterval "AtCoder.Extra.Bisect.maxRightM" l0 r0
 
 -- | \(O(\log n)\) Applies the bisection method on a half-open interval \([l, r)\) and returns the
--- right boundary point, or `Nothing` if no such point exists.
---
+-- left boundary point.
 --
 -- @
--- Y Y Y Y Y N N N N N      Y: user predicate holds
--- --------- *---------> X  N: user predicate does not hold
---           R              R: the right boundary point returned
+-- N N N N N Y Y Y Y Y      Y: user predicate holds
+-- --------* ----------> X  N: user predicate does not hold
+--         L                L: the left boundary point returned
 -- @
 --
 -- ==== __Example__
 -- >>> import Data.Vector.Unboxed qualified as VU
 -- >>> let xs = VU.fromList [10, 10, 20, 20, 30, 30]
 -- >>> let n = VU.length xs
--- >>> bisectR 0 n ((<= 20) . (xs VU.!))
--- Just 4
---
--- >>> bisectR 0 n ((<= 0) . (xs VU.!))
--- Just 0
+-- >>> minLeft 0 n ((>= 20) . (xs VU.!))
+-- 2
 --
--- >>> bisectR 0 n ((<= 100) . (xs VU.!))
--- Nothing
+-- >>> minLeft 0 n ((>= 0) . (xs VU.!))
+-- 0
 --
--- >>> bisectR 0 4 ((<= 20) . (xs VU.!))
--- Nothing
+-- >>> minLeft 0 n ((>= 100) . (xs VU.!))
+-- 6
 --
--- @since 1.1.0.0
-{-# INLINE bisectR #-}
-bisectR :: (HasCallStack) => Int -> Int -> (Int -> Bool) -> Maybe Int
-bisectR l r p = runIdentity $ bisectRM l r (pure . p)
+-- @since 1.3.0.0
+{-# INLINE minLeft #-}
+minLeft ::
+  (HasCallStack) =>
+  -- | \(l\)
+  Int ->
+  -- | \(r\)
+  Int ->
+  -- | \(p\)
+  (Int -> Bool) ->
+  -- | Minimum \(l'\) where \(p\) holds for \([l', r)\)
+  Int
+minLeft l r p = runIdentity $ minLeftM l r (pure . p)
 
--- | \(O(\log n)\) Monadic variant of `bisectR`.
+-- | \(O(\log n)\) Monadic variant of `maxRight`.
 --
--- @since 1.1.0.0
-{-# INLINE bisectRM #-}
-bisectRM :: (HasCallStack, Monad m) => Int -> Int -> (Int -> m Bool) -> m (Maybe Int)
-bisectRM l r p
-  | l >= r = pure Nothing
-  | otherwise =
-      bisectRImpl p l r <&> \case
-        i | i == r -> Nothing
-        i -> Just i
+-- @since 1.3.0.0
+{-# INLINE minLeftM #-}
+minLeftM :: (HasCallStack, Monad m) => Int -> Int -> (Int -> m Bool) -> m Int
+minLeftM l r p = (+ 1) <$> bisectImpl r (l - 1) p
+  where
+    !_ = ACIA.checkInterval "AtCoder.Extra.Bisect.minLeftM" l r
+
+-- | Takes an open interval (l, r) or (r, l).
+{-# INLINE bisectImpl #-}
+bisectImpl :: (HasCallStack, Monad m) => Int -> Int -> (Int -> m Bool) -> m Int
+bisectImpl l0 r0 p = inner l0 r0
+  where
+    inner l r
+      | abs (r - l) <= 1 = pure r
+      | otherwise =
+          p mid >>= \case
+            True -> inner mid r
+            False -> inner l mid
+      where
+        mid = (l + r) `div` 2
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
@@ -56,7 +56,7 @@
     dijkstra,
     trackingDijkstra,
 
-    -- ** Bellman–ford algorithm
+    -- ** Bellman–Ford algorithm
 
     -- | - Vertex type is restricted to one-dimensional `Int`.
     bellmanFord,
@@ -74,7 +74,7 @@
 
     -- ** Path reconstruction
 
-    -- *** Single start point (root)
+    -- *** Single source point (root)
 
     -- | Functions for retrieving a path from a predecessor array where @-1@ represents none.
     constructPathFromRoot,
@@ -871,7 +871,7 @@
 -- -- | Option for `bellmanFord`.
 -- data BellmanFordPolicy = QuitOnNegaitveLoop | ContinueOnNegaitveLoop
 
--- | \(O(nm)\) Bellman–ford algorithm that returns a distance array, or `Nothing` on negative loop
+-- | \(O(nm)\) Bellman–Ford algorithm that returns a distance array, or `Nothing` on negative loop
 -- detection. Vertices are one-dimensional.
 --
 -- ==== __Example__
@@ -907,7 +907,7 @@
   (!dist, !_) <- bellmanFordImpl False nVerts gr undefW source
   pure dist
 
--- | \(O(nm)\) Bellman–ford algorithm that returns a distance array and a predecessor array, or
+-- | \(O(nm)\) Bellman–Ford algorithm that returns a distance array and a predecessor array, or
 -- `Nothing` on negative loop detection. Vertices are one-dimensional.
 --
 -- ==== __Example__
@@ -1328,7 +1328,7 @@
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
--- - There must be a path from the root to the @end@ vertex, otherwise the returned path is not
+-- - There must be a path from the root to the @sink@ vertex, otherwise the returned path is not
 -- connected to the root.
 --
 -- @since 1.2.4.0
@@ -1340,7 +1340,7 @@
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
--- - There must be a path from the root to the @end@ vertex, otherwise the returned path is not
+-- - There must be a path from the root to the @sink@ vertex, otherwise the returned path is not
 -- connected to the root.
 --
 -- @since 1.2.4.0
@@ -1352,11 +1352,11 @@
     f v = Just (v, parents VG.! v)
 
 -- | \(O(n)\) Given a NxN predecessor matrix (created with `trackingFloydWarshall`), retrieves a
--- path from the root to an end vertex.
+-- path from the root to a sink vertex.
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
--- - There must be a path from the root to the @end@ vertex, otherwise the returned path is not
+-- - There must be a path from the root to the @sink@ vertex, otherwise the returned path is not
 -- connected to the root.
 --
 -- @since 1.2.4.0
@@ -1365,20 +1365,20 @@
   (HasCallStack) =>
   -- | Predecessor matrix.
   VU.Vector Int ->
-  -- | Start vertex.
+  -- | Source vertex.
   Int ->
-  -- | End vertex.
+  -- | Sink vertex.
   Int ->
   -- | Path.
   VU.Vector Int
-constructPathFromRootMat parents start = VU.reverse . constructPathToRootMat parents start
+constructPathFromRootMat parents source = VU.reverse . constructPathToRootMat parents source
 
 -- | \(O(n)\) Given a NxN predecessor matrix(created with `trackingFloydWarshall`), retrieves a
 -- path from a vertex to the root.
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
--- - There must be a path from the root to the @end@ vertex, otherwise the returned path is not
+-- - There must be a path from the root to the @sink@ vertex, otherwise the returned path is not
 -- connected to the root.
 --
 -- @since 1.2.4.0
@@ -1387,26 +1387,26 @@
   (HasCallStack) =>
   -- | Predecessor matrix.
   VU.Vector Int ->
-  -- | Start vertex.
+  -- | Source vertex.
   Int ->
-  -- | End vertex.
+  -- | Sink vertex.
   Int ->
   -- | Path.
   VU.Vector Int
-constructPathToRootMat parents start end =
-  let parents' = VU.take n $ VU.drop (n * start) parents
-   in constructPathToRoot parents' end
+constructPathToRootMat parents source sink =
+  let parents' = VU.take n $ VU.drop (n * source) parents
+   in constructPathToRoot parents' sink
   where
     -- Assuming `n < 2^32`, it should always be correct:
     -- 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 an end vertex.
+-- path from the root to a sink vertex.
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
--- - There must be a path from the root to the @end@ vertex, otherwise the returned path is not
+-- - There must be a path from the root to the @nd@ vertex, otherwise the returned path is not
 -- connected to the root.
 --
 -- @since 1.2.4.0
@@ -1415,20 +1415,20 @@
   (HasCallStack, PrimMonad m) =>
   -- | Predecessor matrix.
   VUM.MVector (PrimState m) Int ->
-  -- | Start vertex.
+  -- | Source vertex.
   Int ->
-  -- | End vertex.
+  -- | Sink vertex.
   Int ->
   -- | Path.
   m (VU.Vector Int)
-constructPathFromRootMatM parents start = (VU.reverse <$>) . constructPathToRootMatM parents start
+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.
 --
 -- ==== Constraints
 -- - The path must not make a cycle, otherwise this function loops forever.
--- - There must be a path from the root to the @end@ vertex, otherwise the returned path is not
+-- - There must be a path from the root to the @sink@ vertex, otherwise the returned path is not
 -- connected to the root.
 --
 -- @since 1.2.4.0
@@ -1437,12 +1437,12 @@
   (HasCallStack, PrimMonad m) =>
   -- | Predecessor matrix.
   VUM.MVector (PrimState m) Int ->
-  -- | Start vertex.
+  -- | Source vertex.
   Int ->
-  -- | End vertex.
+  -- | Sink vertex.
   Int ->
   -- | Path.
   m (VU.Vector Int)
-constructPathToRootMatM parents start end = stToPrim $ do
+constructPathToRootMatM parents source sink = stToPrim $ do
   parents' <- VU.unsafeFreeze parents
-  pure $ constructPathToRootMat parents' start end
+  pure $ constructPathToRootMat parents' source sink
diff --git a/src/AtCoder/Extra/Math.hs b/src/AtCoder/Extra/Math.hs
--- a/src/AtCoder/Extra/Math.hs
+++ b/src/AtCoder/Extra/Math.hs
@@ -2,12 +2,7 @@
 --
 -- @since 1.0.0.0
 module AtCoder.Extra.Math
-  ( -- * Re-exports from the internal math module
-    isPrime32,
-    ACIM.invGcd,
-    primitiveRoot32,
-
-    -- * Prime numbers and divisors
+  ( -- * Prime numbers and divisors
     primes,
     isPrime,
     primeFactors,
@@ -15,6 +10,9 @@
     divisors,
     divisorsUnsorted,
 
+    -- * PrimitiveRoot
+    primitiveRoot,
+
     -- * Binary exponentiation
 
     -- | ==== __Examples__
@@ -35,8 +33,6 @@
 where
 
 import AtCoder.Extra.Math.Montgomery64 qualified as M64
-import AtCoder.Internal.Assert qualified as ACIA
-import AtCoder.Internal.Math qualified as ACIM
 import Control.Monad (unless, when)
 import Data.Bit (Bit (..))
 import Data.Bits (bit, countTrailingZeros, (.<<.), (.>>.))
@@ -51,33 +47,6 @@
 import GHC.Stack (HasCallStack)
 import System.Random
 
--- | \(O(k \log^3 n) (k = 3)\). Returns whether the given `Int` value is a prime number.
---
--- ==== Constraints
--- - \(n < 4759123141 (2^{32} < 4759123141)\), otherwise the return value can lie
---   (see [Wikipedia](https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Testing_against_small_sets_of_bases)).
---
---
--- @since 1.1.0.0
-{-# INLINE isPrime32 #-}
-isPrime32 :: (HasCallStack) => Int -> Bool
-isPrime32 x = ACIM.isPrime x
-  where
-    !_ = ACIA.runtimeAssert (x < 4759123141) $ "AtCoder.Extra.Math.isPrime32: given too large number `" ++ show x ++ "`"
-
--- | Returns the primitive root of the given `Int`.
---
--- ==== Constraints
--- - The input must be a prime number.
--- - The input must be less than \(2^32\).
---
--- @since 1.2.0.0
-{-# INLINE primitiveRoot32 #-}
-primitiveRoot32 :: (HasCallStack) => Int -> Int
-primitiveRoot32 x = ACIM.primitiveRoot x
-  where
-    !_ = ACIA.runtimeAssert (x < (1 .>>. 32)) $ "AtCoder.Extra.Math.primitiveRoot32: given too large number `" ++ show x ++ "`"
-
 -- | \(O(n \log \log n)\) Creates an array of prime numbers up to the given limit, using Sieve of
 -- Eratosthenes.
 --
@@ -88,9 +57,13 @@
 -- - The upper limit must be less than or equal to \(2^{30} (\gt 10^9)\), otherwise the returned
 -- prime table is incorrect.
 --
+-- ==== __Example__
+-- >>> primes 100
+-- [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
+--
 -- @since 1.2.6.0
 {-# INLINEABLE primes #-}
-primes :: Int -> VU.Vector Int
+primes :: (HasCallStack) => Int -> VU.Vector Int
 primes upperLimit
   | upperLimit <= 1 = VU.empty
   | otherwise = VU.create $ do
@@ -180,9 +153,13 @@
 -- | \(O(w \log^3 n)\) Miller–Rabin primality test, where \(w = 3\) for \(x \lt 2^{32}\) and
 -- \(w = 7\) for \(x \ge 3^{32}\).
 --
+-- ==== __Example__
+-- >>> isPrime 100055128505716009
+-- True
+--
 -- @since 1.2.6.0
 {-# INLINEABLE isPrime #-}
-isPrime :: Int -> Bool
+isPrime :: (HasCallStack) => Int -> Bool
 isPrime x
   | x <= 1 = False
   -- Up to 11^2:
@@ -279,6 +256,13 @@
 -- ==== Constraints
 -- - \(x \ge 1\)
 --
+-- ==== __Example__
+-- >>> primeFactors 180
+-- [(2,2),(3,2),(5,1)]
+--
+-- >>> primeFactors 123123123123123123
+-- [(3,2),(7,1),(11,1),(13,1),(19,1),(41,1),(52579,1),(333667,1)]
+--
 -- @since 1.2.6.0
 {-# INLINE primeFactors #-}
 primeFactors :: (HasCallStack) => Int -> VU.Vector (Int, Int)
@@ -339,9 +323,13 @@
 -- ==== Constraints
 -- - \(x \ge 1\)
 --
+-- ==== __Example__
+-- >>> divisors 180
+-- [1,2,3,4,5,6,9,10,12,15,18,20,30,36,45,60,90,180]
+--
 -- @since 1.2.6.0
 {-# INLINE divisors #-}
-divisors :: Int -> VU.Vector Int
+divisors :: (HasCallStack) => Int -> VU.Vector Int
 -- TODO: use intro sort?
 divisors = VU.modify VAR.sort . divisorsUnsorted
 
@@ -352,7 +340,7 @@
 --
 -- @since 1.2.6.0
 {-# INLINEABLE divisorsUnsorted #-}
-divisorsUnsorted :: Int -> VU.Vector Int
+divisorsUnsorted :: (HasCallStack) => Int -> VU.Vector Int
 divisorsUnsorted x = VU.create $ do
   vec <- VUM.unsafeNew nDivisors
   VGM.write vec 0 1
@@ -379,7 +367,46 @@
     (!_, !ns) = VU.unzip pns
     nDivisors = VU.foldl' (\ !acc n -> acc * (n + 1)) (1 :: Int) ns
 
--- | Calculates \(x^n\) with custom multiplication operator using the binary exponentiation
+-- | Returns a primitive root of module \(p\), where \(p\) is a prime number.
+--
+-- ==== Constraints
+-- - \(p\) must be a prime number.
+--
+-- ==== __Example__
+-- >>> primitiveRoot 999999999999999989
+-- 833278905416200545
+--
+-- >>> primitiveRoot 100055128505716009
+-- 40765942246299710
+--
+-- @since 1.2.7.0
+{-# INLINEABLE primitiveRoot #-}
+primitiveRoot :: (HasCallStack) => Int -> Int
+primitiveRoot x
+  | not (isPrime x) = error $ "AtCoder.Extra.Math.primitiveRoot: give non-prime value `" ++ show x ++ "`"
+  | x == 2 = 1
+primitiveRoot x = tryRandom $ mkStdGen 123456789
+  where
+    !x64 :: Word64 = fromIntegral x
+    !mont = M64.fromVal x64
+    (!ps, !_) = VU.unzip $ primeFactorsUnsorted (x - 1)
+    -- g s.t. for all p_i. g^n mod p_i /= 1
+    test :: (HasCallStack) => Word64 -> Int -> Bool
+    test g p =
+      let !n = (x - 1) `div` p
+       in (/= 1)
+            . M64.decode mont
+            . power (M64.mulMod mont) n
+            . M64.encode mont
+            $ fromIntegral g
+    tryRandom :: (HasCallStack) => StdGen -> Int
+    tryRandom !gen
+      | VU.all (test rnd) ps = fromIntegral rnd
+      | otherwise = tryRandom gen'
+      where
+        (!rnd, !gen') = uniformR (1, x64 - 1) gen
+
+-- | Calculates \(f^n(x)\) with custom multiplication operator using the binary exponentiation
 -- technique.
 --
 -- The internal implementation is taken from @Data.Semigroup.stimes@, but `power` uses strict
@@ -393,7 +420,7 @@
 --
 -- @since 1.0.0.0
 {-# INLINE power #-}
-power :: (a -> a -> a) -> Int -> a -> a
+power :: (HasCallStack) => (a -> a -> a) -> Int -> a -> a
 power op n0 x1
   | n0 <= 0 = errorWithoutStackTrace "AtCoder.Extra.Math.power: positive multiplier expected"
   | otherwise = f x1 n0
@@ -407,6 +434,22 @@
       | n == 1 = x `op` z
       | otherwise = g (x `op` x) (n .>>. 1) (x `op` z)
 
+-- TODO: powMod for arbitrary modulus value (needs Barrett64)
+
+-- -- | \(O(\log n)\) One-shot \(x^n \bmod m\) calculation.
+-- --
+-- -- @since 1.2.7.0
+-- {-# INLINE powMod #-}
+-- powMod :: (HasCallStack) => Int -> Int -> Int -> Int
+-- powMod x n m =
+--   fromIntegral
+--     . M64.decode mont
+--     . power (M64.mulMod mont) n
+--     . M64.encode mont
+--     $ fromIntegral x
+--   where
+--     !mont = M64.fromVal $! fromIntegral m
+
 -- | Strict variant of @Data.Semigroup.stimes@.
 --
 -- ==== Complexity
@@ -417,7 +460,7 @@
 --
 -- @since 1.0.0.0
 {-# INLINE stimes' #-}
-stimes' :: (Semigroup a) => Int -> a -> a
+stimes' :: (HasCallStack) => (Semigroup a) => Int -> a -> a
 stimes' = power (<>)
 
 -- | Strict variant of @Data.Monoid.mtimes@.
@@ -430,8 +473,8 @@
 --
 -- @since 1.0.0.0
 {-# INLINE mtimes' #-}
-mtimes' :: (Monoid a) => Int -> a -> a
+mtimes' :: (HasCallStack) => (Monoid a) => Int -> a -> a
 mtimes' n x = case compare n 0 of
-  LT -> errorWithoutStackTrace "AtCoder.Extra.Math.mtimes': non-negative multiplier expected"
+  LT -> error "AtCoder.Extra.Math.mtimes': non-negative multiplier expected"
   EQ -> mempty
   GT -> power (<>) n x
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
@@ -66,7 +66,7 @@
 --
 -- @since 1.2.6.0
 {-# NOINLINE new #-}
-new :: forall a. (KnownNat a) => Proxy# a -> Montgomery64
+new :: forall a. (HasCallStack, KnownNat a) => Proxy# a -> Montgomery64
 -- FIXME: test allocated once
 new p = fromVal . fromIntegral $! natVal' p
 
@@ -78,7 +78,7 @@
 --
 -- @since 1.2.6.0
 {-# INLINE fromVal #-}
-fromVal :: Word64 -> Montgomery64
+fromVal :: (HasCallStack) => Word64 -> Montgomery64
 fromVal m =
   let !m128 :: Word128 = fromIntegral m
       !n2 = word128Lo64 $ (-m128) `mod` m128
diff --git a/src/AtCoder/Extra/ModInt64.hs b/src/AtCoder/Extra/ModInt64.hs
--- a/src/AtCoder/Extra/ModInt64.hs
+++ b/src/AtCoder/Extra/ModInt64.hs
@@ -3,8 +3,11 @@
 {-# LANGUAGE MagicHash #-}
 {-# LANGUAGE TypeFamilies #-}
 
--- | @ModInt@ for 64 bit modulus values.
+-- | @ModInt@ for 64 bit modulus values with Montgomery modular multiplication.
 --
+-- ==== Constraints
+-- - The modulus value should be an odd number, otherwise it would be too slow.
+--
 -- @since 1.2.6.0
 module AtCoder.Extra.ModInt64
   ( -- * ModInt64
@@ -49,6 +52,9 @@
 
 -- | `Word64` value that treats the modular arithmetic.
 --
+-- ==== Constraints
+-- - The modulus value should be an odd number, otherwise it would be too slow.
+--
 -- @since 1.2.6.0
 newtype ModInt64 a = ModInt64
   { -- | Montgomery form of the value. Use `val` to retrieve the value.
@@ -225,13 +231,7 @@
 
 -- | @since 1.2.6.0
 instance (KnownNat p) => Integral (ModInt64 p) where
-  -- FIXME: THIS IS COMPLETELY WRONG. Compare with `ModInt`.
   {-# INLINE quotRem #-}
-  -- quotRem x y =
-  --   let !x' = val x
-  --       !y' = val y
-  --       (!q, !r) = x' `quotRem` y'
-  --    in (new q, new r)
   quotRem x y = (x / y, x - x / y * y)
   {-# INLINE toInteger #-}
   toInteger = toInteger . val
diff --git a/src/AtCoder/Extra/SegTree2d.hs b/src/AtCoder/Extra/SegTree2d.hs
--- a/src/AtCoder/Extra/SegTree2d.hs
+++ b/src/AtCoder/Extra/SegTree2d.hs
@@ -68,7 +68,7 @@
 import Control.Monad.ST (ST)
 import Data.Bits
 import Data.Foldable (for_)
-import Data.Maybe (fromJust, fromMaybe)
+import Data.Maybe (fromMaybe)
 import Data.Vector.Algorithms.Intro qualified as VAI
 import Data.Vector.Generic qualified as VG
 import Data.Vector.Generic.Mutable qualified as VGM
@@ -317,7 +317,7 @@
   let nxSt = VU.length dictXSt
   let logSt = countTrailingZeros $ ACIB.bitCeil (nxSt + 1)
   let sizeSt = bit logSt
-  let compressedXs = VU.map (fromJust . lowerBound dictXSt) xs
+  let compressedXs = VU.map (lowerBound dictXSt) xs
 
   -- TODO: what is this?
   let indptrSt = VU.create $ do
@@ -398,13 +398,13 @@
 {-# INLINEABLE prodST #-}
 prodST :: forall s a. (HasCallStack, Monoid a, VU.Unbox a) => SegTree2d s a -> Int -> Int -> Int -> Int -> ST s a
 prodST seg@SegTree2d {..} lx rx ly ry = do
-  let a0 = fromMaybe (VU.length allYSt) $ lowerBound allYSt ly
-  let b0 = fromMaybe (VU.length allYSt) $ lowerBound allYSt ry
+  let a0 = lowerBound allYSt ly
+  let b0 = lowerBound allYSt ry
   dfs mempty 1 0 sizeSt a0 b0
   where
     !_ = ACIA.runtimeAssert (lx <= rx && ly <= ry) "AtCoder.Extra.SegTree2d.prodST: given invalid rectangle"
-    !l0 = fromMaybe (VU.length dictXSt) $ lowerBound dictXSt lx
-    !r0 = fromMaybe (VU.length dictXSt) $ lowerBound dictXSt rx
+    !l0 = lowerBound dictXSt lx
+    !r0 = lowerBound dictXSt rx
     dfs :: a -> Int -> Int -> Int -> Int -> Int -> ST s a
     dfs !res i l r a b
       -- empty rect
@@ -450,13 +450,13 @@
 {-# INLINEABLE countST #-}
 countST :: forall s a. (HasCallStack, Monoid a, VU.Unbox a) => SegTree2d s a -> Int -> Int -> Int -> Int -> ST s Int
 countST SegTree2d {..} lx rx ly ry = do
-  let a0 = fromMaybe (VU.length allYSt) $ lowerBound allYSt ly
-  let b0 = fromMaybe (VU.length allYSt) $ lowerBound allYSt ry
+  let a0 = lowerBound allYSt ly
+  let b0 = lowerBound allYSt ry
   dfs 0 1 0 sizeSt a0 b0
   where
     !_ = ACIA.runtimeAssert (lx <= rx && ly <= ry) "AtCoder.Extra.SegTree2d.countST: given invalid rectangle"
-    !l0 = fromMaybe (VU.length dictXSt) $ lowerBound dictXSt lx
-    !r0 = fromMaybe (VU.length dictXSt) $ lowerBound dictXSt rx
+    !l0 = lowerBound dictXSt lx
+    !r0 = lowerBound dictXSt rx
     dfs :: Int -> Int -> Int -> Int -> Int -> Int -> ST s Int
     dfs (res :: Int) i l r a b
       -- empty rect
diff --git a/src/AtCoder/Extra/Seq.hs b/src/AtCoder/Extra/Seq.hs
--- a/src/AtCoder/Extra/Seq.hs
+++ b/src/AtCoder/Extra/Seq.hs
@@ -102,10 +102,6 @@
 
     -- * Handle (re-exports)
     Handle (..),
-    -- TODO: hide
-    newHandle,
-    nullHandle,
-    invalidateHandle,
 
     -- * Re-exports
     SegAct (..),
@@ -173,7 +169,7 @@
   )
 where
 
-import AtCoder.Extra.Pool (Handle (..), invalidateHandle, newHandle, nullHandle)
+import AtCoder.Extra.Pool (Handle (..), newHandle)
 import AtCoder.Extra.Pool qualified as P
 import AtCoder.Extra.Seq.Raw (Seq (..))
 import AtCoder.Extra.Seq.Raw qualified as Seq
diff --git a/src/AtCoder/Extra/Seq/Map.hs b/src/AtCoder/Extra/Seq/Map.hs
--- a/src/AtCoder/Extra/Seq/Map.hs
+++ b/src/AtCoder/Extra/Seq/Map.hs
@@ -144,7 +144,7 @@
 new n = stToPrim $ do
   seqMap <- Seq.new n
   kMap <- VUM.unsafeNew n
-  rootMap <- Seq.newHandle P.undefIndex
+  rootMap <- P.newHandle P.undefIndex
   pure Map {..}
 
 -- | \(O(n \log n)\) Creates a new `Map` of capacity \(n\) with initial values. Always prefer `build` to
diff --git a/src/AtCoder/Extra/SqrtDecomposition.hs b/src/AtCoder/Extra/SqrtDecomposition.hs
--- a/src/AtCoder/Extra/SqrtDecomposition.hs
+++ b/src/AtCoder/Extra/SqrtDecomposition.hs
@@ -3,13 +3,13 @@
 -- interval query block by block, typically in \(O(\sqrt n)\) time, where a whole block processing
 -- take \(O(1)\) time and partial block processing take \(O(\sqrt n)\) time.
 --
--- For simplicity, in this document, we assume that highder order functions applided to an entier
+-- For simplicity, in this document, we assume that higher order functions applied to an entire
 -- block (@readFull@ and @actFull@) work in \(O(1)\) time, and those applied to a part of block work
 -- in \(O(\sqrt n)\) time. In total, \(q\) query processing takes \(O(q \sqrt n)\) time. Note that
 -- it's a rather large number and often requires performance tuning.
 --
 -- ==== Lazy propagation
--- Typiaclly, an action to a whole block can be delayed; store the aggregation value for the block,
+-- Typically, an action to a whole block can be delayed; store the aggregation value for the block,
 -- delay the internal sequence update, and restore them when part of the block is accessed. Such
 -- lazy propagation should be handled on the user side on partial block access functions
 -- (@foldPart@ or @actPart@) are called.
@@ -65,7 +65,7 @@
       when (remR > 0) $ do
         actPart ir (r - remR) r
 
--- | \(O(\sqrt n)\) Runs user function for each block and concatanate their monoid output.
+-- | \(O(\sqrt n)\) Runs user function for each block and concatenates their monoid output.
 --
 -- ==== Constraints
 -- - \(l \le r\)
@@ -90,7 +90,7 @@
   m a
 foldMapM blockLen = foldMapWithM blockLen (<>)
 
--- | \(O(\sqrt n)\) Runs user function for each block and concatanates their output with user
+-- | \(O(\sqrt n)\) Runs user function for each block and concatenates their output with user
 -- function.
 --
 -- ==== Constraints
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
@@ -268,10 +268,10 @@
 -- :}
 -- [0,1,4,1,0]
 --
--- @since 1.1.0.0
+-- @since 1.3.0.0
 {-# INLINE scan #-}
 scan ::
-  (VU.Unbox w, VG.Vector v a) =>
+  (VU.Unbox w, VU.Unbox a) =>
   -- | The number of vertices.
   Int ->
   -- | Graph as a function.
@@ -285,9 +285,9 @@
   -- | Root vertex.
   Int ->
   -- | Tree scanning result from a root vertex.
-  v a
-scan n tree acc0At toF act root = VG.create $ do
-  dp <- VGM.unsafeNew n
+  VU.Vector a
+scan n tree acc0At toF act root = VU.create $ do
+  dp <- VUM.unsafeNew n
   !_ <- foldImpl tree acc0At toF act root $ \v a -> do
     VGM.unsafeWrite dp v a
   pure dp
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
@@ -3,23 +3,13 @@
 -- @since 1.2.2.0
 module AtCoder.Extra.Vector
   ( argsort,
-    unsafePermuteInPlace,
-    unsafePermuteInPlaceST,
   )
 where
 
-import AtCoder.Internal.Assert qualified as ACIA
-import Control.Monad (unless)
-import Control.Monad.Primitive (PrimMonad, PrimState, stToPrim)
-import Control.Monad.ST (ST)
 import Data.Vector.Algorithms.Intro qualified as VAI
 import Data.Vector.Generic qualified as VG
-import Data.Vector.Generic.Mutable qualified as VGM
 import Data.Vector.Unboxed qualified as VU
 
--- TODO: test `unsafePermuteInPlace`
--- TODO: is `unsafePermuteInPlace` fast enough as specialized one?
-
 -- | \(O(n \log n)\) Returns indices of the vector, stably sorted by their value.
 --
 -- ==== Example
@@ -27,6 +17,8 @@
 -- >>> import Data.Vector.Unboxed qualified as VU
 -- >>> argsort $ VU.fromList [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
 argsort xs =
@@ -39,29 +31,5 @@
     )
     $ VU.generate (VU.length xs) id
 
--- | \(O(n)\) Applies a permutation to a mutable vector in-place.
---
--- ==== Constraints
--- - The index array must be a permutation (0-based).
-{-# INLINE unsafePermuteInPlace #-}
-unsafePermuteInPlace :: (PrimMonad m, VGM.MVector v a) => v (PrimState m) a -> VU.Vector Int -> m ()
-unsafePermuteInPlace vec is = stToPrim $ unsafePermuteInPlaceST vec is
-
--- | \(O(n)\) Applies a permutation to a mutable vector in-place.
---
--- ==== Constraints
--- - The index array must be a permutation (0-based).
-{-# INLINEABLE unsafePermuteInPlaceST #-}
-unsafePermuteInPlaceST :: (VGM.MVector v a) => v s a -> VU.Vector Int -> ST s ()
-unsafePermuteInPlaceST vec is = do
-  let !_ = ACIA.runtimeAssert (VGM.length vec == VG.length is) "AtCoder.Extra.Vector.unsafePermuteInPlaceST: the length of the index array must be equal to the length of the permuted vector"
-  let inner i lastX = do
-        VGM.unsafeWrite vec i lastX
-        unless (i == 0) $ do
-          let i0' = VG.unsafeIndex is i
-          lastX' <- VGM.unsafeRead vec i
-          inner i0' lastX'
+-- TODO: maybe add lexical permutations, combinations and subsequences.
 
-  let i0' = VG.unsafeIndex is 0
-  x0' <- VGM.unsafeRead vec 0
-  inner i0' x0'
diff --git a/src/AtCoder/Extra/WaveletMatrix.hs b/src/AtCoder/Extra/WaveletMatrix.hs
--- a/src/AtCoder/Extra/WaveletMatrix.hs
+++ b/src/AtCoder/Extra/WaveletMatrix.hs
@@ -65,23 +65,23 @@
 import AtCoder.Extra.Bisect
 import AtCoder.Extra.WaveletMatrix.Raw qualified as Rwm
 import Control.Monad
-import Data.Maybe (fromJust, fromMaybe)
 import Data.Vector.Algorithms.Intro qualified as VAI
 import Data.Vector.Generic qualified as VG
 import Data.Vector.Unboxed qualified as VU
+import GHC.Stack (HasCallStack)
 
 -- | A static Wavelet Matrix.
 --
--- @since 1.1.0.0
+-- @since 1.3.0.0
 data WaveletMatrix = WaveletMatrix
   { -- | The internal wavelet matrix, where index compression is not automatically performed.
     --
-    -- @since 1.1.0.0
-    rawWM :: !Rwm.RawWaveletMatrix,
+    -- @since 1.3.0.0
+    rawWm :: !Rwm.RawWaveletMatrix,
     -- | Index compression dictionary.
     --
-    -- @since 1.1.0.0
-    xDictWM :: !(VU.Vector Int)
+    -- @since 1.3.0.0
+    yDictWm :: !(VU.Vector Int)
   }
 
 -- | \(O(n \log n)\) Creates a `WaveletMatrix` from an array \(a\).
@@ -90,9 +90,9 @@
 {-# INLINE build #-}
 build :: VU.Vector Int -> WaveletMatrix
 build ys =
-  let !xDictWM = VU.uniq $ VU.modify (VAI.sortBy compare) ys
-      !ys' = VU.map (fromJust . lowerBound xDictWM) ys
-      !rawWM = Rwm.build (VG.length ys) ys'
+  let !yDictWm = VU.uniq $ VU.modify (VAI.sortBy compare) ys
+      !ys' = VU.map (lowerBound yDictWm) ys
+      !rawWm = Rwm.build (VG.length ys) ys'
    in WaveletMatrix {..}
 
 -- | \(O(\log |S|)\) Returns \(a[k]\) or `Nothing` if the index is out of the bounds. Try to use the
@@ -100,8 +100,8 @@
 --
 -- @since 1.1.0.0
 {-# INLINE access #-}
-access :: WaveletMatrix -> Int -> Maybe Int
-access WaveletMatrix {..} i = (xDictWM VG.!) <$> Rwm.access rawWM i
+access ::  WaveletMatrix -> Int -> Maybe Int
+access WaveletMatrix {..} i = (yDictWm VG.!) <$> Rwm.access rawWm i
 
 -- | \(O(\log |S|)\) Returns the number of \(y\) in \([l, r)\).
 --
@@ -140,12 +140,12 @@
 rankBetween WaveletMatrix {..} l r y1 y2
   | not $ 0 <= l && l < r && r <= n = 0
   | y1' >= y2' = 0
-  | otherwise = Rwm.rankBetween rawWM l r y1' y2'
+  | otherwise = Rwm.rankBetween rawWm l r y1' y2'
   where
     -- Handles the case @yl@ or  @yr@ is not in the dict
-    n = Rwm.lengthRwm rawWM
-    y1' = fromMaybe n (bisectR 0 (VG.length xDictWM) ((< y1) . VG.unsafeIndex xDictWM))
-    y2' = maybe (-1) (+ 1) (bisectL 0 (VG.length xDictWM) ((< y2) . VG.unsafeIndex xDictWM))
+    n = Rwm.lengthRwm rawWm
+    y1' = lowerBound yDictWm y1
+    y2' = lowerBound yDictWm y2
 
 -- | \(O(\log |S|)\) Returns the index of the first \(y\) in \(a\), or `Nothing` if \(y\) is
 -- not found.
@@ -170,11 +170,12 @@
   -- | The index of \(k\)-th \(y\)
   Maybe Int
 selectKth WaveletMatrix {..} k y = do
-  i <- lowerBound xDictWM y
+  let !i = lowerBound yDictWm y
+  guard $ i < VG.length yDictWm
   -- TODO: we don't need such an explicit branch?
-  let !y' = xDictWM VG.! i
+  let !y' = yDictWm VG.! i
   guard $ y' == y
-  Rwm.selectKth rawWM k i
+  Rwm.selectKth rawWm k i
 
 -- | \(O(\log |S|)\) Given an interval \([l, r)\), it returns the index of the first occurrence
 -- (0-based) of \(y\) in the sequence, or `Nothing` if no such occurrence exists.
@@ -200,6 +201,7 @@
 -- @since 1.1.0.0
 {-# INLINEABLE selectKthIn #-}
 selectKthIn ::
+  (HasCallStack) =>
   -- | A wavelet matrix
   WaveletMatrix ->
   -- | \(l\)
@@ -213,11 +215,12 @@
   -- | The index of the \(k\)-th \(y\) in \([l, r)\).
   Maybe Int
 selectKthIn WaveletMatrix {..} l r k y = do
-  i <- lowerBound xDictWM y
+  let !i = lowerBound yDictWm y
+  guard $ i < VG.length yDictWm
   -- TODO: we don't need such an explicit branch?
-  let !y' = xDictWM VG.! i
+  let !y' = yDictWm VG.! i
   guard $ y' == y
-  Rwm.selectKthIn rawWM l r k i
+  Rwm.selectKthIn rawWm l r k i
 
 -- | \(O(\log |S|)\) Given the interval \([l, r)\), returns the index of the \(k\)-th (0-based)
 -- largest value. Note that duplicated values are treated as distinct occurrences.
@@ -225,6 +228,7 @@
 -- @since 1.1.0.0
 {-# INLINEABLE kthLargestIn #-}
 kthLargestIn ::
+  (HasCallStack) =>
   -- | A wavelet matrix
   WaveletMatrix ->
   -- | \(l\)
@@ -236,7 +240,7 @@
   -- | \(k\)-th largest \(y\) in \([l, r)\)
   Maybe Int
 kthLargestIn WaveletMatrix {..} l r k
-  | Just !y <- Rwm.kthLargestIn rawWM l r k = Just $ xDictWM VG.! y
+  | Just !y <- Rwm.kthLargestIn rawWm l r k = Just $ yDictWm VG.! y
   | otherwise = Nothing
 
 -- | \(O(\log |S|)\) Given the interval \([l, r)\), returns both the index and the value of the
@@ -256,7 +260,7 @@
   -- | \((i, y)\) for \(k\)-th largest \(y\) in \([l, r)\)
   Maybe (Int, Int)
 ikthLargestIn WaveletMatrix {..} l r k
-  | Just (!i, !y) <- Rwm.ikthLargestIn rawWM l r k = Just (i, xDictWM VG.! y)
+  | Just (!i, !y) <- Rwm.ikthLargestIn rawWm l r k = Just (i, yDictWm VG.! y)
   | otherwise = Nothing
 
 -- | \(O(\log |S|)\) Given the interval \([l, r)\), returns the index of the \(k\)-th (0-based)
@@ -276,7 +280,7 @@
   -- | \(k\)-th largest \(y\) in \([l, r)\)
   Maybe Int
 kthSmallestIn WaveletMatrix {..} l r k
-  | Just !y <- Rwm.kthSmallestIn rawWM l r k = Just $ xDictWM VG.! y
+  | Just !y <- Rwm.kthSmallestIn rawWm l r k = Just $ yDictWm VG.! y
   | otherwise = Nothing
 
 -- | \(O(\log |S|)\) Given the interval \([l, r)\), returns both the index and the value of the
@@ -295,7 +299,7 @@
   -- | \((i, y)\) for \(k\)-th largest \(y\) in \([l, r)\)
   Maybe (Int, Int)
 ikthSmallestIn WaveletMatrix {..} l r k
-  | Just (!i, !y) <- Rwm.ikthSmallestIn rawWM l r k = Just (i, xDictWM VG.! y)
+  | Just (!i, !y) <- Rwm.ikthSmallestIn rawWm l r k = Just (i, yDictWm VG.! y)
   | otherwise = Nothing
 
 -- | \(O(\log |S|)\)
@@ -304,7 +308,7 @@
 {-# INLINE unsafeKthSmallestIn #-}
 unsafeKthSmallestIn :: WaveletMatrix -> Int -> Int -> Int -> Int
 unsafeKthSmallestIn WaveletMatrix {..} l r k =
-  xDictWM VG.! Rwm.unsafeKthSmallestIn rawWM l r k
+  yDictWm VG.! Rwm.unsafeKthSmallestIn rawWm l r k
 
 -- | \(O(\log |S|)\) Looks up the maximum \(y\) in \([l, r) \times (-\infty, y_0]\).
 --
@@ -328,7 +332,7 @@
   where
     -- clamp
     l' = max 0 l
-    r' = min (Rwm.lengthRwm (rawWM wm)) r
+    r' = min (Rwm.lengthRwm (rawWm wm)) r
     rank_ = rankBetween wm l' r' minBound (y0 + 1)
 
 -- | \(O(\log |S|)\) Looks up the maximum \(y\) in \([l, r) \times (-\infty, y_0)\).
@@ -370,7 +374,7 @@
   where
     -- clamp
     l' = max 0 l
-    r' = min (Rwm.lengthRwm (rawWM wm)) r
+    r' = min (Rwm.lengthRwm (rawWm wm)) r
     rank_ = rankBetween wm l' r' minBound y0
 
 -- | \(O(\log |S|)\) Looks up the minimum \(y\) in \([l, r) \times (y_0, \infty)\).
@@ -396,7 +400,7 @@
 -- @since 1.1.0.0
 {-# INLINE assocsIn #-}
 assocsIn :: WaveletMatrix -> Int -> Int -> [(Int, Int)]
-assocsIn WaveletMatrix {..} l r = Rwm.assocsWith rawWM l r (xDictWM VG.!)
+assocsIn WaveletMatrix {..} l r = Rwm.assocsWith rawWm l r (yDictWm VG.!)
 
 -- | \(O(\min(|S|, L) \log |S|)\) Collects \((y, \mathrm{rank}(y))\) in range \([l, r)\) in
 -- descending order of \(y\). Note that it's only fast when the \(|S|\) is very small.
@@ -404,4 +408,4 @@
 -- @since 1.1.0.0
 {-# INLINE descAssocsIn #-}
 descAssocsIn :: WaveletMatrix -> Int -> Int -> [(Int, Int)]
-descAssocsIn WaveletMatrix {..} l r = Rwm.descAssocsInWith rawWM l r (xDictWM VG.!)
+descAssocsIn WaveletMatrix {..} l r = Rwm.descAssocsInWith rawWm l r (yDictWm VG.!)
diff --git a/src/AtCoder/Extra/WaveletMatrix/BitVector.hs b/src/AtCoder/Extra/WaveletMatrix/BitVector.hs
--- a/src/AtCoder/Extra/WaveletMatrix/BitVector.hs
+++ b/src/AtCoder/Extra/WaveletMatrix/BitVector.hs
@@ -27,7 +27,7 @@
   )
 where
 
-import AtCoder.Extra.Bisect (bisectL)
+import AtCoder.Extra.Bisect (maxRight)
 import Control.Monad.Primitive (PrimMonad (PrimState))
 import Data.Bit (Bit (..))
 import Data.Bits (popCount)
@@ -150,7 +150,8 @@
   Maybe Int
 selectKthIn0 bv l r k
   | k < 0 || nZeros <= k = Nothing
-  | otherwise = bisectL l r $ \i -> rank0 bv i - rankL0 < k + 1
+  -- note that `rank0` takes exclusive index
+  | otherwise = Just . maxRight l r $ \i -> rank0 bv (i + 1) - rankL0 < k + 1
   where
     nZeros = rank0 bv r - rankL0
     rankL0 = rank0 bv l
@@ -173,7 +174,8 @@
   Maybe Int
 selectKthIn1 bv l r k
   | k < 0 || nOnes <= k = Nothing
-  | otherwise = bisectL l r $ \i -> rank1 bv i - rankL1 < k + 1
+  -- note that `rank1` takes exclusive index
+  | otherwise = Just . maxRight l r $ \i -> rank1 bv (i + 1) - rankL1 < k + 1
   where
     nOnes = rank1 bv r - rankL1
     rankL1 = rank1 bv l
diff --git a/src/AtCoder/Extra/WaveletMatrix2d.hs b/src/AtCoder/Extra/WaveletMatrix2d.hs
--- a/src/AtCoder/Extra/WaveletMatrix2d.hs
+++ b/src/AtCoder/Extra/WaveletMatrix2d.hs
@@ -60,7 +60,7 @@
   )
 where
 
-import AtCoder.Extra.Bisect (bisectR, lowerBound)
+import AtCoder.Extra.Bisect (lowerBound)
 import AtCoder.Extra.WaveletMatrix.BitVector qualified as BV
 import AtCoder.Extra.WaveletMatrix.Raw qualified as Rwm
 import AtCoder.Internal.Assert qualified as ACIA
@@ -69,7 +69,6 @@
 import Control.Monad.ST (ST)
 import Data.Bit (Bit (..))
 import Data.Bits (Bits (testBit))
-import Data.Maybe (fromJust, fromMaybe)
 import Data.Vector qualified as V
 import Data.Vector.Algorithms.Intro qualified as VAI
 import Data.Vector.Generic qualified as VG
@@ -85,12 +84,12 @@
 
 -- | Segment Tree on Wavelet Matrix: points on a 2D plane and rectangle products.
 --
--- @since 1.1.0.0
+-- @since 1.3.0.0
 data WaveletMatrix2d s a = WaveletMatrix2d
   { -- | The wavelet matrix that represents points on a 2D plane.
     --
-    -- @since 1.1.0.0
-    rawWmWm2d :: !Rwm.RawWaveletMatrix,
+    -- @since 1.3.0.0
+    rawWm2d :: !Rwm.RawWaveletMatrix,
     -- | (x, y) index compression dictionary.
     --
     -- @since 1.1.0.0
@@ -115,7 +114,7 @@
 -- @since 1.1.0.0
 {-# INLINEABLE new #-}
 new ::
-  (PrimMonad m, Monoid a, VU.Unbox a) =>
+  (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) =>
   -- | Inverse operator of the monoid
   (a -> a) ->
   -- | Input points
@@ -130,8 +129,8 @@
   -- REMARK: Be sure to use `n + 1` because the product function cannot handle the case
   --         `yUpper` is `2^{height}`.
   let (!_, !ysInput) = VU.unzip xyDictWm2d
-  let rawWmWm2d = Rwm.build (n + 1) $ VU.map (fromJust . lowerBound yDictWm2d) ysInput
-  segTreesWm2d <- V.replicateM (Rwm.heightRwm rawWmWm2d) (ST.new n)
+  let rawWm2d = Rwm.build (n + 1) $ VU.map (lowerBound yDictWm2d) ysInput
+  segTreesWm2d <- V.replicateM (Rwm.heightRwm rawWm2d) (ST.new n)
   pure WaveletMatrix2d {..}
 
 -- | \(O(n \log n)\) Creates a `WaveletMatrix2d` with wavelet matrix with segment tree
@@ -140,7 +139,7 @@
 -- @since 1.1.0.0
 {-# INLINEABLE build #-}
 build ::
-  (PrimMonad m, Monoid a, VU.Unbox a) =>
+  (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) =>
   -- | Inverse operator of the monoid
   (a -> a) ->
   -- | Input points with initial values
@@ -161,7 +160,7 @@
 {-# INLINEABLE read #-}
 read :: (HasCallStack, PrimMonad m, VU.Unbox a, Monoid a) => WaveletMatrix2d (PrimState m) a -> (Int, Int) -> m a
 read WaveletMatrix2d {..} (!x, !y) = do
-  ST.read (V.head segTreesWm2d) . fromJust $ lowerBound xyDictWm2d (x, y)
+  ST.read (V.head segTreesWm2d) $ lowerBound xyDictWm2d (x, y)
 
 -- | \(O(\log^2 n)\) Writes the monoid value at \((x, y)\). Access to unknown points are undefined.
 --
@@ -169,19 +168,19 @@
 {-# INLINEABLE write #-}
 write :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => WaveletMatrix2d (PrimState m) a -> (Int, Int) -> a -> m ()
 write WaveletMatrix2d {..} (!x, !y) v = stToPrim $ do
-  let !i_ = fromJust $ lowerBound xyDictWm2d (x, y)
+  let !i_ = lowerBound xyDictWm2d (x, y)
   V.ifoldM'_
     ( \i iRow (!bits, !seg) -> do
         let !i0 = BV.rank0 bits i
         let !i'
               | unBit $ VG.unsafeIndex (BV.bitsBv bits) i =
-                  i + Rwm.nZerosRwm rawWmWm2d VG.! iRow - i0
+                  i + Rwm.nZerosRwm rawWm2d VG.! iRow - i0
               | otherwise = i0
         ST.write seg i' v
         pure i'
     )
     i_
-    $ V.zip (Rwm.bitsRwm rawWmWm2d) segTreesWm2d
+    $ V.zip (Rwm.bitsRwm rawWm2d) segTreesWm2d
 
 -- | \(O(\log^2 n)\) Modifies the monoid value at \((x, y)\). Access to unknown points are
 -- undefined.
@@ -190,19 +189,19 @@
 {-# INLINEABLE modify #-}
 modify :: (HasCallStack, PrimMonad m, Monoid a, VU.Unbox a) => WaveletMatrix2d (PrimState m) a -> (a -> a) -> (Int, Int) -> m ()
 modify WaveletMatrix2d {..} f (!x, !y) = stToPrim $ do
-  let !i_ = fromJust $ lowerBound xyDictWm2d (x, y)
+  let !i_ = lowerBound xyDictWm2d (x, y)
   V.ifoldM'_
     ( \i iRow (!bits, !seg) -> do
         let !i0 = BV.rank0 bits i
         let !i'
               | unBit $ VG.unsafeIndex (BV.bitsBv bits) i =
-                  i + Rwm.nZerosRwm rawWmWm2d VG.! iRow - i0
+                  i + Rwm.nZerosRwm rawWm2d VG.! iRow - i0
               | otherwise = i0
         ST.modify seg f i'
         pure i'
     )
     i_
-    $ V.zip (Rwm.bitsRwm rawWmWm2d) segTreesWm2d
+    $ V.zip (Rwm.bitsRwm rawWm2d) segTreesWm2d
 
 -- | \(O(\log^2 n)\) Returns monoid product \(\Pi_{p \in [x_1, x_2) \times [y_1, y_2)} a_p\).
 --
@@ -215,10 +214,10 @@
   where
     (!xDict, !_) = VU.unzip xyDictWm2d
     -- NOTE: clamping here!
-    xl' = fromMaybe 0 $ bisectR 0 (VG.length xDict) $ (< xl) . VG.unsafeIndex xDict
-    xr' = fromMaybe (VG.length xDict) $ bisectR 0 (VG.length xDict) $ (< xr) . VG.unsafeIndex xDict
-    yl' = fromMaybe 0 $ bisectR 0 (VG.length yDictWm2d) $ (< yl) . VG.unsafeIndex yDictWm2d
-    yr' = fromMaybe (VG.length yDictWm2d) $ bisectR 0 (VG.length yDictWm2d) $ (< yr) . VG.unsafeIndex yDictWm2d
+    xl' = lowerBound xDict xl
+    xr' = lowerBound xDict xr
+    yl' = lowerBound yDictWm2d yl
+    yr' = lowerBound yDictWm2d yr
     !_ = ACIA.checkInterval "AtCoder.Extra.WaveletMatrix.SegTree.prod (compressed x)" xl' xr' (VG.length xDict)
     !_ = ACIA.checkInterval "AtCoder.Extra.WaveletMatrix.SegTree.prod (compressed y)" yl' yr' (VG.length yDictWm2d)
 
@@ -236,10 +235,10 @@
   where
     (!xDict, !_) = VU.unzip xyDictWm2d
     -- NOTE: clamping here!
-    xl' = fromMaybe 0 $ bisectR 0 (VG.length xDict) $ (< xl) . VG.unsafeIndex xDict
-    xr' = fromMaybe (VG.length xDict) $ bisectR 0 (VG.length xDict) $ (< xr) . VG.unsafeIndex xDict
-    yl' = fromMaybe 0 $ bisectR 0 (VG.length yDictWm2d) $ (< yl) . VG.unsafeIndex yDictWm2d
-    yr' = fromMaybe (VG.length yDictWm2d) $ bisectR 0 (VG.length yDictWm2d) $ (< yr) . VG.unsafeIndex yDictWm2d
+    xl' = lowerBound xDict xl
+    xr' = lowerBound xDict xr
+    yl' = lowerBound yDictWm2d yl
+    yr' = lowerBound yDictWm2d yr
 
 -- | \(O(\log^2 n)\) Return the monoid product in \([-\infty, \infty) \times [-\infty, \infty)\).
 --
@@ -271,15 +270,15 @@
               !r0 = BV.rank0 bits r
           -- REMARK: The function cannot handle the case yUpper = N = 2^i. See the constructor for
           -- how it's handled and note that l_ and r_ are compressed indices.
-          if testBit yUpper (Rwm.heightRwm rawWmWm2d - 1 - iRow)
+          if testBit yUpper (Rwm.heightRwm rawWm2d - 1 - iRow)
             then do
               !acc' <- (acc <>) <$> ST.prod seg l0 r0
-              let !l' = l + Rwm.nZerosRwm rawWmWm2d VG.! iRow - l0
-              let !r' = r + Rwm.nZerosRwm rawWmWm2d VG.! iRow - r0
+              let !l' = l + Rwm.nZerosRwm rawWm2d VG.! iRow - l0
+              let !r' = r + Rwm.nZerosRwm rawWm2d VG.! iRow - r0
               pure (acc', l', r')
             else do
               pure (acc, l0, r0)
       )
       (mempty, l_, r_)
-      $ V.zip (Rwm.bitsRwm rawWmWm2d) segTreesWm2d
+      $ V.zip (Rwm.bitsRwm rawWm2d) segTreesWm2d
   pure res
diff --git a/src/AtCoder/Internal/Math.hs b/src/AtCoder/Internal/Math.hs
--- a/src/AtCoder/Internal/Math.hs
+++ b/src/AtCoder/Internal/Math.hs
@@ -107,7 +107,7 @@
   | otherwise =
       let d = innerD $ n - 1
           test a = inner d $ powMod a d n
-       in all test [2, 7, 61 :: Int]
+       in test 2 && test 7 && test 61
   where
     innerD d
       | even d = innerD $ d `div` 2
diff --git a/src/AtCoder/Internal/Queue.hs b/src/AtCoder/Internal/Queue.hs
--- a/src/AtCoder/Internal/Queue.hs
+++ b/src/AtCoder/Internal/Queue.hs
@@ -1,6 +1,6 @@
 {-# LANGUAGE RecordWildCards #-}
 
--- | Fixed-sized queue. Internally it has an \([l, r)\) pair of valid element bounds.
+-- | Fixed-sized double-ended queue. Internally it has an \([l, r)\) pair of valid element bounds.
 --
 -- ==== __Example__
 -- >>> import AtCoder.Internal.Queue qualified as Q
@@ -177,7 +177,7 @@
 import GHC.Stack (HasCallStack)
 import Prelude hiding (length, null)
 
--- | Fixed-sized queue. Internally it has an \([l, r)\) pair of valid element bounds.
+-- | Fixed-sized double-ended queue. Internally it has an \([l, r)\) pair of valid element bounds.
 --
 -- @since 1.0.0.0
 data Queue s a = Queue
@@ -193,7 +193,7 @@
 new :: (PrimMonad m, VU.Unbox a) => Int -> m (Queue (PrimState m) a)
 new n = stToPrim $ newST n
 
--- | \(O(n)\) Creates a `Queue` with capacity \(2n + 1\) where the internal front/back position is
+-- | \(O(n)\) Creates a `Queue` with capacity \(2n + 1\), where the internal front/back position is
 -- initialzed at \(n\).
 --
 -- @since 1.2.4.0
diff --git a/src/AtCoder/LazySegTree.hs b/src/AtCoder/LazySegTree.hs
--- a/src/AtCoder/LazySegTree.hs
+++ b/src/AtCoder/LazySegTree.hs
@@ -656,6 +656,8 @@
             else inner2 l' sm
       | otherwise = pure $ l - sizeLst
 
+-- TODO: add minLeft
+
 -- | \(O(n)\) Yields an immutable copy of the mutable vector.
 --
 -- @since 1.0.0.0
diff --git a/src/AtCoder/SegTree.hs b/src/AtCoder/SegTree.hs
--- a/src/AtCoder/SegTree.hs
+++ b/src/AtCoder/SegTree.hs
@@ -449,6 +449,8 @@
             else inner2 l' sm
       | otherwise = pure $ l - sizeSt
 
+-- TODO: add `minLeft`
+
 -- | \(O(n)\) Yields an immutable copy of the mutable vector.
 --
 -- @since 1.0.0.0
diff --git a/test/Tests/Extra/Bisect.hs b/test/Tests/Extra/Bisect.hs
--- a/test/Tests/Extra/Bisect.hs
+++ b/test/Tests/Extra/Bisect.hs
@@ -8,35 +8,25 @@
 import Test.Tasty.QuickCheck as QC
 
 -- | Takes half-open interval [l, r).
-naivePartition :: Int -> Int -> (Int -> Bool) -> VU.Vector Int -> (Maybe Int, Maybe Int)
-naivePartition l r p xs
-  | l >= r = (Nothing, Nothing)
-  | otherwise = case (VU.null ls, VU.null rs) of
-      (True, True) -> error "unreachable"
-      (False, True) -> (Just l', Nothing)
-      (True, False) -> (Nothing, Just r')
-      _ -> (Just l', Just r')
-  where
-    xs' = VU.take (r - l) . VU.drop l $ xs
-    (!ls, !rs) = VU.partition p xs'
-    l' = l + VU.length ls - 1
-    r' = l' + 1
+naiveMaxRightIn :: Int -> Int -> (Int -> Bool) -> VU.Vector Int -> Int
+naiveMaxRightIn l r p xs =
+  (+ l)
+    . VU.length
+    . VU.takeWhile p
+    . VU.take (r - l)
+    $ VU.drop l xs
 
-naiveLowerBound :: VU.Vector Int -> Int -> Maybe Int
+naiveLowerBound :: VU.Vector Int -> Int -> Int
 naiveLowerBound xs = naiveLowerBoundIn 0 (VU.length xs) xs
 
-naiveLowerBoundIn :: Int -> Int -> VU.Vector Int -> Int -> Maybe Int
-naiveLowerBoundIn l r xs target = case naivePartition l r (< target) xs of
-  (!_, Just i) -> Just i
-  _ -> Nothing
+naiveLowerBoundIn :: Int -> Int -> VU.Vector Int -> Int -> Int
+naiveLowerBoundIn l r xs target = naiveMaxRightIn l r (< target) xs
 
-naiveUpperBound :: VU.Vector Int -> Int -> Maybe Int
+naiveUpperBound :: VU.Vector Int -> Int -> Int
 naiveUpperBound xs = naiveUpperBoundIn 0 (VU.length xs) xs
 
-naiveUpperBoundIn :: Int -> Int -> VU.Vector Int -> Int -> Maybe Int
-naiveUpperBoundIn l r xs target = case naivePartition l r (<= target) xs of
-  (!_, Just i) -> Just i
-  _ -> Nothing
+naiveUpperBoundIn :: Int -> Int -> VU.Vector Int -> Int -> Int
+naiveUpperBoundIn l r xs target = naiveMaxRightIn l r (<= target) xs
 
 boundsQueryGen :: Gen (Int, Int, VU.Vector Int)
 boundsQueryGen = do
@@ -45,8 +35,8 @@
   xs <- VU.fromList . L.sort <$> QC.vectorOf n (QC.chooseInt (-20, 20))
   pure (n, p, xs)
 
-bisectQueryGen :: Gen (Int, Int, VU.Vector Int, [(Int, Int)])
-bisectQueryGen = do
+maxRightQueryGen :: Gen (Int, Int, VU.Vector Int, [(Int, Int)])
+maxRightQueryGen = do
   n <- QC.chooseInt (1, 100)
   p <- QC.chooseInt (-25, 25)
   xs <- VU.fromList . L.sort <$> QC.vectorOf n (QC.chooseInt (-20, 20))
@@ -56,11 +46,11 @@
 prop_lowerBound :: TestTree
 prop_lowerBound = QC.testProperty "lowerBound" $ do
   (!_, !target, !xs) <- boundsQueryGen
-  pure $ naiveLowerBound xs target QC.=== lowerBound xs target
+  pure . QC.counterexample (show (target, xs)) $ naiveLowerBound xs target QC.=== lowerBound xs target
 
 prop_lowerBoundIn :: TestTree
 prop_lowerBoundIn = QC.testProperty "lowerBoundIn" $ do
-  (!_, !target, !xs, !lrs) <- bisectQueryGen
+  (!_, !target, !xs, !lrs) <- maxRightQueryGen
   pure . QC.conjoin $
     map
       ( \(!l, !r) ->
@@ -75,7 +65,7 @@
 
 prop_upperBoundIn :: TestTree
 prop_upperBoundIn = QC.testProperty "upperBoundIn" $ do
-  (!_, !target, !xs, !lrs) <- bisectQueryGen
+  (!_, !target, !xs, !lrs) <- maxRightQueryGen
   pure . QC.conjoin $
     map
       ( \(!l, !r) ->
@@ -83,23 +73,40 @@
       )
       lrs
 
-prop_bisectL :: TestTree
-prop_bisectL = QC.testProperty "bisectL" $ do
-  (!_, !boundary, !xs, !lrs) <- bisectQueryGen
+prop_maxRight :: TestTree
+prop_maxRight = QC.testProperty "maxRight" $ do
+  (!_, !boundary, !xs, !lrs) <- maxRightQueryGen
   pure . QC.conjoin $
     map
       ( \(!l, !r) ->
-          fst (naivePartition l r (<= boundary) xs) == bisectL l r (\i -> xs VG.! i <= boundary)
+          naiveMaxRightIn l r (<= boundary) xs == maxRight l r (\i -> xs VG.! i <= boundary)
       )
       lrs
 
-prop_bisectR :: TestTree
-prop_bisectR = QC.testProperty "bisectR" $ do
-  (!_, !boundary, !xs, !lrs) <- bisectQueryGen
+minLeftQueryGen :: Gen (Int, Int, VU.Vector Int, [(Int, Int)])
+minLeftQueryGen = do
+  n <- QC.chooseInt (1, 100)
+  p <- QC.chooseInt (-25, 25)
+  xs <- VU.fromList . L.sort <$> QC.vectorOf n (QC.chooseInt (-20, 20))
+  let lrs = [(l, r) | l <- [0 .. n], r <- [l .. n]]
+  pure (n, p, xs, lrs)
+
+prop_minLeft :: TestTree
+prop_minLeft = QC.testProperty "minLeft" $ do
+  (!_, !boundary, !xs, !lrs) <- minLeftQueryGen
   pure . QC.conjoin $
     map
       ( \(!l, !r) ->
-          snd (naivePartition l r (<= boundary) xs) == bisectR l r (\i -> xs VG.! i <= boundary)
+          let expected =
+                (r -)
+                  . VU.length
+                  . VU.takeWhile (>= boundary)
+                  . VU.reverse
+                  . VU.take (r - l)
+                  $ VU.drop l xs
+              res = minLeft l r (\i -> xs VG.! i >= boundary)
+           in QC.counterexample (show ((l, r), boundary, xs)) $
+                expected QC.=== res
       )
       lrs
 
@@ -109,6 +116,6 @@
     prop_upperBound,
     prop_lowerBoundIn,
     prop_upperBoundIn,
-    prop_bisectL,
-    prop_bisectR
+    prop_maxRight,
+    prop_minLeft
   ]
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
@@ -11,7 +11,6 @@
 import Test.Tasty.QuickCheck as QC
 import Tests.Util
 import Test.Tasty.HUnit
-import Debug.Trace
 
 data Init = Init
   { n :: {-# UNPACK #-} !Int,
@@ -39,7 +38,6 @@
       x <- QC.chooseInt rng
       y <- QC.chooseInt rng
       pure (x, y)
-    let !_ = traceShow refVec ()
     let kt = Kt.build2 refVec
     pure Init {..}
 
diff --git a/test/Tests/Extra/Math.hs b/test/Tests/Extra/Math.hs
--- a/test/Tests/Extra/Math.hs
+++ b/test/Tests/Extra/Math.hs
@@ -1,14 +1,12 @@
 module Tests.Extra.Math (tests) where
 
 import AtCoder.Extra.Math qualified as ACEM
-import Data.Foldable (for_)
 import Data.List qualified as L
 import Data.Proxy (Proxy (..))
 import Data.Semigroup (Max (..), Min (..), Sum (..), mtimesDefault, stimes)
 import Data.Vector.Unboxed qualified as VU
 import Test.QuickCheck.Property qualified as QC
 import Test.Tasty
-import Test.Tasty.HUnit
 import Test.Tasty.QuickCheck qualified as QC
 import Tests.Util (myForAllShrink)
 
diff --git a/test/Tests/Extra/ModInt64.hs b/test/Tests/Extra/ModInt64.hs
--- a/test/Tests/Extra/ModInt64.hs
+++ b/test/Tests/Extra/ModInt64.hs
@@ -117,6 +117,14 @@
   where
     !m = fromIntegral $ natVal' (proxy# @a)
 
+prop_quotRem2 :: forall a. (KnownNat a) => Proxy# a -> Int -> QC.NonZero Int -> QC.Property
+prop_quotRem2 _ x (QC.NonZero y) =
+  let !a = M.new @a x
+      !d = M.new @a y
+   in d /= 0 QC.==>
+        let (!q, !r) = a `quotRem` d
+         in (d * q + r) QC.=== a
+
 prop_eq :: forall a. (KnownNat a) => Proxy# a -> Word64 -> Word64 -> QC.Property
 prop_eq _ x y = lhs QC.=== rhs
   where
@@ -160,6 +168,14 @@
         QC.testProperty "4" (prop_quotRem (proxy# @M4))
         -- 64 bit
         -- QC.testProperty "5" (prop_quotRem (proxy# @M5))
+      ],
+    testGroup
+      "quotRem2"
+      [ QC.testProperty "1" (prop_quotRem2 (proxy# @M1)),
+        QC.testProperty "2" (prop_quotRem2 (proxy# @M2)),
+        QC.testProperty "3" (prop_quotRem2 (proxy# @M3)),
+        QC.testProperty "4" (prop_quotRem2 (proxy# @M4)),
+        QC.testProperty "5" (prop_quotRem2 (proxy# @M5))
       ],
     testGroup
       "eq"
diff --git a/test/Tests/Extra/Seq.hs b/test/Tests/Extra/Seq.hs
--- a/test/Tests/Extra/Seq.hs
+++ b/test/Tests/Extra/Seq.hs
@@ -229,7 +229,7 @@
 handleAcl seq handle q = case q of
   Reset -> do
     Seq.reset seq
-    Seq.invalidateHandle handle
+    P.invalidateHandle handle
     pure None
   Read k -> S <$> Seq.read seq handle k
   ReadMaybe k -> MS <$> Seq.readMaybe seq handle k
diff --git a/test/Tests/Extra/WaveletMatrix/Raw.hs b/test/Tests/Extra/WaveletMatrix/Raw.hs
--- a/test/Tests/Extra/WaveletMatrix/Raw.hs
+++ b/test/Tests/Extra/WaveletMatrix/Raw.hs
@@ -6,7 +6,6 @@
 import AtCoder.Extra.WaveletMatrix.Raw qualified as WM
 import Control.Exception (evaluate)
 import Data.IntMap qualified as IM
-import Data.Maybe (fromJust)
 import Data.Ord (comparing)
 import Data.Vector.Algorithms.Intro qualified as VAI
 import Data.Vector.Unboxed qualified as VU
@@ -21,7 +20,7 @@
 compress :: VU.Vector Int -> VU.Vector Int
 compress xs =
   let dict = VU.uniq $ VU.modify VAI.sort xs
-   in VU.map (fromJust . lowerBound dict) xs
+   in VU.map (lowerBound dict) xs
 
 data Init = Init
   { capacity :: {-# UNPACK #-} !Int,