diff --git a/Crypto/Number/F2m.hs b/Crypto/Number/F2m.hs
index 93b1f48..6ca2604 100644
--- a/Crypto/Number/F2m.hs
+++ b/Crypto/Number/F2m.hs
@@ -16,7 +16,9 @@ module Crypto.Number.F2m
     , mulF2m
     , squareF2m'
     , squareF2m
+    , powF2m
     , modF2m
+    , sqrtF2m
     , invF2m
     , divF2m
     ) where
@@ -66,8 +68,8 @@ mulF2m :: BinaryPolynomial -- ^ Modulus
 mulF2m fx n1 n2
     |    fx < 0
       || n1 < 0
-      || n2 < 0 = error "mulF2m: negative number represent no binary binary polynomial"
-    | fx == 0   = error "modF2m: cannot multiply modulo zero polynomial"
+      || n2 < 0 = error "mulF2m: negative number represent no binary polynomial"
+    | fx == 0   = error "mulF2m: cannot multiply modulo zero polynomial"
     | otherwise = modF2m fx $ go (if n2 `mod` 2 == 1 then n1 else 0) (log2 n2)
       where
         go n s | s == 0  = n
@@ -96,10 +98,37 @@ squareF2m fx = modF2m fx . squareF2m'
 squareF2m' :: Integer
            -> Integer
 squareF2m' n
-    | n < 0     = error "mulF2m: negative number represent no binary binary polynomial"
+    | n < 0     = error "mulF2m: negative number represent no binary polynomial"
     | otherwise = foldl' (\acc s -> if testBit n s then setBit acc (2 * s) else acc) 0 [0 .. log2 n]
 {-# INLINE squareF2m' #-}
 
+-- | Exponentiation in F₂m by computing @a^b mod fx@.
+--
+-- This implements an exponentiation by squaring based solution. It inherits the
+-- same restrictions as 'squareF2m'. Negative exponents are disallowed.
+powF2m :: BinaryPolynomial -- ^Modulus
+       -> Integer          -- ^a
+       -> Integer          -- ^b
+       -> Integer
+powF2m fx a b
+  | b < 0     = error "powF2m: negative exponents disallowed"
+  | b == 0    = if fx > 1 then 1 else 0
+  | even b    = squareF2m fx x
+  | otherwise = mulF2m fx a (squareF2m' x)
+  where x = powF2m fx a (b `div` 2)
+
+-- | Square rooot in F₂m.
+--
+-- We exploit the fact that @a^(2^m) = a@, or in particular, @a^(2^m - 1) = 1@
+-- from a classical result by Lagrange. Thus the square root is simply @a^(2^(m
+-- - 1))@.
+sqrtF2m :: BinaryPolynomial -- ^Modulus
+        -> Integer          -- ^a
+        -> Integer
+sqrtF2m fx a = go (log2 fx - 1) a
+  where go 0 x = x
+        go n x = go (n - 1) (squareF2m fx x)
+
 -- | Extended GCD algorithm for polynomials. For @a@ and @b@ returns @(g, u, v)@ such that @a * u + b * v == g@.
 --
 -- Reference: https://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#B.C3.A9zout.27s_identity_and_extended_GCD_algorithm
diff --git a/tests/Number/F2m.hs b/tests/Number/F2m.hs
index afa6e50..19acd61 100644
--- a/tests/Number/F2m.hs
+++ b/tests/Number/F2m.hs
@@ -2,6 +2,7 @@ module Number.F2m (tests) where
 
 import Imports hiding ((.&.))
 import Data.Bits
+import Data.Maybe
 import Crypto.Number.Basic (log2)
 import Crypto.Number.F2m
 
@@ -52,6 +53,32 @@ mulTests = testGroup "mulF2m"
 squareTests = testGroup "squareF2m"
     [ testProperty "sqr(a) == a * a"
         $ \(Positive m) (NonNegative a) -> mulF2m m a a == squareF2m m a
+    -- disabled because we require @m@ to be a suitable modulus and there is no
+    -- way to guarantee this
+    -- , testProperty "sqrt(a) * sqrt(a) = a"
+    --     $ \(Positive m) (NonNegative aa) -> let a = sqrtF2m m aa in mulF2m m a a == modF2m m aa
+    , testProperty "sqrt(a) * sqrt(a) = a in GF(2^16)"
+        $ let m = 65581 :: Integer -- x^16 + x^5 + x^3 + x^2 + 1
+              nums = [0 .. 65535 :: Integer]
+          in  nums == [let y = sqrtF2m m x in squareF2m m y | x <- nums]
+    ]
+
+powTests = testGroup "powF2m"
+    [ testProperty "2 is square"
+        $ \(Positive m) (NonNegative a) -> powF2m m a 2 == squareF2m m a
+    , testProperty "1 is identity"
+        $ \(Positive m) (NonNegative a) -> powF2m m a 1 == modF2m m a
+    , testProperty "0 is annihilator"
+        $ \(Positive m) (NonNegative a) -> powF2m m a 0 == modF2m m 1
+    , testProperty "(a * b) ^ c == (a ^ c) * (b ^ c)"
+        $ \(Positive m) (NonNegative a) (NonNegative b) (NonNegative c)
+            -> powF2m m (mulF2m m a b) c == mulF2m m (powF2m m a c) (powF2m m b c)
+    , testProperty "a ^ (b + c) == (a ^ b) * (a ^ c)"
+        $ \(Positive m) (NonNegative a) (NonNegative b) (NonNegative c)
+            -> powF2m m a (b + c) == mulF2m m (powF2m m a b) (powF2m m a c)
+    , testProperty "a ^ (b * c) == (a ^ b) ^ c"
+        $ \(Positive m) (NonNegative a) (NonNegative b) (NonNegative c)
+            -> powF2m m a (b * c) == powF2m m (powF2m m a b) c
     ]
 
 invTests = testGroup "invF2m"
@@ -70,7 +97,7 @@ divTests = testGroup "divF2m"
             -> divF2m m a b == (mulF2m m a <$> invF2m m b)
     , testProperty "a * b / b == a"
         $ \(Positive m) (NonNegative a) (NonNegative b)
-            -> invF2m m b == Nothing || divF2m m (mulF2m m a b) b == Just (modF2m m a)
+            -> isNothing (invF2m m b) || divF2m m (mulF2m m a b) b == Just (modF2m m a)
     ]
 
 tests = testGroup "number.F2m"
@@ -78,6 +105,7 @@ tests = testGroup "number.F2m"
     , modTests
     , mulTests
     , squareTests
+    , powTests
     , invTests
     , divTests
     ]