diff --git a/Crypto/Number/F2m.hs b/Crypto/Number/F2m.hs
index 274b8df..503c309 100644
--- a/Crypto/Number/F2m.hs
+++ b/Crypto/Number/F2m.hs
@@ -9,100 +9,133 @@
 -- not optimal and it doesn't provide protection against timing
 -- attacks. The 'm' parameter is implicitly derived from the irreducible
 -- polynomial where applicable.
+
 module Crypto.Number.F2m
     ( BinaryPolynomial
     , addF2m
     , mulF2m
+    , squareF2m'
     , squareF2m
     , modF2m
     , invF2m
     , divF2m
     ) where
 
-import Data.Bits ((.&.),(.|.),xor,shift,testBit)
-import Crypto.Number.Basic
+import Data.Bits (xor, shift, testBit, setBit)
+import Data.List
 import Crypto.Internal.Imports
+import Crypto.Number.Basic
 
 -- | Binary Polynomial represented by an integer
 type BinaryPolynomial = Integer
 
--- | Addition over F₂m. This is just a synonym of  'xor'.
-addF2m :: Integer -> Integer -> Integer
+-- | Addition over F₂m. This is just a synonym of 'xor'.
+addF2m :: Integer
+       -> Integer
+       -> Integer
 addF2m = xor
 {-# INLINE addF2m #-}
 
--- | Binary polynomial reduction modulo using long division algorithm.
-modF2m :: BinaryPolynomial -- ^ Irreducible binary polynomial
-       -> Integer -> Integer
-modF2m fx = go
-  where
-    lfx = log2 fx
-    go n | s == 0  = n `xor` fx
-         | s < 0 = n
-         | otherwise = go $ n `xor` shift fx s
+-- | Reduction by modulo over F₂m.
+--
+-- This function is undefined for negative arguments, because their bit
+-- representation is platform-dependent. Zero modulus is also prohibited.
+modF2m :: BinaryPolynomial -- ^ Modulus
+       -> Integer
+       -> Integer
+modF2m fx i
+    | fx < 0 || i < 0 = error "modF2m: negative number represent no binary polynomial"
+    | fx == 0         = error "modF2m: cannot divide by zero polynomial"
+    | fx == 1         = 0
+    | otherwise       = go i
       where
-        s = log2 n - lfx
+        lfx = log2 fx
+        go n | s == 0    = n `addF2m` fx
+             | s < 0     = n
+             | otherwise = go $ n `addF2m` shift fx s
+                where s = log2 n - lfx
 {-# INLINE modF2m #-}
 
 -- | Multiplication over F₂m.
 --
---     n1 * n2 (in F(2^m))
-mulF2m :: BinaryPolynomial  -- ^ Irreducible binary polynomial
-       -> Integer -> Integer -> Integer
-mulF2m fx n1 n2 = modF2m fx
-                $ go (if n2 `mod` 2 == 1 then n1 else 0) (log2 n2)
-  where
-    go n s | s == 0  = n
-           | otherwise = if testBit n2 s
-                            then go (n `xor` shift n1 s) (s - 1)
-                            else go n (s - 1)
+-- This function is undefined for negative arguments, because their bit
+-- representation is platform-dependent. Zero modulus is also prohibited.
+mulF2m :: BinaryPolynomial -- ^ Modulus
+       -> Integer
+       -> Integer
+       -> Integer
+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"
+    | otherwise = modF2m fx $ go (if n2 `mod` 2 == 1 then n1 else 0) (log2 n2)
+      where
+        go n s | s == 0  = n
+               | otherwise = if testBit n2 s
+                                then go (n `addF2m` shift n1 s) (s - 1)
+                                else go n (s - 1)
 {-# INLINABLE mulF2m #-}
 
 -- | Squaring over F₂m.
--- TODO: This is still slower than @mulF2m@.
-
--- Multiplication table? C?
-squareF2m :: BinaryPolynomial  -- ^ Irreducible binary polynomial
-          -> Integer -> Integer
-squareF2m fx = modF2m fx . square
+--
+-- This function is undefined for negative arguments, because their bit
+-- representation is platform-dependent. Zero modulus is also prohibited.
+squareF2m :: BinaryPolynomial -- ^ Modulus
+          -> Integer
+          -> Integer
+squareF2m fx = modF2m fx . squareF2m'
 {-# INLINE squareF2m #-}
 
-square :: Integer -> Integer
-square n1 = go n1 ln1
-  where
-    ln1 = log2 n1
-    go n s | s == 0 = n
-           | otherwise = go (x .|. y) (s - 1)
-      where
-        x = shift (shift n (2 * (s - ln1) - 1)) (2 * (ln1 - s) + 2)
-        y = n .&. (shift 1 (2 * (ln1 - s) + 1) - 1)
-{-# INLINE square #-}
-
--- | Inversion of @n over F₂m using extended Euclidean algorithm.
+-- | Squaring over F₂m without reduction by modulo.
 --
--- If @n doesn't have an inverse, Nothing is returned.
-invF2m :: BinaryPolynomial -- ^ Irreducible binary polynomial
-       -> Integer -> Maybe Integer
-invF2m _  0 = Nothing
-invF2m fx n
-    | n >= fx   = Nothing
-    | otherwise = go n fx 1 0
-    where
-      go u v g1 g2
-          | u == 1    = Just $ modF2m fx g1
-          | j < 0     = go u  (v  `xor` shift  u (-j)) g1 (g2 `xor` shift g1 (-j))
-          | otherwise = go (u  `xor` shift v  j) v (g1 `xor` shift g2 j) g2
-        where
-          j = log2 u - log2 v
+-- The implementation utilizes the fact that for binary polynomial S(x) we have
+-- S(x)^2 = S(x^2). In other words, insert a zero bit between every bits of argument: 1101 -> 1010001.
+--
+-- This function is undefined for negative arguments, because their bit
+-- representation is platform-dependent.
+squareF2m' :: Integer
+           -> Integer
+squareF2m' n
+    | n < 0     = error "mulF2m: negative number represent no binary binary polynomial"
+    | otherwise = foldl' (\acc s -> if testBit n s then setBit acc (2 * s) else acc) 0 [0 .. log2 n]
+{-# INLINE squareF2m' #-}
+
+-- | 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
+gcdF2m :: Integer
+       -> Integer
+       -> (Integer, Integer, Integer)
+gcdF2m a b = go (a, b, 1, 0, 0, 1)
+  where
+    go (g, 0, u, _, v, _)
+        = (g, u, v)
+    go (r0, r1, s0, s1, t0, t1)
+        = go (r1, r0 `addF2m` shift r1 j, s1, s0 `addF2m` shift s1 j, t1, t0 `addF2m` shift t1 j)
+            where j = max 0 (log2 r0 - log2 r1)
+
+-- | Modular inversion over F₂m.
+-- If @n@ doesn't have an inverse, 'Nothing' is returned.
+--
+-- This function is undefined for negative arguments, because their bit
+-- representation is platform-dependent. Zero modulus is also prohibited.
+invF2m :: BinaryPolynomial -- ^ Modulus
+       -> Integer
+       -> Maybe Integer
+invF2m fx n = if g == 1 then Just (modF2m fx u) else Nothing
+  where
+    (g, u, _) = gcdF2m n fx
 {-# INLINABLE invF2m #-}
 
 -- | Division over F₂m. If the dividend doesn't have an inverse it returns
 -- 'Nothing'.
 --
--- Compute n1 / n2
-divF2m :: BinaryPolynomial  -- ^ Irreducible binary polynomial
-       -> Integer  -- ^ Dividend
-       -> Integer  -- ^ Quotient
-       -> Maybe Integer
+-- This function is undefined for negative arguments, because their bit
+-- representation is platform-dependent. Zero modulus is also prohibited.
+divF2m :: BinaryPolynomial -- ^ Modulus
+       -> Integer          -- ^ Dividend
+       -> Integer          -- ^ Divisor
+       -> Maybe Integer    -- ^ Quotient
 divF2m fx n1 n2 = mulF2m fx n1 <$> invF2m fx n2
 {-# INLINE divF2m #-}
diff --git a/Crypto/PubKey/ECC/Prim.hs b/Crypto/PubKey/ECC/Prim.hs
index 139d28a..a3a8324 100644
--- a/Crypto/PubKey/ECC/Prim.hs
+++ b/Crypto/PubKey/ECC/Prim.hs
@@ -26,6 +26,13 @@ scalarGenerate curve = generateBetween 1 (n - 1)
 
 --TODO: Extract helper function for `fromMaybe PointO...`
 
+-- | Elliptic Curve point negation:
+-- @pointNegate c p@ returns point @q@ such that @pointAdd c p q == PointO@.
+pointNegate :: Curve -> Point -> Point
+pointNegate _           PointO     = PointO
+pointNegate CurveFP{}  (Point x y) = Point x (-y)
+pointNegate CurveF2m{} (Point x y) = Point x (x `addF2m` y)
+
 -- | Elliptic Curve point addition.
 --
 -- /WARNING:/ Vulnerable to timing attacks.
@@ -33,22 +40,21 @@ pointAdd :: Curve -> Point -> Point -> Point
 pointAdd _ PointO PointO = PointO
 pointAdd _ PointO q = q
 pointAdd _ p PointO = p
-pointAdd c@(CurveFP (CurvePrime pr _)) p@(Point xp yp) q@(Point xq yq)
-    | p == Point xq (-yq) = PointO
-    | p == q = pointDouble c p
-    | otherwise = fromMaybe PointO $ do
-                      s <- divmod (yp - yq) (xp - xq) pr
-                      let xr = (s ^ (2::Int) - xp - xq) `mod` pr
-                          yr = (s * (xp - xr) - yp) `mod` pr
-                      return $ Point xr yr
-pointAdd c@(CurveF2m (CurveBinary fx cc)) p@(Point xp yp) q@(Point xq yq)
-    | p == Point xq (xq `addF2m` yq) = PointO
-    | p == q = pointDouble c p
-    | otherwise = fromMaybe PointO $ do
-                     s <- divF2m fx (yp `addF2m` yq) (xp `addF2m` xq)
-                     let xr = mulF2m fx s s `addF2m` s `addF2m` xp `addF2m` xq `addF2m` a
-                         yr = mulF2m fx s (xp `addF2m` xr) `addF2m` xr `addF2m` yp
-                     return $ Point xr yr
+pointAdd c p q
+  | p == q = pointDouble c p
+  | p == pointNegate c q = PointO
+pointAdd (CurveFP (CurvePrime pr _)) (Point xp yp) (Point xq yq)
+    = fromMaybe PointO $ do
+        s <- divmod (yp - yq) (xp - xq) pr
+        let xr = (s ^ (2::Int) - xp - xq) `mod` pr
+            yr = (s * (xp - xr) - yp) `mod` pr
+        return $ Point xr yr
+pointAdd (CurveF2m (CurveBinary fx cc)) (Point xp yp) (Point xq yq)
+    = fromMaybe PointO $ do
+        s <- divF2m fx (yp `addF2m` yq) (xp `addF2m` xq)
+        let xr = mulF2m fx s s `addF2m` s `addF2m` xp `addF2m` xq `addF2m` a
+            yr = mulF2m fx s (xp `addF2m` xr) `addF2m` xr `addF2m` yp
+        return $ Point xr yr
   where a = ecc_a cc
 
 -- | Elliptic Curve point doubling.
@@ -95,8 +101,8 @@ pointBaseMul c n = pointMul c n (ecc_g $ common_curve c)
 -- /WARNING:/ Vulnerable to timing attacks.
 pointMul :: Curve -> Integer -> Point -> Point
 pointMul _ _ PointO = PointO
-pointMul c n p@(Point xp yp)
-    | n <  0 = pointMul c (-n) (Point xp (-yp))
+pointMul c n p
+    | n <  0 = pointMul c (-n) (pointNegate c p)
     | n == 0 = PointO
     | n == 1 = p
     | odd n = pointAdd c p (pointMul c (n - 1) p)
diff --git a/benchs/Number/F2m.hs b/benchs/Number/F2m.hs
new file mode 100644
index 0000000..408a192
--- /dev/null
+++ b/benchs/Number/F2m.hs
@@ -0,0 +1,53 @@
+{-# LANGUAGE PackageImports #-}
+
+module Main where
+
+import Criterion.Main
+import System.Random
+
+import "cryptonite" Crypto.Number.Basic (log2)
+import "cryptonite" Crypto.Number.F2m
+
+genInteger :: Int -> Int -> Integer
+genInteger salt bits
+    = head
+    . dropWhile ((< bits) . log2)
+    . scanl (\a r -> a * 2^31 + abs r) 0
+    . randoms
+    . mkStdGen
+    $ salt + bits
+
+benchMod :: Int -> Benchmark
+benchMod bits = bench (show bits) $ nf (modF2m m) a
+  where
+    m = genInteger 0 bits
+    a = genInteger 1 (2 * bits)
+
+benchMul :: Int -> Benchmark
+benchMul bits = bench (show bits) $ nf (mulF2m m a) b
+  where
+    m = genInteger 0 bits
+    a = genInteger 1 bits
+    b = genInteger 2 bits
+
+benchSquare :: Int -> Benchmark
+benchSquare bits = bench (show bits) $ nf (squareF2m m) a
+  where
+    m = genInteger 0 bits
+    a = genInteger 1 bits
+
+benchInv :: Int -> Benchmark
+benchInv bits = bench (show bits) $ nf (invF2m m) a
+  where
+    m = genInteger 0 bits
+    a = genInteger 1 bits
+
+bitsList :: [Int]
+bitsList = [64, 128, 256, 512, 1024, 2048]
+
+main = defaultMain
+    [ bgroup    "modF2m" $ map benchMod    bitsList
+    , bgroup    "mulF2m" $ map benchMul    bitsList
+    , bgroup "squareF2m" $ map benchSquare bitsList
+    , bgroup    "invF2m" $ map benchInv    bitsList
+    ]
diff --git a/cryptonite.cabal b/cryptonite.cabal
index 2450f0a..c2dc811 100644
--- a/cryptonite.cabal
+++ b/cryptonite.cabal
@@ -308,8 +308,10 @@ Test-Suite test-cryptonite
                      KAT_Camellia
                      KAT_Curve25519
                      KAT_DES
+                     KAT_Ed448
                      KAT_Ed25519
                      KAT_CMAC
+                     KAT_HKDF
                      KAT_HMAC
                      KAT_MiyaguchiPreneel
                      KAT_PBKDF2
@@ -323,6 +325,10 @@ Test-Suite test-cryptonite
                      KAT_RC4
                      KAT_Scrypt
                      KAT_TripleDES
+                     ChaChaPoly1305
+                     Number
+                     Number.F2m
+                     Padding
                      Poly1305
                      Salsa
                      Utils
diff --git a/tests/Number/F2m.hs b/tests/Number/F2m.hs
new file mode 100644
index 0000000..afa6e50
--- /dev/null
+++ b/tests/Number/F2m.hs
@@ -0,0 +1,83 @@
+module Number.F2m (tests) where
+
+import Imports hiding ((.&.))
+import Data.Bits
+import Crypto.Number.Basic (log2)
+import Crypto.Number.F2m
+
+addTests = testGroup "addF2m"
+    [ testProperty "commutative"
+        $ \a b -> a `addF2m` b == b `addF2m` a
+    , testProperty "associative"
+        $ \a b c -> (a `addF2m` b) `addF2m` c == a `addF2m` (b `addF2m` c)
+    , testProperty "0 is neutral"
+        $ \a -> a `addF2m` 0 == a
+    , testProperty "nullable"
+        $ \a -> a `addF2m` a == 0
+    , testProperty "works per bit"
+        $ \a b -> (a `addF2m` b) .&. b == (a .&. b) `addF2m` b
+    ]
+
+modTests = testGroup "modF2m"
+    [ testProperty "idempotent"
+        $ \(Positive m) (NonNegative a) -> modF2m m a == modF2m m (modF2m m a)
+    , testProperty "upper bound"
+        $ \(Positive m) (NonNegative a) -> modF2m m a < 2 ^ log2 m
+    , testProperty "reach upper"
+        $ \(Positive m) -> let a = 2 ^ log2 m - 1 in modF2m m (m `addF2m` a) == a
+    , testProperty "lower bound"
+        $ \(Positive m) (NonNegative a) -> modF2m m a >= 0
+    , testProperty "reach lower"
+        $ \(Positive m) -> modF2m m m == 0
+    , testProperty "additive"
+        $ \(Positive m) (NonNegative a) (NonNegative b)
+            -> modF2m m a `addF2m` modF2m m b == modF2m m (a `addF2m` b)
+    ]
+
+mulTests = testGroup "mulF2m"
+    [ testProperty "commutative"
+        $ \(Positive m) (NonNegative a) (NonNegative b) -> mulF2m m a b == mulF2m m b a
+    , testProperty "associative"
+        $ \(Positive m) (NonNegative a) (NonNegative b) (NonNegative c)
+            -> mulF2m m (mulF2m m a b) c == mulF2m m a (mulF2m m b c)
+    , testProperty "1 is neutral"
+        $ \(Positive m) (NonNegative a) -> mulF2m m a 1 == modF2m m a
+    , testProperty "0 is annihilator"
+        $ \(Positive m) (NonNegative a) -> mulF2m m a 0 == 0
+    , testProperty "distributive"
+        $ \(Positive m) (NonNegative a) (NonNegative b) (NonNegative c)
+            -> mulF2m m a (b `addF2m` c) == mulF2m m a b `addF2m` mulF2m m a c
+    ]
+
+squareTests = testGroup "squareF2m"
+    [ testProperty "sqr(a) == a * a"
+        $ \(Positive m) (NonNegative a) -> mulF2m m a a == squareF2m m a
+    ]
+
+invTests = testGroup "invF2m"
+    [ testProperty "1 / a * a == 1"
+        $ \(Positive m) (NonNegative a)
+            -> maybe True (\c -> mulF2m m c a == modF2m m 1) (invF2m m a)
+    , testProperty "1 / a == a (mod a^2-1)"
+        $ \(NonNegative a) -> a < 2 || invF2m (squareF2m' a `addF2m` 1) a == Just a
+    ]
+
+divTests = testGroup "divF2m"
+    [ testProperty "1 / a == inv a"
+        $ \(Positive m) (NonNegative a) -> divF2m m 1 a == invF2m m a
+    , testProperty "a / b == a * inv b"
+        $ \(Positive m) (NonNegative a) (NonNegative b)
+            -> 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)
+    ]
+
+tests = testGroup "number.F2m"
+    [ addTests
+    , modTests
+    , mulTests
+    , squareTests
+    , invTests
+    , divTests
+    ]
diff --git a/tests/Tests.hs b/tests/Tests.hs
index 9d2c017..1ab9ea5 100644
--- a/tests/Tests.hs
+++ b/tests/Tests.hs
@@ -4,6 +4,7 @@ module Main where
 import Imports
 
 import qualified Number
+import qualified Number.F2m
 import qualified BCrypt
 import qualified Hash
 import qualified Poly1305
@@ -33,6 +34,7 @@ import qualified Padding
 
 tests = testGroup "cryptonite"
     [ Number.tests
+    , Number.F2m.tests
     , Hash.tests
     , Padding.tests
     , testGroup "ConstructHash"