diff --git a/Data/Mod.hs b/Data/Mod.hs
--- a/Data/Mod.hs
+++ b/Data/Mod.hs
@@ -1,12 +1,12 @@
 -- |
 -- Module:      Data.Mod
--- Copyright:   (c) 2017-2020 Andrew Lelechenko
+-- Copyright:   (c) 2017-2022 Andrew Lelechenko
 -- Licence:     MIT
 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
 --
 -- <https://en.wikipedia.org/wiki/Modular_arithmetic Modular arithmetic>,
 -- promoting moduli to the type level, with an emphasis on performance.
--- Originally part of <https://hackage.haskell.org/package/arithmoi arithmoi> package.
+-- Originally part of the <https://hackage.haskell.org/package/arithmoi arithmoi> package.
 --
 -- This module supports moduli of arbitrary size.
 -- Use "Data.Mod.Word" to achieve better performance,
@@ -34,10 +34,11 @@
 import Control.DeepSeq
 import Control.Monad
 import Data.Bits
+import Data.Mod.Compat (timesWord2#, remWord2#)
+import Data.Ratio
 import Data.Word (Word8)
 #ifdef MIN_VERSION_semirings
 import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)
-import Data.Ratio
 import Data.Semiring (Semiring(..), Ring(..))
 #endif
 #ifdef MIN_VERSION_vector
@@ -49,51 +50,58 @@
 import qualified Data.Vector.Unboxed         as U
 import qualified Data.Vector.Primitive       as P
 import Foreign (copyBytes)
-import GHC.IO.Unsafe (unsafeDupablePerformIO)
 #endif
 import Foreign.Storable (Storable(..))
-import GHC.Exts
+import GHC.Exts hiding (timesWord2#, quotRemWord2#)
 import GHC.Generics
-import GHC.Integer.GMP.Internals
+import GHC.IO (IO(..))
 import GHC.Natural (Natural(..), powModNatural)
+import GHC.Num.BigNat
+import GHC.Num.Integer
 import GHC.TypeNats (Nat, KnownNat, natVal, natVal')
+import Text.Read (Read(readPrec))
 
 -- | This data type represents
 -- <https://en.wikipedia.org/wiki/Modular_arithmetic#Integers_modulo_n integers modulo m>,
 -- equipped with useful instances.
 --
 -- For example, 3 :: 'Mod' 10 stands for the class of integers
--- congruent to \( 3 \bmod 10 \colon \ldots {}−17, −7, 3, 13, 23 \ldots \)
+-- congruent to \( 3 \bmod 10 \colon \ldots {−17}, −7, 3, 13, 23 \ldots \)
 --
 -- >>> :set -XDataKinds
 -- >>> 3 + 8 :: Mod 10 -- 3 + 8 = 11 ≡ 1 (mod 10)
--- (1 `modulo` 10)
+-- 1
 --
--- __Warning:__ division by residue, which is not
--- <https://en.wikipedia.org/wiki/Coprime_integers coprime>
--- with the modulo, throws 'DivideByZero'.
--- Consider using 'invertMod' for non-prime moduli.
+-- __Note:__ 'Mod' 0 has no inhabitants, eventhough \( \mathbb{Z}/0\mathbb{Z} \) is technically isomorphic to \( \mathbb{Z} \).
 newtype Mod (m :: Nat) = Mod
   { unMod :: Natural
   -- ^ The canonical representative of the residue class,
-  -- always between 0 and \( m - 1 \) inclusively.
+  -- always between 0 and \( m - 1 \) (inclusively).
   --
   -- >>> :set -XDataKinds
   -- >>> -1 :: Mod 10
-  -- (9 `modulo` 10)
+  -- 9
   }
   deriving (Eq, Ord, Generic)
 
 instance NFData (Mod m)
 
-instance KnownNat m => Show (Mod m) where
-  show m = "(" ++ show (unMod m) ++ " `modulo` " ++ show (natVal m) ++ ")"
+instance Show (Mod m) where
+  show (Mod x) = show x
 
+-- | Wrapping behaviour, similar to
+-- the existing @instance@ 'Read' 'Int'.
+instance KnownNat m => Read (Mod m) where
+  readPrec = fromInteger <$> readPrec
+
+instance KnownNat m => Real (Mod m) where
+  toRational (Mod x) = toRational x
+
 instance KnownNat m => Enum (Mod m) where
   succ x = if x == maxBound then throw Overflow  else coerce (succ @Natural) x
   pred x = if x == minBound then throw Underflow else coerce (pred @Natural) x
 
-  toEnum   = (fromIntegral :: Int -> Mod m)
+  toEnum   = fromIntegral :: Int -> Mod m
   fromEnum = (fromIntegral :: Natural -> Int) . unMod
 
   enumFrom x       = enumFromTo x maxBound
@@ -103,25 +111,22 @@
   enumFromThenTo = coerce (enumFromThenTo @Natural)
 
 instance KnownNat m => Bounded (Mod m) where
-  minBound = Mod 0
-  maxBound = let mx = Mod (natVal mx - 1) in mx
+  minBound = mx
+    where
+      mx = if natVal mx > 0 then Mod 0 else throw DivideByZero
+  maxBound = mx
+    where
+      mx = if m > 0 then Mod (m - 1) else throw DivideByZero
+      m = natVal mx
 
-bigNatToNat :: BigNat -> Natural
+bigNatToNat :: BigNat# -> Natural
 bigNatToNat r# =
-  if isTrue# (sizeofBigNat# r# <=# 1#) then NatS# (bigNatToWord r#) else NatJ# r#
-
-subIfGe :: BigNat -> BigNat -> Natural
-subIfGe z# m# = case z# `compareBigNat` m# of
-  LT -> NatJ# z#
-  EQ -> NatS# 0##
-  GT -> bigNatToNat $ z# `minusBigNat` m#
+  if isTrue# (bigNatSize# r# <=# 1#) then NatS# (bigNatToWord# r#) else NatJ# (BN# r#)
 
-#if !MIN_VERSION_base(4,12,0)
-addWordC# :: Word# -> Word# -> (# Word#, Int# #)
-addWordC# x# y# = (# z#, word2Int# c# #)
-  where
-    !(# c#, z# #) = x# `plusWord2#` y#
-#endif
+subIfGe :: BigNat# -> BigNat# -> Natural
+subIfGe z# m# = case z# `bigNatSub` m# of
+  (# (# #) | #) -> NatJ# (BN# z#)
+  (# | zm# #)   -> bigNatToNat zm#
 
 addMod :: Natural -> Natural -> Natural -> Natural
 addMod (NatS# m#) (NatS# x#) (NatS# y#) =
@@ -129,13 +134,13 @@
   where
     !(# z#, c# #) = x# `addWordC#` y#
 addMod NatS#{} _ _ = brokenInvariant
-addMod (NatJ# m#) (NatS# x#) (NatS# y#) =
-  if isTrue# c# then subIfGe (wordToBigNat2 1## z#) m# else NatS# z#
+addMod (NatJ# (BN# m#)) (NatS# x#) (NatS# y#) =
+  if isTrue# c# then subIfGe (bigNatFromWord2# 1## z#) m# else NatS# z#
   where
     !(# z#, c# #) = x# `addWordC#` y#
-addMod (NatJ# m#) (NatS# x#) (NatJ# y#) = subIfGe (y# `plusBigNatWord` x#) m#
-addMod (NatJ# m#) (NatJ# x#) (NatS# y#) = subIfGe (x# `plusBigNatWord` y#) m#
-addMod (NatJ# m#) (NatJ# x#) (NatJ# y#) = subIfGe (x# `plusBigNat`     y#) m#
+addMod (NatJ# (BN# m#)) (NatS# x#) (NatJ# (BN# y#)) = subIfGe (y# `bigNatAddWord#` x#) m#
+addMod (NatJ# (BN# m#)) (NatJ# (BN# x#)) (NatS# y#) = subIfGe (x# `bigNatAddWord#` y#) m#
+addMod (NatJ# (BN# m#)) (NatJ# (BN# x#)) (NatJ# (BN# y#)) = subIfGe (x# `bigNatAdd` y#) m#
 
 subMod :: Natural -> Natural -> Natural -> Natural
 subMod (NatS# m#) (NatS# x#) (NatS# y#) =
@@ -143,45 +148,50 @@
   where
     z# = x# `minusWord#` y#
 subMod NatS#{} _ _ = brokenInvariant
-subMod (NatJ# m#) (NatS# x#) (NatS# y#) =
+subMod (NatJ# (BN# m#)) (NatS# x#) (NatS# y#) =
   if isTrue# (x# `geWord#` y#)
     then NatS# (x# `minusWord#` y#)
-    else bigNatToNat $ m# `minusBigNatWord` (y# `minusWord#` x#)
-subMod (NatJ# m#) (NatS# x#) (NatJ# y#) =
-  bigNatToNat $ (m# `minusBigNat` y#) `plusBigNatWord` x#
-subMod NatJ#{} (NatJ# x#) (NatS# y#) =
-  bigNatToNat $ x# `minusBigNatWord` y#
-subMod (NatJ# m#) (NatJ# x#) (NatJ# y#) = case x# `compareBigNat` y# of
-  LT -> bigNatToNat $ (m# `minusBigNat` y#) `plusBigNat` x#
-  EQ -> NatS# 0##
-  GT -> bigNatToNat $ x# `minusBigNat` y#
+    else bigNatToNat (m# `bigNatSubWordUnsafe#` (y# `minusWord#` x#))
+subMod (NatJ# (BN# m#)) (NatS# x#) (NatJ# (BN# y#)) =
+  bigNatToNat (m# `bigNatSubUnsafe` y# `bigNatAddWord#` x#)
+subMod NatJ#{} (NatJ# (BN# x#)) (NatS# y#) =
+  bigNatToNat (x# `bigNatSubWordUnsafe#` y#)
+subMod (NatJ# (BN# m#)) (NatJ# (BN# x#)) (NatJ# (BN# y#)) =
+  case x# `bigNatSub` y# of
+    (# (# #) | #) -> bigNatToNat (m# `bigNatSubUnsafe` y# `bigNatAdd` x#)
+    (# | xy# #) -> bigNatToNat xy#
 
 negateMod :: Natural -> Natural -> Natural
 negateMod _ (NatS# 0##) = NatS# 0##
 negateMod (NatS# m#) (NatS# x#) = NatS# (m# `minusWord#` x#)
 negateMod NatS#{} _ = brokenInvariant
-negateMod (NatJ# m#) (NatS# x#) = bigNatToNat $ m# `minusBigNatWord` x#
-negateMod (NatJ# m#) (NatJ# x#) = bigNatToNat $ m# `minusBigNat`     x#
+negateMod (NatJ# (BN# m#)) (NatS# x#) = bigNatToNat (m# `bigNatSubWordUnsafe#` x#)
+negateMod (NatJ# (BN# m#)) (NatJ# (BN# x#)) = bigNatToNat (m# `bigNatSubUnsafe` x#)
 
+halfWord :: Word
+halfWord = 1 `shiftL` (finiteBitSize (0 :: Word) `shiftR` 1)
+
 mulMod :: Natural -> Natural -> Natural -> Natural
-mulMod (NatS# m#) (NatS# x#) (NatS# y#) = NatS# r#
+mulMod (NatS# m#) (NatS# x#) (NatS# y#)
+  | W# m# <= halfWord = NatS# (timesWord# x# y# `remWord#` m#)
+  | otherwise = NatS# r#
   where
-    !(# z1#, z2# #) = timesWord2# x# y#
-    !(# _, r# #) = quotRemWord2# z1# z2# m#
+    !(# hi#, lo# #) = timesWord2# x# y#
+    !r# = remWord2# lo# hi# m#
 mulMod NatS#{} _ _ = brokenInvariant
-mulMod (NatJ# m#) (NatS# x#) (NatS# y#) =
-  bigNatToNat $ wordToBigNat2 z1# z2# `remBigNat` m#
+mulMod (NatJ# (BN# m#)) (NatS# x#) (NatS# y#) =
+  bigNatToNat (bigNatFromWord2# z1# z2# `bigNatRem` m#)
   where
     !(# z1#, z2# #) = timesWord2# x# y#
-mulMod (NatJ# m#) (NatS# x#) (NatJ# y#) =
-  bigNatToNat $ (y# `timesBigNatWord` x#) `remBigNat` m#
-mulMod (NatJ# m#) (NatJ# x#) (NatS# y#) =
-  bigNatToNat $ (x# `timesBigNatWord` y#) `remBigNat` m#
-mulMod (NatJ# m#) (NatJ# x#) (NatJ# y#) =
-  bigNatToNat $ (x# `timesBigNat` y#) `remBigNat` m#
+mulMod (NatJ# (BN# m#)) (NatS# x#) (NatJ# (BN# y#)) =
+  bigNatToNat ((y# `bigNatMulWord#` x#) `bigNatRem` m#)
+mulMod (NatJ# (BN# m#)) (NatJ# (BN# x#)) (NatS# y#) =
+  bigNatToNat ((x# `bigNatMulWord#` y#) `bigNatRem` m#)
+mulMod (NatJ# (BN# m#)) (NatJ# (BN# x#)) (NatJ# (BN# y#)) =
+  bigNatToNat ((x# `bigNatMul` y#) `bigNatRem` m#)
 
 brokenInvariant :: a
-brokenInvariant = error "argument is larger than modulo"
+brokenInvariant = error "argument is larger than modulus"
 
 instance KnownNat m => Num (Mod m) where
   mx@(Mod !x) + (Mod !y) = Mod $ addMod (natVal mx) x y
@@ -210,11 +220,17 @@
   {-# INLINE plus #-}
   times = (*)
   {-# INLINE times #-}
-  zero  = Mod 0
+  zero  = mx
+    where
+      mx = if natVal mx > 0 then Mod 0 else throw DivideByZero
   {-# INLINE zero #-}
   one   = mx
     where
-      mx = if natVal mx > 1 then Mod 1 else Mod 0
+      mx = case m `compare` 1 of
+        LT -> throw DivideByZero
+        EQ -> Mod 0
+        GT -> Mod 1
+      m = natVal mx
   {-# INLINE one #-}
   fromNatural x = mx
     where
@@ -225,7 +241,87 @@
   negate = Prelude.negate
   {-# INLINE negate #-}
 
--- | See the warning about division above.
+-- | 'Mod' @m@ is not even an
+-- <https://en.wikipedia.org/wiki/Integral_domain integral domain> for
+-- <https://en.wikipedia.org/wiki/Composite_number composite> @m@,
+-- much less a <https://en.wikipedia.org/wiki/GCD_domain GCD domain>.
+-- However, 'Data.Euclidean.gcd' and 'Data.Euclidean.lcm' are still meaningful
+-- even for composite @m@, corresponding to a sum and an intersection of
+-- <https://en.wikipedia.org/wiki/Ideal_(ring_theory) ideals>.
+--
+-- The instance is lawful only for
+-- <https://en.wikipedia.org/wiki/Prime_number prime> @m@, otherwise
+-- @'Data.Euclidean.divide' x y@ tries to return any @Just z@ such that @x == y * z@.
+--
+instance KnownNat m => GcdDomain (Mod m) where
+  divide (Mod 0) _ = Just (Mod 0)
+  divide _ (Mod 0) = Nothing
+  divide mx@(Mod x) (Mod y) = case mry of
+    Just ry -> if xr == 0 then Just (Mod xq * Mod ry) else Nothing
+    Nothing -> Nothing
+    where
+      m = natVal mx
+      gmy = Prelude.gcd m y
+      (xq, xr) = Prelude.quotRem x gmy
+      mry = invertModInternal (y `Prelude.quot` gmy)  (m `Prelude.quot` gmy)
+
+  gcd (Mod x) (Mod y) = g
+    where
+      m = natVal g
+      g = Mod $ if m > 1 then Prelude.gcd (Prelude.gcd m x) y else 0
+  lcm (Mod x) (Mod y) = l
+    where
+      m = natVal l
+      l = Mod $ if m > 1 then Prelude.lcm (Prelude.gcd m x) (Prelude.gcd m y) else 0
+  coprime x y = Data.Euclidean.gcd x y == one
+
+-- | 'Mod' @m@ is not even an
+-- <https://en.wikipedia.org/wiki/Integral_domain integral domain> for
+-- <https://en.wikipedia.org/wiki/Composite_number composite> @m@,
+-- much less a <https://en.wikipedia.org/wiki/Euclidean_domain Euclidean domain>.
+--
+-- The instance is lawful only for
+-- <https://en.wikipedia.org/wiki/Prime_number prime> @m@, otherwise
+-- we try to do our best:
+-- @'Data.Euclidean.quot' x y@ returns any @z@ such that @x == y * z@,
+-- 'Data.Euclidean.rem' is not always 0, and both can throw 'DivideByZero'.
+--
+instance KnownNat m => Euclidean (Mod m) where
+  degree = unMod
+  {-# INLINABLE degree #-}
+
+  quotRem (Mod 0) _ = (Mod 0, Mod 0)
+  quotRem _ (Mod 0) = throw DivideByZero
+  quotRem mx@(Mod x) (Mod y) = case mry of
+    Just ry -> (Mod xq * Mod ry, Mod xr)
+    Nothing -> throw DivideByZero
+    where
+      m = natVal mx
+      gmy = Prelude.gcd m y
+      (xq, xr) = Prelude.quotRem x gmy
+      mry = invertModInternal (y `Prelude.quot` gmy)  (m `Prelude.quot` gmy)
+
+-- | 'Mod' @m@ is not even an
+-- <https://en.wikipedia.org/wiki/Integral_domain integral domain> for
+-- <https://en.wikipedia.org/wiki/Composite_number composite> @m@,
+-- much less a <https://en.wikipedia.org/wiki/Field_(mathematics) field>.
+--
+-- The instance is lawful only for
+-- <https://en.wikipedia.org/wiki/Prime_number prime> @m@, otherwise
+-- division by a residue, which is not
+-- <https://en.wikipedia.org/wiki/Coprime_integers coprime>
+-- with the modulus, throws 'DivideByZero'.
+-- Consider using 'invertMod' for non-prime moduli.
+--
+instance KnownNat m => Field (Mod m)
+
+#endif
+
+-- | Division by a residue, which is not
+-- <https://en.wikipedia.org/wiki/Coprime_integers coprime>
+-- with the modulus, throws 'DivideByZero'.
+-- Consider using 'invertMod' for non-prime moduli.
+--
 instance KnownNat m => Fractional (Mod m) where
   fromRational r = case denominator r of
     1   -> num
@@ -238,55 +334,38 @@
     Just y  -> y
   {-# INLINE recip #-}
 
--- | See the warning about division above.
-instance KnownNat m => GcdDomain (Mod m) where
-  divide x y = Just (x / y)
-  gcd        = const $ const 1
-  lcm        = const $ const 1
-  coprime    = const $ const True
-
--- | See the warning about division above.
-instance KnownNat m => Euclidean (Mod m) where
-  degree      = const 0
-  quotRem x y = (x / y, 0)
-  quot        = (/)
-  rem         = const $ const 0
-
--- | See the warning about division above.
-instance KnownNat m => Field (Mod m)
-
-#endif
-
 -- | If an argument is
 -- <https://en.wikipedia.org/wiki/Coprime_integers coprime>
--- with the modulo, return its modular inverse.
+-- with the modulus, return its modular inverse.
 -- Otherwise return 'Nothing'.
 --
 -- >>> :set -XDataKinds
 -- >>> invertMod 3 :: Mod 10 -- 3 * 7 = 21 ≡ 1 (mod 10)
--- Just (7 `modulo` 10)
+-- Just 7
 -- >>> invertMod 4 :: Mod 10 -- 4 and 10 are not coprime
 -- Nothing
 invertMod :: KnownNat m => Mod m -> Maybe (Mod m)
-invertMod mx
-  = if y <= 0
-    then Nothing
-    else Just $ Mod $ fromInteger y
-  where
-    y = recipModInteger (toInteger (unMod mx)) (toInteger (natVal mx))
+invertMod x = Mod <$> invertModInternal (unMod x) (natVal x)
 {-# INLINABLE invertMod #-}
 
+invertModInternal
+  :: Natural -- Value
+  -> Natural -- Modulo
+  -> Maybe Natural
+invertModInternal x m = case integerRecipMod# (toInteger x) m of
+  (# | () #) -> Nothing
+  (# y | #)  -> Just y
+{-# INLINABLE invertModInternal #-}
+
 -- | Drop-in replacement for 'Prelude.^' with much better performance.
 -- Negative powers are allowed, but may throw 'DivideByZero', if an argument
--- is not <https://en.wikipedia.org/wiki/Coprime_integers coprime> with the modulo.
---
--- Building with @-O@ triggers a rewrite rule 'Prelude.^' = '^%'.
+-- is not <https://en.wikipedia.org/wiki/Coprime_integers coprime> with the modulus.
 --
 -- >>> :set -XDataKinds
 -- >>> 3 ^% 4 :: Mod 10    -- 3 ^ 4 = 81 ≡ 1 (mod 10)
--- (1 `modulo` 10)
+-- 1
 -- >>> 3 ^% (-1) :: Mod 10 -- 3 * 7 = 21 ≡ 1 (mod 10)
--- (7 `modulo` 10)
+-- 7
 -- >>> 4 ^% (-1) :: Mod 10 -- 4 and 10 are not coprime
 -- (*** Exception: divide by zero
 (^%) :: (KnownNat m, Integral a) => Mod m -> a -> Mod m
@@ -311,8 +390,6 @@
   KnownNat m => Mod m -> Word    -> Mod m #-}
 
 {-# RULES
-"powMod"               forall (x :: KnownNat m => Mod m) p. x ^ p = x ^% p
-
 "powMod/2/Integer"     forall x. x ^% (2 :: Integer) = let u = x in u*u
 "powMod/3/Integer"     forall x. x ^% (3 :: Integer) = let u = x in u*u*u
 "powMod/2/Int"         forall x. x ^% (2 :: Int)     = let u = x in u*u
@@ -331,10 +408,12 @@
   64 -> 3 -- 2^3 bytes in word
   _  -> error "lgWordSize: unknown architecture"
 
+-- | No validation checks are performed;
+-- reading untrusted data may corrupt internal invariants.
 instance KnownNat m => Storable (Mod m) where
   sizeOf _ = case natVal' (proxy# :: Proxy# m) of
     NatS#{}  -> sizeOf (0 :: Word)
-    NatJ# m# -> I# (sizeofBigNat# m#) `shiftL` lgWordSize
+    NatJ# (BN# m#) -> I# (bigNatSize# m#) `shiftL` lgWordSize
   {-# INLINE sizeOf #-}
 
   alignment _ = alignment (0 :: Word)
@@ -344,10 +423,10 @@
     NatS#{} -> do
       W# w# <- peek (Ptr addr#)
       pure . Mod $! NatS# w#
-    NatJ# m# -> do
+    NatJ# (BN# m#) -> do
       let !(I# lgWordSize#) = lgWordSize
-          sz# = sizeofBigNat# m# `iShiftL#` lgWordSize#
-      bn <- importBigNatFromAddr addr# (int2Word# sz#) 0#
+          sz# = bigNatSize# m# `iShiftL#` lgWordSize#
+      BN# bn <- IO (\token -> case bigNatFromAddrLE# (int2Word# sz#) addr# token of (# newToken, bn# #) -> (# newToken, BN# bn# #))
       pure . Mod $! bigNatToNat bn
   {-# INLINE peek #-}
 
@@ -355,21 +434,23 @@
     NatS#{} -> case x of
       NatS# x# -> poke (Ptr addr#) (W# x#)
       _        -> brokenInvariant
-    NatJ# m# -> case x of
+    NatJ# (BN# m#) -> case x of
       NatS# x# -> do
         poke (Ptr addr#) (W# x#)
         forM_ [1 .. sz - 1] $ \off ->
           pokeElemOff (Ptr addr#) off (0 :: Word)
-      NatJ# bn -> do
-        l <- exportBigNatToAddr bn addr# 0#
+      NatJ# (BN# bn) -> do
+        l <- IO (\token -> case bigNatToAddrLE# bn addr# token of (# newToken, l# #) -> (# newToken, W# l# #))
         forM_ [(fromIntegral :: Word -> Int) l .. (sz `shiftL` lgWordSize) - 1] $ \off ->
           pokeElemOff (Ptr addr#) off (0 :: Word8)
       where
-        sz = I# (sizeofBigNat# m#)
+        sz = I# (bigNatSize# m#)
   {-# INLINE poke #-}
 
 #ifdef MIN_VERSION_vector
 
+-- | No validation checks are performed;
+-- reading untrusted data may corrupt internal invariants.
 instance KnownNat m => P.Prim (Mod m) where
   sizeOf# x    = let !(I# sz#) = sizeOf x    in sz#
   {-# INLINE sizeOf# #-}
@@ -381,10 +462,10 @@
     NatS#{} -> Mod (NatS# w#)
       where
         !(W# w#) = P.indexByteArray# arr# i'
-    NatJ# m# -> Mod $ bigNatToNat $ importBigNatFromByteArray arr# (int2Word# i#) (int2Word# sz#) 0#
+    NatJ# (BN# m#) -> Mod $ bigNatToNat (runRW# (\token -> case bigNatFromByteArrayLE# (int2Word# sz#) arr# (int2Word# i#) token of (# _, bn# #) -> bn#))
       where
         !(I# lgWordSize#) = lgWordSize
-        sz# = sizeofBigNat# m# `iShiftL#` lgWordSize#
+        sz# = bigNatSize# m# `iShiftL#` lgWordSize#
         i# = i' *# sz#
   {-# INLINE indexByteArray# #-}
 
@@ -392,62 +473,63 @@
     NatS#{} -> Mod (NatS# w#)
       where
         !(W# w#) = P.indexOffAddr# arr# i'
-    NatJ# m# -> Mod $ bigNatToNat $ unsafeDupablePerformIO $ importBigNatFromAddr (arr# `plusAddr#` i#) (int2Word# sz#) 0#
+    NatJ# (BN# m#) -> Mod $ bigNatToNat (runRW# (\token -> case bigNatFromAddrLE# (int2Word# sz#) (arr# `plusAddr#` i#) token of (# _, bn# #) -> bn#))
       where
         !(I# lgWordSize#) = lgWordSize
-        sz# = sizeofBigNat# m# `iShiftL#` lgWordSize#
+        sz# = bigNatSize# m# `iShiftL#` lgWordSize#
         i# = i' *# sz#
   {-# INLINE indexOffAddr# #-}
 
   readByteArray# marr !i' token = case natVal' (proxy# :: Proxy# m) of
     NatS#{} -> case P.readByteArray# marr i' token of
       (# newToken, W# w# #) -> (# newToken, Mod (NatS# w#) #)
-    NatJ# m# -> case unsafeFreezeByteArray# marr token of
-      (# newToken, arr #) -> (# newToken, Mod (bigNatToNat (importBigNatFromByteArray arr (int2Word# i#) (int2Word# sz#) 0#)) #)
+    NatJ# (BN# m#) -> case unsafeFreezeByteArray# marr token of
+      (# newToken, arr #) -> case bigNatFromByteArrayLE# (int2Word# sz#) arr (int2Word# i#) newToken of
+        (# veryNewToken, bn# #) -> (# veryNewToken,Mod (bigNatToNat bn#) #)
       where
         !(I# lgWordSize#) = lgWordSize
-        sz# = sizeofBigNat# m# `iShiftL#` lgWordSize#
+        sz# = bigNatSize# m# `iShiftL#` lgWordSize#
         i# = i' *# sz#
   {-# INLINE readByteArray# #-}
 
   readOffAddr# marr !i' token = case natVal' (proxy# :: Proxy# m) of
     NatS#{} -> case P.readOffAddr# marr i' token of
       (# newToken, W# w# #) -> (# newToken, Mod (NatS# w#) #)
-    NatJ# m# -> case internal (unsafeIOToPrim (importBigNatFromAddr (marr `plusAddr#` i#) (int2Word# sz#) 0#) :: ST s BigNat) token of
+    NatJ# (BN# m#) -> case bigNatFromAddrLE# (int2Word# sz#) (marr `plusAddr#` i#) token of
       (# newToken, bn #) -> (# newToken, Mod (bigNatToNat bn) #)
       where
         !(I# lgWordSize#) = lgWordSize
-        sz# = sizeofBigNat# m# `iShiftL#` lgWordSize#
+        sz# = bigNatSize# m# `iShiftL#` lgWordSize#
         i# = i' *# sz#
   {-# INLINE readOffAddr# #-}
 
   writeByteArray# marr !i' !(Mod x) token = case natVal' (proxy# :: Proxy# m) of
     NatS#{} -> case x of
       NatS# x# -> P.writeByteArray# marr i' (W# x#) token
-      _        -> error "argument is larger than modulo"
-    NatJ# m# -> case x of
+      _        -> error "argument is larger than modulus"
+    NatJ# (BN# m#) -> case x of
       NatS# x# -> case P.writeByteArray# marr i# (W# x#) token of
         newToken -> P.setByteArray# marr (i# +# 1#) (sz# -# 1#) (0 :: Word) newToken
-      NatJ# bn -> case internal (unsafeIOToPrim (exportBigNatToMutableByteArray bn (unsafeCoerce# marr) (int2Word# (i# `iShiftL#` lgWordSize#)) 0#) :: ST s Word) token of
-        (# newToken, W# l# #) -> P.setByteArray# marr (i# `iShiftL#` lgWordSize# +# word2Int# l#) (sz# `iShiftL#` lgWordSize# -# word2Int# l#) (0 :: Word8) newToken
+      NatJ# (BN# bn) -> case bigNatToMutableByteArrayLE# bn (unsafeCoerce# marr) (int2Word# (i# `iShiftL#` lgWordSize#)) token of
+        (# newToken, l# #) -> P.setByteArray# marr (i# `iShiftL#` lgWordSize# +# word2Int# l#) (sz# `iShiftL#` lgWordSize# -# word2Int# l#) (0 :: Word8) newToken
       where
         !(I# lgWordSize#) = lgWordSize
-        !sz@(I# sz#) = I# (sizeofBigNat# m#)
+        !sz@(I# sz#) = I# (bigNatSize# m#)
         !(I# i#)     = I# i' * sz
   {-# INLINE writeByteArray# #-}
 
   writeOffAddr# marr !i' !(Mod x) token = case natVal' (proxy# :: Proxy# m) of
     NatS#{} -> case x of
       NatS# x# -> P.writeOffAddr# marr i' (W# x#) token
-      _        -> error "argument is larger than modulo"
-    NatJ# m# -> case x of
+      _        -> error "argument is larger than modulus"
+    NatJ# (BN# m#) -> case x of
       NatS# x# -> case P.writeOffAddr# marr i# (W# x#) token of
         newToken -> P.setOffAddr# marr (i# +# 1#) (sz# -# 1#) (0 :: Word) newToken
-      NatJ# bn -> case internal (unsafeIOToPrim (exportBigNatToAddr bn (marr `plusAddr#` (i# `iShiftL#` lgWordSize#)) 0#) :: ST s Word) token of
-        (# newToken, W# l# #) -> P.setOffAddr# marr (i# `iShiftL#` lgWordSize# +# word2Int# l#) (sz# `iShiftL#` lgWordSize# -# word2Int# l#) (0 :: Word8) newToken
+      NatJ# (BN# bn) -> case bigNatToAddrLE# bn (marr `plusAddr#` (i# `iShiftL#` lgWordSize#)) token of
+        (# newToken, l# #) -> P.setOffAddr# marr (i# `iShiftL#` lgWordSize# +# word2Int# l#) (sz# `iShiftL#` lgWordSize# -# word2Int# l#) (0 :: Word8) newToken
       where
         !(I# lgWordSize#) = lgWordSize
-        !sz@(I# sz#) = I# (sizeofBigNat# m#)
+        !sz@(I# sz#) = I# (bigNatSize# m#)
         !(I# i#)   = I# i' * sz
   {-# INLINE writeOffAddr# #-}
 
@@ -455,12 +537,12 @@
   setByteArray# marr off len mx@(Mod x) token = case natVal' (proxy# :: Proxy# m) of
     NatS#{} -> case x of
       NatS# x# -> P.setByteArray# marr off len (W# x#) token
-      _        -> error "argument is larger than modulo"
-    NatJ# m# -> case P.writeByteArray# marr off mx token of
+      _        -> error "argument is larger than modulus"
+    NatJ# (BN# m#) -> case P.writeByteArray# marr off mx token of
       newToken -> doSet (sz `iShiftL#` lgWordSize#) newToken
       where
         !(I# lgWordSize#) = lgWordSize
-        sz = sizeofBigNat# m#
+        sz = bigNatSize# m#
         off' = (off *# sz) `iShiftL#` lgWordSize#
         len' = (len *# sz) `iShiftL#` lgWordSize#
         doSet i tkn
@@ -473,12 +555,12 @@
   setOffAddr# marr off len mx@(Mod x) token = case natVal' (proxy# :: Proxy# m) of
     NatS#{} -> case x of
       NatS# x# -> P.setOffAddr# marr off len (W# x#) token
-      _        -> error "argument is larger than modulo"
-    NatJ# m# -> case P.writeOffAddr# marr off mx token of
+      _        -> error "argument is larger than modulus"
+    NatJ# (BN# m#) -> case P.writeOffAddr# marr off mx token of
       newToken -> doSet (sz `iShiftL#` lgWordSize#) newToken
       where
         !(I# lgWordSize#) = lgWordSize
-        sz = sizeofBigNat# m#
+        sz = bigNatSize# m#
         off' = (off *# sz) `iShiftL#` lgWordSize#
         len' = (len *# sz) `iShiftL#` lgWordSize#
         doSet i tkn -- = tkn
@@ -489,17 +571,21 @@
   {-# INLINE setOffAddr# #-}
 
 -- | Unboxed vectors of 'Mod' cause more nursery allocations
--- than boxed ones, but reduce pressure on garbage collector,
+-- than boxed ones, but reduce pressure on the garbage collector,
 -- especially for large vectors.
 newtype instance U.MVector s (Mod m) = ModMVec (P.MVector s (Mod m))
 
 -- | Unboxed vectors of 'Mod' cause more nursery allocations
--- than boxed ones, but reduce pressure on garbage collector,
+-- than boxed ones, but reduce pressure on the garbage collector,
 -- especially for large vectors.
 newtype instance U.Vector    (Mod m) = ModVec  (P.Vector (Mod m))
 
+-- | No validation checks are performed;
+-- reading untrusted data may corrupt internal invariants.
 instance KnownNat m => U.Unbox (Mod m)
 
+-- | No validation checks are performed;
+-- reading untrusted data may corrupt internal invariants.
 instance KnownNat m => M.MVector U.MVector (Mod m) where
   {-# INLINE basicLength #-}
   {-# INLINE basicUnsafeSlice #-}
@@ -527,6 +613,8 @@
   basicUnsafeMove (ModMVec v1) (ModMVec v2) = M.basicUnsafeMove v1 v2
   basicUnsafeGrow (ModMVec v) n = ModMVec `liftM` M.basicUnsafeGrow v n
 
+-- | No validation checks are performed;
+-- reading untrusted data may corrupt internal invariants.
 instance KnownNat m => G.Vector U.Vector (Mod m) where
   {-# INLINE basicUnsafeFreeze #-}
   {-# INLINE basicUnsafeThaw #-}
diff --git a/Data/Mod/Compat.hs b/Data/Mod/Compat.hs
new file mode 100644
--- /dev/null
+++ b/Data/Mod/Compat.hs
@@ -0,0 +1,46 @@
+-- | See https://gitlab.haskell.org/ghc/ghc/-/issues/22933
+--   and https://gitlab.haskell.org/ghc/ghc/-/issues/22966
+--
+
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE CApiFFI #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE MagicHash #-}
+{-# LANGUAGE UnboxedTuples #-}
+{-# LANGUAGE UnliftedFFITypes #-}
+
+module Data.Mod.Compat
+  ( timesWord2#
+  , remWord2#
+  ) where
+
+#ifdef aarch64_HOST_ARCH
+import GHC.Exts (Word(..), Word#, timesWord#)
+
+timesWord2# :: Word# -> Word# -> (# Word#, Word# #)
+timesWord2# x y = (# z, timesWord# x y #)
+  where
+    !(W# z) = c_umulh (W# x) (W# y)
+{-# INLINE timesWord2# #-}
+
+foreign import capi unsafe "aarch64.h umulh" c_umulh :: Word -> Word -> Word
+
+remWord2# :: Word# -> Word# -> Word# -> Word#
+remWord2# lo hi m = r
+  where
+    !(W# r) = c_umodh (W# lo) (W# hi) (W# m)
+{-# INLINE remWord2# #-}
+
+foreign import capi unsafe "aarch64.h umodh" c_umodh :: Word -> Word -> Word -> Word
+
+#else
+
+import GHC.Exts (Word#, timesWord2#, quotRemWord2#)
+
+remWord2# :: Word# -> Word# -> Word# -> Word#
+remWord2# lo hi m = r
+  where
+    !(# _, r #) = quotRemWord2# hi lo m
+{-# INLINE remWord2# #-}
+
+#endif
diff --git a/Data/Mod/Word.hs b/Data/Mod/Word.hs
--- a/Data/Mod/Word.hs
+++ b/Data/Mod/Word.hs
@@ -1,25 +1,26 @@
 -- |
 -- Module:      Data.Mod.Word
--- Copyright:   (c) 2017-2020 Andrew Lelechenko
+-- Copyright:   (c) 2017-2022 Andrew Lelechenko
 -- Licence:     MIT
 -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
 --
 -- <https://en.wikipedia.org/wiki/Modular_arithmetic Modular arithmetic>,
 -- promoting moduli to the type level, with an emphasis on performance.
--- Originally part of <https://hackage.haskell.org/package/arithmoi arithmoi> package.
+-- Originally part of the <https://hackage.haskell.org/package/arithmoi arithmoi> package.
 --
 -- This module supports only moduli, which fit into 'Word'.
--- Use (slower) "Data.Mod" to handle arbitrary-sized moduli.
+-- Use the (slower) "Data.Mod" module for handling arbitrary-sized moduli.
 
 {-# LANGUAGE BangPatterns               #-}
 {-# LANGUAGE CPP                        #-}
+{-# LANGUAGE DataKinds                  #-}
 {-# LANGUAGE DeriveGeneric              #-}
+{-# LANGUAGE DerivingStrategies         #-}
 {-# LANGUAGE GeneralizedNewtypeDeriving #-}
 {-# LANGUAGE MagicHash                  #-}
 {-# LANGUAGE MultiParamTypeClasses      #-}
 {-# LANGUAGE TypeApplications           #-}
 {-# LANGUAGE TypeFamilies               #-}
-{-# LANGUAGE TypeInType                 #-}
 {-# LANGUAGE UnboxedTuples              #-}
 
 module Data.Mod.Word
@@ -33,9 +34,10 @@
 import Control.Exception
 import Control.DeepSeq
 import Data.Bits
+import Data.Mod.Compat (timesWord2#, remWord2#)
+import Data.Ratio
 #ifdef MIN_VERSION_semirings
 import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)
-import Data.Ratio
 import Data.Semiring (Semiring(..), Ring(..))
 #endif
 #ifdef MIN_VERSION_vector
@@ -46,11 +48,13 @@
 import qualified Data.Vector.Unboxed         as U
 #endif
 import Foreign.Storable (Storable)
-import GHC.Exts
+import GHC.Exts hiding (timesWord2#, quotRemWord2#)
 import GHC.Generics
-import GHC.Integer.GMP.Internals
 import GHC.Natural (Natural(..))
+import GHC.Num.BigNat
+import GHC.Num.Integer
 import GHC.TypeNats (Nat, KnownNat, natVal)
+import Text.Read (Read(readPrec))
 
 -- | This data type represents
 -- <https://en.wikipedia.org/wiki/Modular_arithmetic#Integers_modulo_n integers modulo m>,
@@ -61,32 +65,41 @@
 --
 -- >>> :set -XDataKinds
 -- >>> 3 + 8 :: Mod 10 -- 3 + 8 = 11 ≡ 1 (mod 10)
--- (1 `modulo` 10)
+-- 1
 --
--- __Warning:__ division by residue, which is not
--- <https://en.wikipedia.org/wiki/Coprime_integers coprime>
--- with the modulo, throws 'DivideByZero'.
--- Consider using 'invertMod' for non-prime moduli.
+-- __Note:__ 'Mod' 0 has no inhabitants, eventhough \( \mathbb{Z}/0\mathbb{Z} \) is technically isomorphic to \( \mathbb{Z} \).
 newtype Mod (m :: Nat) = Mod
   { unMod :: Word
   -- ^ The canonical representative of the residue class,
-  -- always between 0 and \( m - 1 \) inclusively.
+  -- always between 0 and \( m - 1 \) (inclusively).
   --
   -- >>> :set -XDataKinds
   -- >>> -1 :: Mod 10
-  -- (9 `modulo` 10)
+  -- 9
   }
+  deriving (Eq, Ord, Generic)
+  deriving Storable
+  -- ^ No validation checks are performed;
+  -- reading untrusted data may corrupt internal invariants.
 #ifdef MIN_VERSION_vector
-  deriving (Eq, Ord, Generic, Storable, Prim)
-#else
-  deriving (Eq, Ord, Generic, Storable)
+  deriving Prim
+  -- ^ No validation checks are performed;
+  -- reading untrusted data may corrupt internal invariants.
 #endif
 
 instance NFData (Mod m)
 
-instance KnownNat m => Show (Mod m) where
-  show m = "(" ++ show (unMod m) ++ " `modulo` " ++ show (natVal m) ++ ")"
+instance Show (Mod m) where
+  show (Mod x) = show x
 
+-- | Wrapping behaviour, similar to
+-- the existing @instance@ 'Read' 'Int'.
+instance KnownNat m => Read (Mod m) where
+  readPrec = fromInteger <$> readPrec
+
+instance KnownNat m => Real (Mod m) where
+  toRational (Mod x) = toRational x
+
 instance KnownNat m => Enum (Mod m) where
   succ x = if x == maxBound then throw Overflow  else coerce (succ @Word) x
   pred x = if x == minBound then throw Underflow else coerce (pred @Word) x
@@ -101,66 +114,73 @@
   enumFromThenTo = coerce (enumFromThenTo @Word)
 
 instance KnownNat m => Bounded (Mod m) where
-  minBound = Mod 0
-  maxBound = let mx = Mod (fromIntegral (natVal mx) - 1) in mx
-
-#if !MIN_VERSION_base(4,12,0)
-addWordC# :: Word# -> Word# -> (# Word#, Int# #)
-addWordC# x# y# = (# z#, word2Int# c# #)
-  where
-    !(# c#, z# #) = x# `plusWord2#` y#
-#endif
+  minBound = mx
+    where
+      mx = if natVal mx > 0 then Mod 0 else throw DivideByZero
+  maxBound = mx
+    where
+      mx = if m > 0 then Mod (fromIntegral (m - 1)) else throw DivideByZero
+      m = natVal mx
 
 addMod :: Natural -> Word -> Word -> Word
 addMod (NatS# m#) (W# x#) (W# y#) =
   if isTrue# c# || isTrue# (z# `geWord#` m#) then W# (z# `minusWord#` m#) else W# z#
   where
     !(# z#, c# #) = x# `addWordC#` y#
-addMod NatJ#{} _ _ = tooLargeModulo
+addMod NatJ#{} _ _ = tooLargeModulus
 
 subMod :: Natural -> Word -> Word -> Word
 subMod (NatS# m#) (W# x#) (W# y#) =
   if isTrue# (x# `geWord#` y#) then W# z# else W# (z# `plusWord#` m#)
   where
     z# = x# `minusWord#` y#
-subMod NatJ#{} _ _ = tooLargeModulo
+subMod NatJ#{} _ _ = tooLargeModulus
 
 negateMod :: Natural -> Word -> Word
 negateMod _ (W# 0##) = W# 0##
 negateMod (NatS# m#) (W# x#) = W# (m# `minusWord#` x#)
-negateMod NatJ#{} _ = tooLargeModulo
+negateMod NatJ#{} _ = tooLargeModulus
 
+halfWord :: Word
+halfWord = 1 `shiftL` (finiteBitSize (0 :: Word) `shiftR` 1)
+
 mulMod :: Natural -> Word -> Word -> Word
-mulMod (NatS# m#) (W# x#) (W# y#) = W# r#
+mulMod (NatS# m#) (W# x#) (W# y#)
+  | W# m# <= halfWord = W# (timesWord# x# y# `remWord#` m#)
+  | otherwise = W# r#
   where
-    !(# z1#, z2# #) = timesWord2# x# y#
-    !(# _, r# #) = quotRemWord2# z1# z2# m#
-mulMod NatJ#{} _ _ = tooLargeModulo
+    !(# hi#, lo# #) = timesWord2# x# y#
+    !r# = remWord2# lo# hi# m#
+mulMod NatJ#{} _ _ = tooLargeModulus
 
 fromIntegerMod :: Natural -> Integer -> Word
 fromIntegerMod (NatS# 0##) !_ = throw DivideByZero
-fromIntegerMod (NatS# m#) (S# x#) =
+fromIntegerMod (NatS# m#) (IS x#) =
   if isTrue# (x# >=# 0#)
     then W# (int2Word# x# `remWord#` m#)
     else negateMod (NatS# m#) (W# (int2Word# (negateInt# x#) `remWord#` m#))
-fromIntegerMod (NatS# m#) (Jp# x#) =
-  W# (x# `remBigNatWord` m#)
-fromIntegerMod (NatS# m#) (Jn# x#) =
-  negateMod (NatS# m#) (W# (x# `remBigNatWord` m#))
-fromIntegerMod NatJ#{} _ = tooLargeModulo
+fromIntegerMod (NatS# m#) (IP x#) =
+  W# (x# `bigNatRemWord#` m#)
+fromIntegerMod (NatS# m#) (IN x#) =
+  negateMod (NatS# m#) (W# (x# `bigNatRemWord#` m#))
+fromIntegerMod NatJ#{} _ = tooLargeModulus
 
 #ifdef MIN_VERSION_semirings
 
 fromNaturalMod :: Natural -> Natural -> Word
 fromNaturalMod (NatS# 0##) !_ = throw DivideByZero
 fromNaturalMod (NatS# m#) (NatS# x#) = W# (x# `remWord#` m#)
-fromNaturalMod (NatS# m#) (NatJ# x#) = W# (x# `remBigNatWord` m#)
-fromNaturalMod NatJ#{} _ = tooLargeModulo
+fromNaturalMod (NatS# m#) (NatJ# (BN# x#)) = W# (x# `bigNatRemWord#` m#)
+fromNaturalMod NatJ#{} _ = tooLargeModulus
 
+getModulus :: Natural -> Word
+getModulus (NatS# m#) = W# m#
+getModulus NatJ#{} = tooLargeModulus
+
 #endif
 
-tooLargeModulo :: a
-tooLargeModulo = error "modulo does not fit into a machine word"
+tooLargeModulus :: a
+tooLargeModulus = error "modulus does not fit into a machine word"
 
 instance KnownNat m => Num (Mod m) where
   mx@(Mod !x) + (Mod !y) = Mod $ addMod (natVal mx) x y
@@ -189,11 +209,17 @@
   {-# INLINE plus #-}
   times = (*)
   {-# INLINE times #-}
-  zero  = Mod 0
+  zero  = mx
+    where
+      mx = if natVal mx > 0 then Mod 0 else throw DivideByZero
   {-# INLINE zero #-}
   one   = mx
     where
-      mx = if natVal mx > 1 then Mod 1 else Mod 0
+      mx = case m `compare` 1 of
+        LT -> throw DivideByZero
+        EQ -> Mod 0
+        GT -> Mod 1
+      m = natVal mx
   {-# INLINE one #-}
   fromNatural x = mx
     where
@@ -204,7 +230,85 @@
   negate = P.negate
   {-# INLINE negate #-}
 
--- | See the warning about division above.
+-- | 'Mod' @m@ is not even an
+-- <https://en.wikipedia.org/wiki/Integral_domain integral domain> for
+-- <https://en.wikipedia.org/wiki/Composite_number composite> @m@,
+-- much less a <https://en.wikipedia.org/wiki/GCD_domain GCD domain>.
+-- However, 'Data.Euclidean.gcd' and 'Data.Euclidean.lcm' are still meaningful
+-- even for composite @m@, corresponding to a sum and an intersection of
+-- <https://en.wikipedia.org/wiki/Ideal_(ring_theory) ideals>.
+--
+-- The instance is lawful only for
+-- <https://en.wikipedia.org/wiki/Prime_number prime> @m@, otherwise
+-- @'Data.Euclidean.divide' x y@ tries to return any @Just z@ such that @x == y * z@.
+--
+instance KnownNat m => GcdDomain (Mod m) where
+  divide (Mod 0) !_ = Just (Mod 0)
+  divide _ (Mod 0) = Nothing
+  divide mx@(Mod x) (Mod y) = case mry of
+    Just ry -> if xr == 0 then Just (Mod xq * Mod ry) else Nothing
+    Nothing -> Nothing
+    where
+      m = getModulus (natVal mx)
+      gmy = P.gcd m y
+      (xq, xr) = P.quotRem x gmy
+      mry = invertModWord (y `P.quot` gmy)  (m `P.quot` gmy)
+
+  gcd (Mod !x) (Mod !y) = g
+    where
+      m = getModulus (natVal g)
+      g = Mod $ if m > 1 then P.gcd (P.gcd m x) y else 0
+  lcm (Mod !x) (Mod !y) = l
+    where
+      m = getModulus (natVal l)
+      l = Mod $ if m > 1 then P.lcm (P.gcd m x) (P.gcd m y) else 0
+  coprime x y = Data.Euclidean.gcd x y == one
+
+-- | 'Mod' @m@ is not even an
+-- <https://en.wikipedia.org/wiki/Integral_domain integral domain> for
+-- <https://en.wikipedia.org/wiki/Composite_number composite> @m@,
+-- much less a <https://en.wikipedia.org/wiki/Euclidean_domain Euclidean domain>.
+--
+-- The instance is lawful only for
+-- <https://en.wikipedia.org/wiki/Prime_number prime> @m@, otherwise
+-- we try to do our best:
+-- @'Data.Euclidean.quot' x y@ returns any @z@ such that @x == y * z@,
+-- 'Data.Euclidean.rem' is not always 0, and both can throw 'DivideByZero'.
+--
+instance KnownNat m => Euclidean (Mod m) where
+  degree = fromIntegral . unMod
+
+  quotRem (Mod 0) !_ = (Mod 0, Mod 0)
+  quotRem _ (Mod 0) = throw DivideByZero
+  quotRem mx@(Mod x) (Mod y) = case mry of
+    Just ry -> (Mod xq * Mod ry, Mod xr)
+    Nothing -> throw DivideByZero
+    where
+      m = getModulus (natVal mx)
+      gmy = P.gcd m y
+      (xq, xr) = P.quotRem x gmy
+      mry = invertModWord (y `P.quot` gmy)  (m `P.quot` gmy)
+
+-- | 'Mod' @m@ is not even an
+-- <https://en.wikipedia.org/wiki/Integral_domain integral domain> for
+-- <https://en.wikipedia.org/wiki/Composite_number composite> @m@,
+-- much less a <https://en.wikipedia.org/wiki/Field_(mathematics) field>.
+--
+-- The instance is lawful only for
+-- <https://en.wikipedia.org/wiki/Prime_number prime> @m@, otherwise
+-- division by a residue, which is not
+-- <https://en.wikipedia.org/wiki/Coprime_integers coprime>
+-- with the modulus, throws 'DivideByZero'.
+-- Consider using 'invertMod' for non-prime moduli.
+--
+instance KnownNat m => Field (Mod m)
+
+#endif
+
+-- | Division by a residue, which is not
+-- <https://en.wikipedia.org/wiki/Coprime_integers coprime>
+-- with the modulus, throws 'DivideByZero'.
+-- Consider using 'invertMod' for non-prime moduli.
 instance KnownNat m => Fractional (Mod m) where
   fromRational r = case denominator r of
     1   -> num
@@ -217,44 +321,25 @@
     Just y  -> y
   {-# INLINE recip #-}
 
--- | See the warning about division above.
-instance KnownNat m => GcdDomain (Mod m) where
-  divide x y = Just (x / y)
-  gcd        = const $ const 1
-  lcm        = const $ const 1
-  coprime    = const $ const True
-
--- | See the warning about division above.
-instance KnownNat m => Euclidean (Mod m) where
-  degree      = const 0
-  quotRem x y = (x / y, 0)
-  quot        = (/)
-  rem         = const $ const 0
-
--- | See the warning about division above.
-instance KnownNat m => Field (Mod m)
-
-#endif
-
 -- | If an argument is
 -- <https://en.wikipedia.org/wiki/Coprime_integers coprime>
--- with the modulo, return its modular inverse.
+-- with the modulus, return its modular inverse.
 -- Otherwise return 'Nothing'.
 --
 -- >>> :set -XDataKinds
 -- >>> invertMod 3 :: Mod 10 -- 3 * 7 = 21 ≡ 1 (mod 10)
--- Just (7 `modulo` 10)
+-- Just 7
 -- >>> invertMod 4 :: Mod 10 -- 4 and 10 are not coprime
 -- Nothing
 invertMod :: KnownNat m => Mod m -> Maybe (Mod m)
-invertMod mx@(Mod x) = case natVal mx of
-  NatJ#{}   -> tooLargeModulo
+invertMod mx@(Mod !x) = case natVal mx of
+  NatJ#{}   -> tooLargeModulus
   NatS# 0## -> Nothing
   NatS# m#  -> Mod <$> invertModWord x (W# m#)
 
 invertModWord :: Word -> Word -> Maybe Word
 invertModWord x m@(W# m#)
-  -- If both x and k are even, no inverse exists
+  -- If both x and m are even, no inverse exists
   | even x, isTrue# (k# `gtWord#` 0##) = Nothing
   | otherwise = case invertModWordOdd x m' of
     Nothing -> Nothing
@@ -335,34 +420,26 @@
 
 -- | Drop-in replacement for 'Prelude.^' with a bit better performance.
 -- Negative powers are allowed, but may throw 'DivideByZero', if an argument
--- is not <https://en.wikipedia.org/wiki/Coprime_integers coprime> with the modulo.
---
--- Building with @-O@ triggers a rewrite rule 'Prelude.^' = '^%'.
+-- is not <https://en.wikipedia.org/wiki/Coprime_integers coprime> with the modulus.
 --
 -- >>> :set -XDataKinds
 -- >>> 3 ^% 4 :: Mod 10    -- 3 ^ 4 = 81 ≡ 1 (mod 10)
--- (1 `modulo` 10)
+-- 1
 -- >>> 3 ^% (-1) :: Mod 10 -- 3 * 7 = 21 ≡ 1 (mod 10)
--- (7 `modulo` 10)
+-- 7
 -- >>> 4 ^% (-1) :: Mod 10 -- 4 and 10 are not coprime
 -- (*** Exception: divide by zero
 (^%) :: (KnownNat m, Integral a) => Mod m -> a -> Mod m
-mx@(Mod (W# x#)) ^% a = case natVal mx of
-  NatJ#{} -> tooLargeModulo
-  NatS# m#
+mx@(Mod !x) ^% a = case natVal mx of
+  NatJ#{} -> tooLargeModulus
+  m@(NatS# _)
     | a < 0 -> case invertMod mx of
-      Nothing            -> throw DivideByZero
-      Just (Mod (W# y#)) -> Mod $ W# (f y# (- a) 1##)
-    | otherwise          -> Mod $ W# (f x# a 1##)
+      Nothing      -> throw DivideByZero
+      Just (Mod y) -> Mod $ f y (-a) 1
+    | otherwise    -> Mod $ f x a 1
     where
-      f :: Integral a => Word# -> a -> Word# -> Word#
-      f _  0 acc# = acc#
-      f b# e acc# = f bb# (e `P.quot` 2) (if odd e then ba# else acc#)
-        where
-          !(# bb1#, bb2# #) = timesWord2# b# b#
-          !(#    _, bb#  #) = quotRemWord2# bb1# bb2# m#
-          !(# ba1#, ba2# #) = timesWord2# b# acc#
-          !(#    _, ba#  #) = quotRemWord2# ba1# ba2# m#
+      f !_ 0 acc = acc
+      f b  e acc = f (mulMod m b b) (e `P.quot` 2) (if odd e then mulMod m b acc else acc)
 {-# INLINABLE [1] (^%) #-}
 
 {-# SPECIALISE [1] (^%) ::
@@ -372,8 +449,6 @@
   KnownNat m => Mod m -> Word    -> Mod m #-}
 
 {-# RULES
-"powMod"               forall (x :: KnownNat m => Mod m) p. x ^ p = x ^% p
-
 "powMod/2/Integer"     forall x. x ^% (2 :: Integer) = let u = x in u*u
 "powMod/3/Integer"     forall x. x ^% (3 :: Integer) = let u = x in u*u*u
 "powMod/2/Int"         forall x. x ^% (2 :: Int)     = let u = x in u*u
@@ -388,8 +463,12 @@
 newtype instance U.MVector s (Mod m) = MV_Mod (P.MVector s Word)
 newtype instance U.Vector    (Mod m) = V_Mod  (P.Vector    Word)
 
+-- | No validation checks are performed;
+-- reading untrusted data may corrupt internal invariants.
 instance U.Unbox (Mod m)
 
+-- | No validation checks are performed;
+-- reading untrusted data may corrupt internal invariants.
 instance M.MVector U.MVector (Mod m) where
   {-# INLINE basicLength #-}
   {-# INLINE basicUnsafeSlice #-}
@@ -417,6 +496,8 @@
   basicUnsafeMove (MV_Mod v1) (MV_Mod v2) = M.basicUnsafeMove v1 v2
   basicUnsafeGrow (MV_Mod v) n = MV_Mod <$> M.basicUnsafeGrow v n
 
+-- | No validation checks are performed;
+-- reading untrusted data may corrupt internal invariants.
 instance G.Vector U.Vector (Mod m) where
   {-# INLINE basicUnsafeFreeze #-}
   {-# INLINE basicUnsafeThaw #-}
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2017-2020 Andrew Lelechenko
+Copyright (c) 2017-2022 Andrew Lelechenko
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
  associated documentation files (the "Software"), to deal in the Software without restriction,
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -1,21 +1,21 @@
-# mod [![Build Status](https://github.com/Bodigrim/mod/workflows/Haskell-CI/badge.svg)](https://github.com/Bodigrim/mod/actions?query=workflow%3AHaskell-CI) [![Hackage](http://img.shields.io/hackage/v/mod.svg)](https://hackage.haskell.org/package/mod) [![Hackage CI](https://matrix.hackage.haskell.org/api/v2/packages/mod/badge)](https://matrix.hackage.haskell.org/package/mod) [![Stackage LTS](http://stackage.org/package/mod/badge/lts)](http://stackage.org/lts/package/mod) [![Stackage Nightly](http://stackage.org/package/mod/badge/nightly)](http://stackage.org/nightly/package/mod)
+# mod [![Hackage](https://img.shields.io/hackage/v/mod.svg)](https://hackage.haskell.org/package/mod) [![Stackage LTS](https://www.stackage.org/package/mod/badge/lts)](https://www.stackage.org/lts/package/mod) [![Stackage Nightly](https://www.stackage.org/package/mod/badge/nightly)](https://www.stackage.org/nightly/package/mod)
 
 [Modular arithmetic](https://en.wikipedia.org/wiki/Modular_arithmetic),
 promoting moduli to the type level, with an emphasis on performance.
-Originally a part of [arithmoi](https://hackage.haskell.org/package/arithmoi) package.
+Originally a part of the [arithmoi](https://hackage.haskell.org/package/arithmoi) package.
 
 ```haskell
 > :set -XDataKinds
 > 4 + 5 :: Mod 7
-(2 `modulo` 7)
+2
 > 4 - 5 :: Mod 7
-(6 `modulo` 7)
+6
 > 4 * 5 :: Mod 7
-(6 `modulo` 7)
+6
 > 4 / 5 :: Mod 7
-(5 `modulo` 7)
+5
 > 4 ^ 5 :: Mod 7
-(2 `modulo` 7)
+2
 ```
 
 ## Competitors
@@ -37,12 +37,12 @@
 * __Addition.__
   All competing implementations of
   the modular addition involve divisions, while `mod` completely avoids
-  this costly operation. It makes difference even for small numbers;
+  this costly operation. This makes a difference even for small numbers;
   e. g., `sum [1..10^7]` becomes 5x faster. For larger integers the speed up
   is even more significant, because the computational complexity of division is not linear.
 
 * __Small `(*)`.__
-  When a modulo fits a machine word (which is quite a common case on 64-bit architectures),
+  When a modulus fits in a machine word (which is quite a common case on 64-bit architectures),
   `mod` implements the modular multiplication as a couple of CPU instructions
   and neither allocates intermediate arbitrary-precision values,
   nor calls `libgmp` at all. For computations like `product [1..10^7]`
@@ -64,7 +64,7 @@
   because it allows to specify the underlying representation of a modular residue,
   e. g., `Mod Integer 100`, `Mod Int 100`, `Mod Word8 100`. We argue that this is
   a dangerous freedom, vulnerable to overflows.
-  For instance, `20 ^ 2 :: Mod Word8 100` returns `44` instead of expected `0`.
+  For instance, `20 ^ 2 :: Mod Word8 100` returns `44` instead of the expected `0`.
   Even less expected is that `50 :: Mod Word8 300` appears to be `6`
   (remember that type-level numbers are always `Natural`).
 
@@ -78,7 +78,7 @@
 like indices of arrays in `vector-sized`.
 It features a `Num` instance only for the sake of overloading numeric literals.
 There is no lawful way to define `Num` except modular arithmetic,
-but from `finite-typelits` viewpoint this is a by-product.
+but from `finite-typelits`' viewpoint this is a by-product.
 
 ## Citius, altius, fortius!
 
@@ -95,16 +95,16 @@
 
 | Discipline  | `Data.Mod.Word`  | `Data.Mod`  | `modular` | `modular-arithmetic` | `finite-typelits` | `finite-field`
 | :---------- | :--------------: | :---------: | :-------: | :------------------: | :---------------: | :------------:
-| Sum         |   0.25x           |    1x       |  11.4x    |      5.7x            |  8.9x             | 8.6x
-| Product     |   0.95x           |    1x       |  9.6x     |      4.8x            |  7.0x             | 7.0x
-| Inversion   |   0.95x           |    1x       |  N/A      |      2.6x            |  N/A              | 3.0x
-| Power       |   0.90x           |    1x       |  6.9x     |      3.8x            |  5.0x             | 4.9x
+| Sum         |   0.44x          |    1x       |  16.6x    |      8.9x            |  14.7x            | 14.2x
+| Product     |   0.95x          |    1x       |  7.8x     |      4.5x            |  7.0x             | 7.0x
+| Inversion   |   0.54x          |    1x       |  N/A      |      3.2x            |  N/A              | 1.8x
+| Power       |   0.29x          |    1x       |  2.0x     |      1.2x            |  1.4x             | 1.5x
 
 ## What's next?
 
 This package was cut out of [`arithmoi`](https://hackage.haskell.org/package/arithmoi)
-to provide a modular arithmetic
-with a light dependency footprint. This goal certainly limits the scope of API
+to provide modular arithmetic
+with a light dependency footprint. This goal certainly limits the scope of the API
 to the bare minimum. If you need more advanced tools
 (the Chinese remainder theorem, cyclic groups, modular equations, etc.)
-please refer to [Math.NumberTheory.Moduli](https://hackage.haskell.org/package/arithmoi/docs/Math-NumberTheory-Moduli.html).
+please refer to the [Math.NumberTheory.Moduli](https://hackage.haskell.org/package/arithmoi/docs/Math-NumberTheory-Moduli.html) module.
diff --git a/Setup.hs b/Setup.hs
deleted file mode 100644
--- a/Setup.hs
+++ /dev/null
@@ -1,5 +0,0 @@
-module Main where
-
-import Distribution.Simple
-
-main = defaultMain
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -60,9 +60,10 @@
     measure name p = bench name $ whnf (sumN p) lim
     {-# INLINE measure #-}
 
-    sumN :: (Eq (t P), Num (t P)) => Proxy t -> Int -> t P
+    sumN :: forall t. (Eq (t P), Num (t P)) => Proxy t -> Int -> t P
     sumN = const $ \n -> go 0 (fromIntegral n)
       where
+        go :: t P -> t P -> t P
         go !acc 0 = acc
         go acc n = go (acc + n) (n - 1)
     {-# INLINE sumN #-}
@@ -71,6 +72,7 @@
     sumNModular :: Int -> Numeric.Modular.Mod P
     sumNModular = \n -> go 0 (fromIntegral n)
       where
+        go :: Numeric.Modular.Mod P -> Numeric.Modular.Mod P -> Numeric.Modular.Mod P
         go acc@(forceModular -> !_) 0 = acc
         go acc n = go (acc + n) (n - 1)
     {-# INLINE sumNModular #-}
@@ -101,9 +103,10 @@
     measure name p = bench name $ whnf (productN p) lim
     {-# INLINE measure #-}
 
-    productN :: (Eq (t P), Num (t P)) => Proxy t -> Int -> t P
+    productN :: forall t. (Eq (t P), Num (t P)) => Proxy t -> Int -> t P
     productN = const $ \n -> go 1 (fromIntegral n)
       where
+        go :: t P -> t P -> t P
         go !acc 0 = acc
         go acc n = go (acc * n) (n - 1)
     {-# INLINE productN #-}
@@ -112,6 +115,7 @@
     productNModular :: Int -> Numeric.Modular.Mod P
     productNModular = \n -> go 1 (fromIntegral n)
       where
+        go :: Numeric.Modular.Mod P -> Numeric.Modular.Mod P -> Numeric.Modular.Mod P
         go acc@(forceModular -> !_) 0 = acc
         go acc n = go (acc * n) (n - 1)
     {-# INLINE productNModular #-}
@@ -136,17 +140,18 @@
     measure name p = bench name $ whnf (invertN p) lim
     {-# INLINE measure #-}
 
-    invertN :: (Eq (t P), Fractional (t P)) => Proxy t -> Int -> t P
+    invertN :: forall t. (Eq (t P), Fractional (t P)) => Proxy t -> Int -> t P
     invertN = const $ \n -> go 0 (fromIntegral n)
       where
+        go :: t P -> t P -> t P
         go !acc 0 = acc
         go acc n = go (acc + recip n) (n - 1)
     {-# INLINE invertN #-}
 
 benchPower :: Benchmark
 benchPower = bgroup "Power"
-  [ measure "Data.Mod" (Proxy @Data.Mod.Mod)
-  , cmp $ measure "Data.Mod.Word" (Proxy @Data.Mod.Word.Mod)
+  [ bench "Data.Mod" $ nf powerNMod lim
+  , cmp $ bench "Data.Mod.Word" $ nf powerNModWord lim
 #ifdef MIN_VERSION_finite_field
   , cmp $ measure "finite-field" (Proxy @Data.FiniteField.PrimeField.PrimeField)
 #endif
@@ -164,21 +169,41 @@
     cmp = bcompare "$NF == \"Data.Mod\" && $(NF-1) == \"Power\""
     lim = 1000000
 
+    powerNMod :: Int -> Data.Mod.Mod P
+    powerNMod = go 0
+      where
+        go :: Data.Mod.Mod P -> Int -> Data.Mod.Mod P
+        go !acc 0 = acc
+        go acc n = go (acc + 2 Data.Mod.^% n) (n - 1)
+    {-# INLINE powerNMod #-}
+
+    powerNModWord :: Int -> Data.Mod.Word.Mod P
+    powerNModWord = go 0
+      where
+        go :: Data.Mod.Word.Mod P -> Int -> Data.Mod.Word.Mod P
+        go !acc 0 = acc
+        go acc n = go (acc + 2 Data.Mod.Word.^% n) (n - 1)
+    {-# INLINE powerNModWord #-}
+
+#if defined(MIN_VERSION_finite_field) || defined(MIN_VERSION_modular_arithmetic)
     measure :: (Eq (t P), Num (t P)) => String -> Proxy t -> Benchmark
     measure name p = bench name $ whnf (powerN p) lim
     {-# INLINE measure #-}
 
-    powerN :: (Eq (t P), Num (t P)) => Proxy t -> Int -> t P
+    powerN :: forall t. (Eq (t P), Num (t P)) => Proxy t -> Int -> t P
     powerN = const $ go 0
       where
+        go :: t P -> Int -> t P
         go !acc 0 = acc
         go acc n = go (acc + 2 ^ n) (n - 1)
     {-# INLINE powerN #-}
+#endif
 
 #ifdef MIN_VERSION_modular
     powerNModular :: Int -> Numeric.Modular.Mod P
     powerNModular = go 0
       where
+        go :: Numeric.Modular.Mod P -> Int -> Numeric.Modular.Mod P
         go acc@(forceModular -> !_) 0 = acc
         go acc n = go (acc + 2 ^ n) (n - 1)
     {-# INLINE powerNModular #-}
diff --git a/cbits/aarch64.c b/cbits/aarch64.c
new file mode 100644
--- /dev/null
+++ b/cbits/aarch64.c
@@ -0,0 +1,9 @@
+#include <stdint.h>
+
+uint64_t umulh(uint64_t x, uint64_t y) {
+  return ((unsigned __int128)x * y) >> 64;
+}
+
+uint64_t umodh(uint64_t lo, uint64_t hi, uint64_t m) {
+  return (((unsigned __int128)hi << 64) + lo) % m;
+}
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,3 +1,12 @@
+# 0.2.0.0
+
+* Breaking change: redesign `GcdDomain` and `Euclidean` instances.
+* Add `instance Read` and `instance Real`.
+* Migrate from `integer-gmp` to `ghc-bignum`.
+* Remove `(^) -> (^%)` rewrite rule, it does not fire.
+* Plug loopholes to inhabit `Mod 0`.
+* Work around performance degradation on ARM.
+
 # 0.1.2.2
 
 * Work around an issue with [`fromIntegral`](https://gitlab.haskell.org/ghc/ghc/-/issues/19411) in GHC 9.0.1.
diff --git a/mod.cabal b/mod.cabal
--- a/mod.cabal
+++ b/mod.cabal
@@ -1,10 +1,10 @@
 name:          mod
-version:       0.1.2.2
+version:       0.2.0.0
 cabal-version: >=1.10
 build-type:    Simple
 license:       MIT
 license-file:  LICENSE
-copyright:     2017-2020 Andrew Lelechenko
+copyright:     2017-2022 Andrew Lelechenko
 maintainer:    Andrew Lelechenko <andrew.lelechenko@gmail.com>
 homepage:      https://github.com/Bodigrim/mod
 bug-reports:   https://github.com/Bodigrim/mod/issues
@@ -12,10 +12,10 @@
 description:
   <https://en.wikipedia.org/wiki/Modular_arithmetic Modular arithmetic>,
   promoting moduli to the type level, with an emphasis on performance.
-  Originally part of <https://hackage.haskell.org/package/arithmoi arithmoi> package.
+  Originally part of the <https://hackage.haskell.org/package/arithmoi arithmoi> package.
 category:      Math, Number Theory
 author:        Andrew Lelechenko <andrew.lelechenko@gmail.com>
-tested-with:   GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.3 GHC ==8.10.4 GHC ==9.0.1
+tested-with:   GHC ==9.0.2 GHC ==9.2.5 GHC ==9.4.4 GHC ==9.6.1
 extra-source-files:
   changelog.md
   README.md
@@ -34,9 +34,9 @@
 
 library
   build-depends:
-    base >=4.10 && <5,
+    base >=4.15 && <5,
     deepseq,
-    integer-gmp <1.2
+    ghc-bignum
   if flag(semirings)
     build-depends:
       semirings >= 0.5
@@ -47,9 +47,15 @@
   exposed-modules:
     Data.Mod
     Data.Mod.Word
+  other-modules:
+    Data.Mod.Compat
   default-language: Haskell2010
-  ghc-options: -Wall -O2 -Wno-deprecations
+  ghc-options: -Wall -O2 -Wno-deprecations -Wcompat
 
+  if(arch(aarch64))
+    c-sources: cbits/aarch64.c
+    include-dirs: cbits
+
 test-suite mod-tests
   build-depends:
     base >=4.10 && <5,
@@ -59,6 +65,7 @@
     tasty-quickcheck >=0.9 && <0.11
   if flag(semirings)
     build-depends:
+      containers,
       quickcheck-classes >=0.6.3,
       semirings >= 0.5
   if flag(vector)
@@ -70,19 +77,19 @@
   main-is: Test.hs
   default-language: Haskell2010
   hs-source-dirs: test
-  ghc-options: -Wall -threaded -rtsopts
+  ghc-options: -Wall -threaded -rtsopts -Wcompat
 
 benchmark mod-bench
   build-depends:
     base,
     mod,
-    -- finite-field,
+    -- finite-field >= 0.9,
     -- finite-typelits,
     -- modular,
-    -- modular-arithmetic,
+    -- modular-arithmetic >= 2,
     tasty-bench >= 0.2.5
   type: exitcode-stdio-1.0
   main-is: Bench.hs
   default-language: Haskell2010
   hs-source-dirs: bench
-  ghc-options: -Wall -O2
+  ghc-options: -Wall -O2 -Wcompat
diff --git a/test/Test.hs b/test/Test.hs
--- a/test/Test.hs
+++ b/test/Test.hs
@@ -2,24 +2,28 @@
 {-# LANGUAGE DataKinds           #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications    #-}
+{-# LANGUAGE TypeOperators       #-}
 
-{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# OPTIONS_GHC -Wno-orphans -Wno-type-defaults #-}
 
 module Main (main) where
 
+import Control.Exception (evaluate, try, ArithException(..))
 import Data.Bits
 import Data.Mod
 import qualified Data.Mod.Word as Word
 import Data.Proxy
 import Data.Semigroup
 import Foreign.Storable (Storable(..))
-import GHC.TypeNats (KnownNat, SomeNat(..), natVal, someNatVal)
+import GHC.TypeNats (KnownNat, SomeNat(..), natVal, someNatVal, type (^), type (+), type (-))
 import Test.Tasty
 import Test.Tasty.QuickCheck
 import Test.QuickCheck.Classes.Base
 
 #ifdef MIN_VERSION_semirings
-import Data.Semiring (Ring)
+import qualified Data.Euclidean as E
+import Data.Semiring (Ring, Semiring(..))
+import qualified Data.Set as S
 import Test.QuickCheck.Classes (semiringLaws, ringLaws)
 #endif
 
@@ -29,76 +33,118 @@
 import Test.QuickCheck.Classes (muvectorLaws, primLaws)
 #endif
 
+#define testModLabeled(lbl, n) \
+  testGroup ("Mod " ++ lbl) $ \
+    testProperty "fromInteger" \
+      (fromIntegerProp (Proxy :: Proxy (n))) : \
+    testProperty "invertMod"   (invertModProp   @(n)) : \
+    testProperty "powMod"      (powModProp      @(n)) : \
+    map lawsToTest (laws (Proxy :: Proxy (Mod (n))))
+
+#define testMod(n) testModLabeled(show (n :: Integer), n)
+
+#define testModWordLabeled(lbl, n) \
+  testGroup ("Word.Mod" ++ lbl) $ \
+    testProperty "fromInteger" \
+      (fromIntegerWordProp (Proxy :: Proxy (n))) : \
+    testProperty "powMod"    (powModWordProp    @(n)) : \
+    testProperty "invertMod" (invertModWordProp @(n)) : \
+    map lawsToTest (laws (Proxy :: Proxy (Word.Mod (n))))
+
+#define testModWord(n) testModWordLabeled(show (n :: Integer), n)
+
 main :: IO ()
-main = defaultMain $ testGroup "All"
+main = defaultMain $ testGroup "All" $
   [ testGroup "Mod 1" $
     testProperty "fromInteger"
       (fromIntegerProp (Proxy :: Proxy 1)) :
     map lawsToTest (laws1 (Proxy :: Proxy (Mod 1)))
-  , testGroup "Mod 2310" $
-    testProperty "fromInteger"
-      (fromIntegerProp (Proxy :: Proxy 2310)) :
-    testProperty "invertMod"   (invertModProp   @2310) :
-    testProperty "powMod"      (powModProp      @2310) :
-    map lawsToTest (laws (Proxy :: Proxy (Mod 2310)))
-  , testGroup "Mod 18446744073709551615" $
-    testProperty "fromInteger"
-      (fromIntegerProp (Proxy :: Proxy 18446744073709551615)) :
-    testProperty "invertMod"   (invertModProp   @18446744073709551615) :
-    testProperty "powMod"      (powModProp      @18446744073709551615) :
-    map lawsToTest (laws (Proxy :: Proxy (Mod 18446744073709551615)))
-  , testGroup "Mod 18446744073709551626" $
-    testProperty "fromInteger"
-      (fromIntegerProp (Proxy :: Proxy 18446744073709551626)) :
-    testProperty "powMod"      (powModProp      @18446744073709551626) :
-    testProperty "invertMod"   (invertModProp   @18446744073709551626) :
-    map lawsToTest (laws (Proxy :: Proxy (Mod 18446744073709551626)))
-  , testGroup "Mod 123456789012345678901234567890" $
-    testProperty "fromInteger"
-      (fromIntegerProp (Proxy :: Proxy 123456789012345678901234567890)) :
-    testProperty "powMod"      (powModProp      @123456789012345678901234567890) :
-    testProperty "invertMod"   (invertModProp   @123456789012345678901234567890) :
-    map lawsToTest (laws (Proxy :: Proxy (Mod 123456789012345678901234567890)))
+  , testMod(2310)
+  , testMod(2^16-1)
+  , testMod(2^16)
+  , testMod(2^16+1)
+  , testMod(2^32-1)
+  , testMod(2^32)
+  , testMod(2^32+1)
+  , testMod(2^64-1)
+  , testMod(2^64)
+  , testMod(2^64+1)
+  , testMod(123456789012345678901234567890)
+  , testModLabeled("2^40000", 2^40000)
   , testGroup "Random Mod"
     [ testProperty "fromInteger" fromIntegerRandomProp
     , testProperty "invertMod"   invertModRandomProp
     , testProperty "powMod"      powModRandomProp
     , testProperty "powMod on sum" powModRandomAdditiveProp
     , testProperty "powMod special case" powModCase
+#ifdef MIN_VERSION_semirings
+    , testProperty "divide"  dividePropRandom
+    , testProperty "gcd"     gcdIsPrincipalIdealRandom
+    , testProperty "lcm"     lcmIsIntersectionOfIdealsRandom
+    , testProperty "coprime" coprimeGeneratorsRandom
+    , testProperty "quotRem" quotRemPropRandom
+    , testProperty "degree"  degreePropRandom
+#endif
     ]
+  , testGroup "Mod 0"
+    [ testProperty "0"            (isDivideByZero 0)
+    , testProperty "1"            (isDivideByZero 1)
+    , testProperty "minBound"     (isDivideByZero minBound)
+    , testProperty "maxBound"     (isDivideByZero maxBound)
+    , testProperty "toEnum"       (isDivideByZero (toEnum 0))
+    , testProperty "fromRational" (isDivideByZero (fromRational 0))
+#ifdef MIN_VERSION_semirings
+    , testProperty "zero"         (isDivideByZero zero)
+    , testProperty "one"          (isDivideByZero one)
+    , testProperty "fromNatural"  (isDivideByZero (fromNatural 0))
+#endif
+    ]
 
   , testGroup "Word.Mod 1" $
     testProperty "fromInteger"
       (fromIntegerWordProp (Proxy :: Proxy 1)) :
     map lawsToTest (laws1 (Proxy :: Proxy (Word.Mod 1)))
-  , testGroup "Word.Mod 2310" $
-    testProperty "fromInteger"
-      (fromIntegerWordProp (Proxy :: Proxy 2310)) :
-    testProperty "powMod"    (powModWordProp    @2310) :
-    testProperty "invertMod" (invertModWordProp @2310) :
-    map lawsToTest (laws (Proxy :: Proxy (Word.Mod 2310)))
-  , if finiteBitSize (0 :: Word) == 64 then
-      testGroup "Word.Mod 18446744073709551615" $
-      testProperty "fromInteger"
-        (fromIntegerWordProp (Proxy :: Proxy 18446744073709551615)) :
-      testProperty "powMod"    (powModWordProp    @18446744073709551615) :
-      testProperty "invertMod" (invertModWordProp @18446744073709551615) :
-      map lawsToTest (laws (Proxy :: Proxy (Word.Mod 18446744073709551615)))
-    else
-      testGroup "Word.Mod 4294967295" $
-      testProperty "fromInteger"
-        (fromIntegerWordProp (Proxy :: Proxy 4294967295)) :
-      testProperty "powMod"    (powModWordProp    @4294967295) :
-      testProperty "invertMod" (invertModWordProp @4294967295) :
-      map lawsToTest (laws (Proxy :: Proxy (Word.Mod 4294967295)))
-  , testGroup "Random Word.Mod"
+  , testMod(2310)
+  , testMod(2^16-1)
+  , testMod(2^16)
+  , testMod(2^16+1)
+  , testMod(2^32-1)
+  ] ++ if finiteBitSize (0 :: Word) /= 64 then [] else
+  [ testMod(2^32)
+  , testMod(2^32+1)
+  , testMod(2^64-1)
+  , testMod(2^64)
+  , testMod(2^64+1)
+  ] ++
+  [ testGroup "Random Word.Mod"
     [ testProperty "fromInteger" fromIntegerWordRandomProp
     , testProperty "invertMod"   invertModWordRandomProp
     , testProperty "invertMod near maxBound" invertModWordRandomPropNearMaxBound
     , testProperty "powMod"      powModWordRandomProp
     , testProperty "powMod on sum" powModWordRandomAdditiveProp
     , testProperty "powMod special case" powModWordCase
+#ifdef MIN_VERSION_semirings
+    , testProperty "divide"  divideWordPropRandom
+    , testProperty "gcd"     gcdIsPrincipalIdealWordRandom
+    , testProperty "lcm"     lcmIsIntersectionOfIdealsWordRandom
+    , testProperty "coprime" coprimeGeneratorsWordRandom
+    , testProperty "quotRem" quotRemWordPropRandom
+    , testProperty "degree"  degreeWordPropRandom
+#endif
     ]
+  , testGroup "Word.Mod 0"
+    [ testProperty "0"            (isDivideByZeroWord 0)
+    , testProperty "1"            (isDivideByZeroWord 1)
+    , testProperty "minBound"     (isDivideByZeroWord minBound)
+    , testProperty "maxBound"     (isDivideByZeroWord maxBound)
+    , testProperty "toEnum"       (isDivideByZeroWord (toEnum 0))
+    , testProperty "fromRational" (isDivideByZeroWord (fromRational 0))
+#ifdef MIN_VERSION_semirings
+    , testProperty "zero"         (isDivideByZeroWord zero)
+    , testProperty "one"          (isDivideByZeroWord one)
+    , testProperty "fromNatural"  (isDivideByZeroWord (fromNatural 0))
+#endif
+    ]
   ]
 
 #ifdef MIN_VERSION_semirings
@@ -267,7 +313,7 @@
   where
     n = 1 `shiftL` 64 :: Integer
 
-newtype Huge a = Huge { getHuge :: a }
+newtype Huge a = Huge { _getHuge :: a }
   deriving (Show)
 
 instance (Bits a, Num a, Arbitrary a) => Arbitrary (Huge a) where
@@ -276,3 +322,139 @@
     ds <- vector l
     return $ Huge $ foldl1 (\acc n -> acc `shiftL` 63 + n) ds
   shrink (Huge n) = Huge <$> shrink n
+
+-------------------------------------------------------------------------------
+-- DivideByZero
+
+isDivideByZero :: Mod 0 -> Property
+isDivideByZero x = ioProperty ((=== Left DivideByZero) <$> try (evaluate x))
+
+isDivideByZeroWord :: Word.Mod 0 -> Property
+isDivideByZeroWord x = ioProperty ((=== Left DivideByZero) <$> try (evaluate x))
+
+-------------------------------------------------------------------------------
+-- Ideals
+
+#ifdef MIN_VERSION_semirings
+
+dividePropRandom :: Positive (Small Integer) -> Positive Integer -> Positive Integer -> Property
+dividePropRandom (Positive (Small m)) (Positive x) (Positive y) = case someNatVal (fromInteger m) of
+  SomeNat (Proxy :: Proxy m) -> divideProp (fromInteger x :: Mod m) (fromInteger y)
+
+divideProp :: KnownNat m => Mod m -> Mod m -> Property
+divideProp x y = case E.divide x y of
+  Just z -> x === y * z
+  Nothing -> filter ((== x) . (* y)) [minBound .. maxBound] === []
+
+gcdIsPrincipalIdealRandom :: Positive (Small Integer) -> Integer -> Integer -> Property
+gcdIsPrincipalIdealRandom (Positive (Small m)) x y = case someNatVal (fromInteger m) of
+  SomeNat (Proxy :: Proxy m) -> gcdIsPrincipalIdeal (fromInteger x :: Mod m) (fromInteger y)
+
+gcdIsPrincipalIdeal :: KnownNat m => Mod m -> Mod m -> Property
+gcdIsPrincipalIdeal x y = addIdeals (genIdeal x) (genIdeal y) === genIdeal (E.gcd x y)
+  where
+    genIdeal t = S.fromList $ map (* t) [minBound .. maxBound]
+    addIdeals us vs = S.fromList [ u + v | u <- S.toList us, v <- S.toList vs ]
+
+lcmIsIntersectionOfIdealsRandom :: Positive (Small Integer) -> Integer -> Integer -> Property
+lcmIsIntersectionOfIdealsRandom (Positive (Small m)) x y = case someNatVal (fromInteger m) of
+  SomeNat (Proxy :: Proxy m) -> lcmIsIntersectionOfIdeals (fromInteger x :: Mod m) (fromInteger y)
+
+lcmIsIntersectionOfIdeals :: KnownNat m => Mod m -> Mod m -> Property
+lcmIsIntersectionOfIdeals x y = S.intersection (genIdeal x) (genIdeal y) === genIdeal (E.lcm x y)
+  where
+    genIdeal t = S.fromList $ map (* t) [minBound .. maxBound]
+
+coprimeGeneratorsRandom :: Positive (Small Integer) -> Integer -> Integer -> Property
+coprimeGeneratorsRandom (Positive (Small m)) x y = case someNatVal (fromInteger m) of
+  SomeNat (Proxy :: Proxy m) -> coprimeGenerators (fromInteger x :: Mod m) (fromInteger y)
+
+coprimeGenerators :: KnownNat m => Mod m -> Mod m -> Property
+coprimeGenerators x y = E.coprime x y === (addIdeals (genIdeal x) (genIdeal y) == S.fromList [minBound .. maxBound])
+  where
+    genIdeal t = S.fromList $ map (* t) [minBound .. maxBound]
+    addIdeals us vs = S.fromList [ u + v | u <- S.toList us, v <- S.toList vs ]
+
+quotRemPropRandom :: Positive (Small Integer) -> Positive Integer -> Positive Integer -> Property
+quotRemPropRandom (Positive (Small m)) (Positive x) (Positive y) = case someNatVal (fromInteger m) of
+  SomeNat (Proxy :: Proxy m) -> quotRemProp (fromInteger x :: Mod m) (fromInteger y)
+
+quotRemProp :: KnownNat m => Mod m -> Mod m -> Property
+quotRemProp x y = case E.divide x y of
+  Just z -> E.quotRem x y === (z, 0)
+  Nothing -> y /= 0 ==> let (q, r) = E.quotRem x y in
+    counterexample (show (q, r)) $ x === q * y + r
+
+degreePropRandom :: Positive (Small Integer) -> Positive Integer -> Positive Integer -> Property
+degreePropRandom (Positive (Small m)) (Positive x) (Positive y) = case someNatVal (fromInteger m) of
+  SomeNat (Proxy :: Proxy m) -> degreeProp (fromInteger x :: Mod m) (fromInteger y)
+
+degreeProp :: KnownNat m => Mod m -> Mod m -> Property
+degreeProp x y = ioProperty $ do
+  ret <- try (evaluate (E.quotRem x y))
+  pure $ case ret of
+    Left DivideByZero -> property True
+    Left{}            -> property False
+    Right (_, r)      -> r === 0 .||. property (E.degree r < E.degree y)
+
+divideWordPropRandom :: Positive Word -> Word -> Word -> Property
+divideWordPropRandom (Positive m) x y = case someNatVal (fromIntegral m) of
+  SomeNat (Proxy :: Proxy m) -> divideWordProp (fromIntegral x :: Word.Mod m) (fromIntegral y)
+
+divideWordProp :: KnownNat m => Word.Mod m -> Word.Mod m -> Property
+divideWordProp x y = case E.divide x y of
+  Just z -> x === y * z
+  Nothing -> filter ((== x) . (* y)) [minBound .. maxBound] === []
+
+gcdIsPrincipalIdealWordRandom :: Positive Word -> Word -> Word -> Property
+gcdIsPrincipalIdealWordRandom (Positive m) x y = case someNatVal (fromIntegral m) of
+  SomeNat (Proxy :: Proxy m) -> gcdIsPrincipalIdealWord (fromIntegral x :: Word.Mod m) (fromIntegral y)
+
+gcdIsPrincipalIdealWord :: KnownNat m => Word.Mod m -> Word.Mod m -> Property
+gcdIsPrincipalIdealWord x y = addIdeals (genIdeal x) (genIdeal y) === genIdeal (E.gcd x y)
+  where
+    genIdeal t = S.fromList $ map (* t) [minBound .. maxBound]
+    addIdeals us vs = S.fromList [ u + v | u <- S.toList us, v <- S.toList vs ]
+
+lcmIsIntersectionOfIdealsWordRandom :: Positive Word -> Word -> Word -> Property
+lcmIsIntersectionOfIdealsWordRandom (Positive m) x y = case someNatVal (fromIntegral m) of
+  SomeNat (Proxy :: Proxy m) -> lcmIsIntersectionOfIdealsWord (fromIntegral x :: Word.Mod m) (fromIntegral y)
+
+lcmIsIntersectionOfIdealsWord :: KnownNat m => Word.Mod m -> Word.Mod m -> Property
+lcmIsIntersectionOfIdealsWord x y = S.intersection (genIdeal x) (genIdeal y) === genIdeal (E.lcm x y)
+  where
+    genIdeal t = S.fromList $ map (* t) [minBound .. maxBound]
+
+coprimeGeneratorsWordRandom :: Positive Word -> Word -> Word -> Property
+coprimeGeneratorsWordRandom (Positive m) x y = case someNatVal (fromIntegral m) of
+  SomeNat (Proxy :: Proxy m) -> coprimeGeneratorsWord (fromIntegral x :: Word.Mod m) (fromIntegral y)
+
+coprimeGeneratorsWord :: KnownNat m => Word.Mod m -> Word.Mod m -> Property
+coprimeGeneratorsWord x y = E.coprime x y === (addIdeals (genIdeal x) (genIdeal y) == S.fromList [minBound .. maxBound])
+  where
+    genIdeal t = S.fromList $ map (* t) [minBound .. maxBound]
+    addIdeals us vs = S.fromList [ u + v | u <- S.toList us, v <- S.toList vs ]
+
+quotRemWordPropRandom :: Positive Word -> Word -> Word -> Property
+quotRemWordPropRandom (Positive m) x y = case someNatVal (fromIntegral m) of
+  SomeNat (Proxy :: Proxy m) -> quotRemWordProp (fromIntegral x :: Word.Mod m) (fromIntegral y)
+
+quotRemWordProp :: KnownNat m => Word.Mod m -> Word.Mod m -> Property
+quotRemWordProp x y = case E.divide x y of
+  Just z -> E.quotRem x y === (z, 0)
+  Nothing -> y /= 0 ==> let (q, r) = E.quotRem x y in
+    counterexample (show (q, r)) $ x === q * y + r
+
+degreeWordPropRandom :: Positive Word -> Word -> Word -> Property
+degreeWordPropRandom (Positive m) x y = case someNatVal (fromIntegral m) of
+  SomeNat (Proxy :: Proxy m) -> degreeWordProp (fromIntegral x :: Word.Mod m) (fromIntegral y)
+
+degreeWordProp :: KnownNat m => Word.Mod m -> Word.Mod m -> Property
+degreeWordProp x y = ioProperty $ do
+  ret <- try (evaluate (E.quotRem x y))
+  pure $ case ret of
+    Left DivideByZero -> property True
+    Left{}            -> property False
+    Right (_, r)      -> r === 0 .||. property (E.degree r < E.degree y)
+
+#endif