diff --git a/Math/NumberTheory/ArithmeticFunctions/Inverse.hs b/Math/NumberTheory/ArithmeticFunctions/Inverse.hs
index 6f5bfe7b4..b46143219 100644
--- a/Math/NumberTheory/ArithmeticFunctions/Inverse.hs
+++ b/Math/NumberTheory/ArithmeticFunctions/Inverse.hs
@@ -16,7 +16,9 @@
 
 module Math.NumberTheory.ArithmeticFunctions.Inverse
   ( inverseTotient
+  , inverseJordan
   , inverseSigma
+  , inverseSigmaK
   , -- * Wrappers
     MinWord(..)
   , MaxWord(..)
@@ -44,7 +46,7 @@ import Numeric.Natural
 
 import Math.NumberTheory.ArithmeticFunctions
 import Math.NumberTheory.Logarithms
-import Math.NumberTheory.Roots (integerRoot)
+import Math.NumberTheory.Roots (exactRoot, integerRoot)
 import Math.NumberTheory.Primes
 import Math.NumberTheory.Utils.DirichletSeries (DirichletSeries)
 import qualified Math.NumberTheory.Utils.DirichletSeries as DS
@@ -70,33 +72,37 @@ atomicSeries point (ArithmeticFunction f g) (PrimePowers p ks) =
   DS.fromDistinctAscList (map (\k -> (g (f p k), point (unPrime p ^ k))) ks)
 
 -- | See section 5.1 of the paper.
-invTotient
-  :: forall a. (UniqueFactorisation a, Eq a)
-  => [(Prime a, Word)]
+invJordan
+  :: forall a. (Integral a, UniqueFactorisation a, Eq a)
+  => Word
+  -- ^ Value of 'k' in 'jordan_k'
+  -> [(Prime a, Word)]
   -- ^ Factorisation of a value of the totient function
   -> [PrimePowers a]
   -- ^ Possible prime factors of an argument of the totient function
-invTotient fs = map (\p -> PrimePowers p (doPrime p)) ps
+invJordan k fs = map (\p -> PrimePowers p (doPrime p)) ps
   where
     divs :: [a]
     divs = runFunctionOnFactors divisorsListA fs
 
     ps :: [Prime a]
-    ps = mapMaybe (isPrime . (+ 1)) divs
+    ps = mapMaybe (\d -> exactRoot k (d + 1) >>= isPrime) divs
 
     doPrime :: Prime a -> [Word]
     doPrime p = case lookup p fs of
       Nothing -> [1]
-      Just k  -> [1 .. k+1]
+      Just w  -> [1 .. w+1]
 
 -- | See section 5.2 of the paper.
 invSigma
   :: forall a. (Euclidean a, Integral a, UniqueFactorisation a, Enum (Prime a), Bits a)
-  => [(Prime a, Word)]
+  => Word
+  -- ^ Value of 'k' in 'sigma_k'
+  -> [(Prime a, Word)]
   -- ^ Factorisation of a value of the sum-of-divisors function
   -> [PrimePowers a]
   -- ^ Possible prime factors of an argument of the sum-of-divisors function
-invSigma fs
+invSigma k fs
   = map (\(x, ys) -> PrimePowers x (S.toList ys))
   $ M.assocs
   $ M.unionWith (<>) pksSmall pksLarge
@@ -108,7 +114,7 @@ invSigma fs
     divs = runFunctionOnFactors divisorsListA fs
 
     n :: a
-    n = product $ map (\(p, k) -> unPrime p ^ k) fs
+    n = factorBack fs
 
     -- There are two possible strategies to find possible prime factors
     -- of an argument of the sum-of-divisors function.
@@ -141,17 +147,17 @@ invSigma fs
     doPrime :: Prime a -> Set Word
     doPrime p' = let p = unPrime p' in S.fromDistinctAscList
       [ e
-      | e <- [1 .. intToWord (integerLogBase (toInteger p) (toInteger n))]
-      , n `rem` ((p ^ (e + 1) - 1) `quot` (p - 1)) == 0
+      | e <- [1 .. intToWord (integerLogBase (toInteger (p ^ k)) (toInteger n))]
+      , n `rem` ((p ^ (k * (e + 1)) - 1) `quot` (p ^ k - 1)) == 0
       ]
 
     pksLarge :: Map (Prime a) (Set Word)
     pksLarge = M.unionsWith (<>)
       [ maybe mempty (flip M.singleton (S.singleton e)) (isPrime p)
       | d <- divs
-      , e <- [1 .. intToWord (integerLogBase (toInteger lim) (toInteger d))]
-      , let p = integerRoot e (d - 1)
-      , p ^ (e + 1) - 1 == d * (p - 1)
+      , e <- [1 .. intToWord (quot (integerLogBase (toInteger lim) (toInteger d)) (wordToInt k)) ]
+      , let p = integerRoot (e * k) (d - 1)
+      , p ^ (k * (e + 1)) - 1 == d * (p ^ k - 1)
       ]
 
 -- | Instead of multiplying all atomic series and filtering out everything,
@@ -264,16 +270,55 @@ invertFunction point f invF n
 -- >>> unMaxWord (inverseTotient MaxWord 120)
 -- 462
 inverseTotient
-  :: (Semiring b, Euclidean a, UniqueFactorisation a, Ord a)
+  :: (Semiring b, Integral a, Euclidean a, UniqueFactorisation a)
   => (a -> b)
   -> a
   -> b
-inverseTotient point = invertFunction point totientA invTotient
+inverseTotient = inverseJordan 1
 {-# SPECIALISE inverseTotient :: Semiring b => (Int -> b) -> Int -> b #-}
 {-# SPECIALISE inverseTotient :: Semiring b => (Word -> b) -> Word -> b #-}
 {-# SPECIALISE inverseTotient :: Semiring b => (Integer -> b) -> Integer -> b #-}
 {-# SPECIALISE inverseTotient :: Semiring b => (Natural -> b) -> Natural -> b #-}
 
+-- | The inverse for 'jordan_k' function, where 'k' is a nonnegative integer.
+--
+-- Generalizes the 'inverseTotient' function, which is 'inverseJordan' when 'k'
+-- is '1'.
+--
+-- The return value is parameterized by a 'Semiring', which allows
+-- various applications by providing different (multiplicative) embeddings.
+-- E. g., list all preimages (see a helper 'asSetOfPreimages'):
+--
+-- >>> import qualified Data.Set as S
+-- >>> import Data.Semigroup
+-- >>> S.mapMonotonic getProduct (inverseJordan 2 (S.singleton . Product) 192)
+-- fromList [15,16]
+--
+-- Similarly to 'inverseTotient', it is possible to count and sum preimages, or
+-- get the maximum/minimum preimage.
+--
+-- Note: it is the __user's responsibility__ to use an appropriate type for
+-- 'inverseSigmaK'. Even low values of 'k' with 'Int' or 'Word' will return
+-- invalid results due to over/underflow, or throw exceptions (i.e. division by
+-- zero).
+--
+-- >>> asSetOfPreimages (inverseJordan 15) (jordan 15 19) :: S.Set Int
+-- fromList []
+--
+-- >>> asSetOfPreimages (inverseJordan 15) (jordan 15 19) :: S.Set Integer
+-- fromList [19]
+inverseJordan
+  :: (Semiring b, Integral a, Euclidean a, UniqueFactorisation a)
+  => Word
+  -> (a -> b)
+  -> a
+  -> b
+inverseJordan k point = invertFunction point (jordanA k) (invJordan k)
+{-# SPECIALISE inverseJordan :: Semiring b => Word -> (Int -> b) -> Int -> b #-}
+{-# SPECIALISE inverseJordan :: Semiring b => Word -> (Word -> b) -> Word -> b #-}
+{-# SPECIALISE inverseJordan :: Semiring b => Word -> (Integer -> b) -> Integer -> b #-}
+{-# SPECIALISE inverseJordan :: Semiring b => Word -> (Natural -> b) -> Natural -> b #-}
+
 -- | The inverse for 'sigma' 1 function.
 --
 -- The return value is parameterized by a 'Semiring', which allows
@@ -306,12 +351,48 @@ inverseSigma
   => (a -> b)
   -> a
   -> b
-inverseSigma point = invertFunction point (sigmaA 1) invSigma
+inverseSigma = inverseSigmaK 1
 {-# SPECIALISE inverseSigma :: Semiring b => (Int -> b) -> Int -> b #-}
 {-# SPECIALISE inverseSigma :: Semiring b => (Word -> b) -> Word -> b #-}
 {-# SPECIALISE inverseSigma :: Semiring b => (Integer -> b) -> Integer -> b #-}
 {-# SPECIALISE inverseSigma :: Semiring b => (Natural -> b) -> Natural -> b #-}
 
+-- | The inverse for 'sigma_k' function, where 'k' is a nonnegative integer.
+--
+-- Generalizes the 'inverseSigma' function, which is 'inverseSigmaK' when 'k'
+-- is '1'.
+--
+-- The return value is parameterized by a 'Semiring', which allows
+-- various applications by providing different (multiplicative) embeddings.
+-- E. g., list all preimages (see a helper 'asSetOfPreimages'):
+--
+-- >>> import qualified Data.Set as S
+-- >>> import Data.Semigroup
+-- >>> S.mapMonotonic getProduct (inverseSigmaK 2 (S.singleton . Product) 850)
+-- fromList [24, 26]
+--
+-- Similarly to 'inverseSigma', it is possible to count and sum preimages, or
+-- get the maximum/minimum preimage.
+--
+-- Note: it is the __user's responsibility__ to use an appropriate type for
+-- 'inverseSigmaK'. Even low values of 'k' with 'Int' or 'Word' will return
+-- invalid results due to over/underflow, or throw exceptions (i.e. division by
+-- zero).
+--
+-- >>> asSetOfPreimages (inverseSigmaK 17) (sigma 17 13) :: S.Set Int
+-- fromList *** Exception: divide by zero
+inverseSigmaK
+  :: (Semiring b, Euclidean a, UniqueFactorisation a, Integral a, Enum (Prime a), Bits a)
+  => Word
+  -> (a -> b)
+  -> a
+  -> b
+inverseSigmaK k point = invertFunction point (sigmaA k) (invSigma k)
+{-# SPECIALISE inverseSigmaK :: Semiring b => Word -> (Int -> b) -> Int -> b #-}
+{-# SPECIALISE inverseSigmaK :: Semiring b => Word -> (Word -> b) -> Word -> b #-}
+{-# SPECIALISE inverseSigmaK :: Semiring b => Word -> (Integer -> b) -> Integer -> b #-}
+{-# SPECIALISE inverseSigmaK :: Semiring b => Word -> (Natural -> b) -> Natural -> b #-}
+
 --------------------------------------------------------------------------------
 -- Wrappers
 
diff --git a/Math/NumberTheory/Quadratic/EisensteinIntegers.hs b/Math/NumberTheory/Quadratic/EisensteinIntegers.hs
index 071a1de79..9e0f2308a 100644
--- a/Math/NumberTheory/Quadratic/EisensteinIntegers.hs
+++ b/Math/NumberTheory/Quadratic/EisensteinIntegers.hs
@@ -303,8 +303,7 @@ divideByPrime p p' np k = go k 0
             where
                 (d1, z') = go1 c 0 z
                 d2 = c - d1
-                z'' = head $ drop (wordToInt d2)
-                    $ iterate (\g -> fromMaybe err $ (g * unPrime p) `quotEvenI` np) z'
+                z'' = iterate (\g -> fromMaybe err $ (g * unPrime p) `quotEvenI` np) z' !! max 0 (wordToInt d2)
 
         go1 :: Word -> Word -> EisensteinInteger -> (Word, EisensteinInteger)
         go1 0 d z = (d, z)
diff --git a/Math/NumberTheory/Zeta/Hurwitz.hs b/Math/NumberTheory/Zeta/Hurwitz.hs
index bf0f798a1..bca1d822c 100644
--- a/Math/NumberTheory/Zeta/Hurwitz.hs
+++ b/Math/NumberTheory/Zeta/Hurwitz.hs
@@ -75,7 +75,7 @@ zetaHurwitz eps a = zipWith3 (\s i t -> s + i + t) ss is ts
                       (\powOfA int -> powOfA * fromInteger int)
                       powsOfAPlusN
                       [-1, 0..]
-         in zipWith (/) (repeat aPlusN) denoms
+         in map ((/) aPlusN) denoms
 
     -- [      1      |             ]
     -- [ ----------- | s <- [0 ..] ]
@@ -109,9 +109,8 @@ zetaHurwitz eps a = zipWith3 (\s i t -> s + i + t) ss is ts
         (skipEvens powsOfAPlusN)
 
     fracs :: [a]
-    fracs = zipWith
-            (\sec pochh -> sum $ zipWith (\s p -> s * fromInteger p) sec pochh)
-            (repeat second)
+    fracs = map
+            (\pochh -> sum $ zipWith (\s p -> s * fromInteger p) second pochh)
             pochhammers
 
     -- Infinite list of @T@ values in 4.8.5 formula, for every @s@ in
diff --git a/benchmark/Math/NumberTheory/InverseBench.hs b/benchmark/Math/NumberTheory/InverseBench.hs
index 70c0ae2c7..61b747a15 100644
--- a/benchmark/Math/NumberTheory/InverseBench.hs
+++ b/benchmark/Math/NumberTheory/InverseBench.hs
@@ -21,7 +21,7 @@ fact = product [1..13]
 tens :: Num a => a
 tens = 10 ^ 18
 
-countInverseTotient :: (Ord a, Euclidean a, UniqueFactorisation a) => a -> Word
+countInverseTotient :: (Ord a, Integral a, Euclidean a, UniqueFactorisation a) => a -> Word
 countInverseTotient = inverseTotient (const 1)
 
 countInverseSigma :: (Integral a, Euclidean a, UniqueFactorisation a, Enum (Prime a), Bits a) => a -> Word
diff --git a/test-suite/Math/NumberTheory/ArithmeticFunctions/InverseTests.hs b/test-suite/Math/NumberTheory/ArithmeticFunctions/InverseTests.hs
index 623f7d4c9..da11d0c27 100644
--- a/test-suite/Math/NumberTheory/ArithmeticFunctions/InverseTests.hs
+++ b/test-suite/Math/NumberTheory/ArithmeticFunctions/InverseTests.hs
@@ -9,8 +9,11 @@
 --
 
 {-# LANGUAGE FlexibleContexts      #-}
+{-# LANGUAGE RankNTypes            #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 
+{-# OPTIONS_GHC -fconstraint-solver-iterations=0 #-}
+
 {-# OPTIONS_GHC -fno-warn-type-defaults #-}
 
 module Math.NumberTheory.ArithmeticFunctions.InverseTests
@@ -19,16 +22,21 @@ module Math.NumberTheory.ArithmeticFunctions.InverseTests
 
 import Test.Tasty
 import Test.Tasty.HUnit
+import Test.Tasty.SmallCheck as SC hiding (test)
+import Test.Tasty.QuickCheck as QC hiding (Positive)
 
 import Data.Bits (Bits)
 import Data.Euclidean
+import Data.Semiring (Semiring)
 import qualified Data.Set as S
+import Numeric.Natural (Natural)
 
 import Math.NumberTheory.ArithmeticFunctions
 import Math.NumberTheory.ArithmeticFunctions.Inverse
 import Math.NumberTheory.Primes
 import Math.NumberTheory.Recurrences
 import Math.NumberTheory.TestUtils
+import Math.NumberTheory.TestUtils.Wrappers (Power (..))
 
 -------------------------------------------------------------------------------
 -- Totient
@@ -36,9 +44,28 @@ import Math.NumberTheory.TestUtils
 totientProperty1 :: forall a. (Euclidean a, Integral a, UniqueFactorisation a) => Positive a -> Bool
 totientProperty1 (Positive x) = x `S.member` asSetOfPreimages inverseTotient (totient x)
 
+jordanProperty1
+  :: (Euclidean a, Integral a, UniqueFactorisation a)
+  => Power Word
+  -> Positive a
+  -> Bool
+jordanProperty1 (Power k') (Positive x) =
+  -- 'k' shouldn't be large to avoid slow tests.
+  let k = 2 + k' `Prelude.mod` 20
+  in x `S.member` asSetOfPreimages (inverseJordan k) (jordan k x)
+
 totientProperty2 :: (Euclidean a, Integral a, UniqueFactorisation a) => Positive a -> Bool
 totientProperty2 (Positive x) = all (== x) (S.map totient (asSetOfPreimages inverseTotient x))
 
+jordanProperty2
+  :: (Euclidean a, Integral a, UniqueFactorisation a, Ord a)
+  => Power Word
+  -> Positive a
+  -> Bool
+jordanProperty2 (Power k') (Positive x) =
+  let k = 2 + k' `Prelude.mod` 20
+  in all (== x) (S.map (jordan k) (asSetOfPreimages (inverseJordan k) x))
+
 -- | http://oeis.org/A055506
 totientCountFactorial :: [Word]
 totientCountFactorial =
@@ -132,15 +159,74 @@ totientSpecialCases3 = zipWith mkAssert (tail factorial) totientMaxFactorial
     totientMax :: Word -> Word
     totientMax = unMaxWord . inverseTotient MaxWord
 
+jordans5 :: [Word]
+jordans5 =
+  [ 1
+  , 31
+  , 242
+  , 992
+  , 3124
+  , 7502
+  , 16806
+  , 31744
+  , 58806
+  , 96844
+  , 161050
+  , 240064
+  , 371292
+  , 520986
+  , 756008
+  , 1015808
+  , 1419856
+  , 1822986
+  , 2476098
+  , 3099008
+  , 4067052
+  , 4992550
+  , 6436342
+  , 7682048
+  , 9762500
+  , 11510052
+  , 14289858
+  , 16671552
+  , 20511148
+  ]
+
+jordanSpecialCase1 :: [Assertion]
+jordanSpecialCase1 = zipWith mkAssert ixs jordans5
+  where
+    mkAssert a b = assertEqual "should be equal" (S.singleton a) (asSetOfPreimages (inverseJordan 5) b)
+    ixs = [1 .. 29]
+
 -------------------------------------------------------------------------------
 -- Sigma
 
 sigmaProperty1 :: forall a. (Euclidean a, UniqueFactorisation a, Integral a, Enum (Prime a), Bits a) => Positive a -> Bool
 sigmaProperty1 (Positive x) = x `S.member` asSetOfPreimages inverseSigma (sigma 1 x)
 
+sigmaKProperty1
+  :: forall a
+   . (Euclidean a, UniqueFactorisation a, Integral a, Enum (Prime a), Bits a)
+  => Power Word
+  -> Positive a
+  -> Bool
+sigmaKProperty1 (Power k') (Positive x) =
+  -- 'k' shouldn't be large to avoid slow tests.
+  let k = 2 + k' `Prelude.mod` 20
+  in x `S.member` asSetOfPreimages (inverseSigmaK k) (sigma k x)
+
 sigmaProperty2 :: (Euclidean a, UniqueFactorisation a, Integral a, Enum (Prime a), Bits a) => Positive a -> Bool
 sigmaProperty2 (Positive x) = all (== x) (S.map (sigma 1) (asSetOfPreimages inverseSigma x))
 
+sigmaKProperty2
+  :: (Euclidean a, UniqueFactorisation a, Integral a, Enum (Prime a), Bits a)
+  => Power Word
+  -> Positive a
+  -> Bool
+sigmaKProperty2 (Power k') (Positive x) =
+  let k = 2 + k' `Prelude.mod` 20
+  in all (== x) (S.map (sigma k) (asSetOfPreimages (inverseSigmaK k) x))
+
 -- | http://oeis.org/A055486
 sigmaCountFactorial :: [Word]
 sigmaCountFactorial =
@@ -234,9 +320,60 @@ sigmaSpecialCase4 :: Assertion
 sigmaSpecialCase4 = assertBool "200 should be in inverseSigma(sigma(200))" $
   sigmaProperty1 $ Positive (200 :: Word)
 
+sigmas5 :: [Word]
+sigmas5 =
+  [ 1
+  , 33
+  , 244
+  , 1057
+  , 3126
+  , 8052
+  , 16808
+  , 33825
+  , 59293
+  , 103158
+  , 161052
+  , 257908
+  , 371294
+  , 554664
+  , 762744
+  , 1082401
+  , 1419858
+  , 1956669
+  , 2476100
+  , 3304182
+  , 4101152
+  , 5314716
+  , 6436344
+  , 8253300
+  , 9768751
+  , 12252702
+  , 14408200
+  , 17766056
+  , 20511150
+  ]
+
+sigmaSpecialCase5 :: [Assertion]
+sigmaSpecialCase5 = zipWith mkAssert ixs sigmas5
+ where
+  mkAssert a b = assertEqual "should be equal" (S.singleton a) (asSetOfPreimages (inverseSigmaK 5) b)
+  ixs = [1 .. 29]
+
 -------------------------------------------------------------------------------
 -- TestTree
 
+-- Tests for 'Int', 'Word' are omitted because 'inverseSigmaK/inverseJordan'
+-- tests would quickly oveflow in these types.
+testIntegralPropertyNoLargeInverse
+  :: forall bool. (SC.Testable IO bool, QC.Testable bool)
+  => String -> (forall a. (Euclidean a, Semiring a, Integral a, Bits a, UniqueFactorisation a, Show a, Enum (Prime a)) => Power Word -> Positive a -> bool) -> TestTree
+testIntegralPropertyNoLargeInverse name f = testGroup name
+  [ SC.testProperty "smallcheck Integer" (f :: Power Word -> Positive Integer -> bool)
+  , SC.testProperty "smallcheck Natural" (f :: Power Word -> Positive Natural -> bool)
+  , QC.testProperty "quickcheck Integer" (f :: Power Word -> Positive Integer -> bool)
+  , QC.testProperty "quickcheck Natural" (f :: Power Word -> Positive Natural -> bool)
+  ]
+
 testSuite :: TestTree
 testSuite = testGroup "Inverse"
   [ testGroup "Totient"
@@ -260,4 +397,18 @@ testSuite = testGroup "Inverse"
     , testGroup "max"
       (zipWith (\i a -> testCase ("factorial " ++ show i) a) [1..] sigmaSpecialCases3)
     ]
+
+  , testGroup "Jordan"
+    [ testIntegralPropertyNoLargeInverse "forward"  jordanProperty1
+    , testIntegralPropertyNoLargeInverse "backward" jordanProperty2
+    , testGroup "inverseJordan"
+      (zipWith (\i test -> testCase ("inverseJordan 5" ++ show i) test) jordans5 jordanSpecialCase1)
+    ]
+
+  , testGroup  "SigmaK"
+    [ testIntegralPropertyNoLargeInverse "forward"  sigmaKProperty1
+    , testIntegralPropertyNoLargeInverse "backward" sigmaKProperty2
+    , testGroup "inverseSigma"
+      (zipWith (\i test -> testCase ("inverseSigma 5" ++ show i) test) sigmas5 sigmaSpecialCase5)
+    ]
   ]
diff --git a/test-suite/Math/NumberTheory/TestUtils.hs b/test-suite/Math/NumberTheory/TestUtils.hs
index 182a78f9e..4a0045561 100644
--- a/test-suite/Math/NumberTheory/TestUtils.hs
+++ b/test-suite/Math/NumberTheory/TestUtils.hs
@@ -12,7 +12,6 @@
 {-# LANGUAGE DataKinds                  #-}
 {-# LANGUAGE FlexibleContexts           #-}
 {-# LANGUAGE FlexibleInstances          #-}
-{-# LANGUAGE KindSignatures             #-}
 {-# LANGUAGE MultiParamTypeClasses      #-}
 {-# LANGUAGE RankNTypes                 #-}
 {-# LANGUAGE ScopedTypeVariables        #-}
@@ -43,6 +42,9 @@ module Math.NumberTheory.TestUtils
   -- * Export for @Zeta@ tests
   , assertEqualUpToEps
 
+  -- * Export for Inverse tests
+  , TestableIntegral
+
   , lawsToTest
   ) where