diff --git a/Crypto/ECC/Edwards25519.hs b/Crypto/ECC/Edwards25519.hs
index c6dffd5..ebba8ca 100644
--- a/Crypto/ECC/Edwards25519.hs
+++ b/Crypto/ECC/Edwards25519.hs
@@ -37,6 +37,10 @@
 --
 --     - or cofactor must be cleared by multiplying points by 8
 --
+--     Utility functions are provided to implement this.  Testing
+--     subgroup membership with 'pointHasPrimeOrder' is 50-time slower
+--     than call 'pointMulByCofactor'.
+--
 -- 2. Scalar arithmetic is always reduced modulo L, allowing fixed
 -- length and constant execution time, but this reduction is valid
 -- only when points are in the prime-order subgroup.
@@ -56,6 +60,7 @@ module Crypto.ECC.Edwards25519
     -- * Points
     , pointDecode
     , pointEncode
+    , pointHasPrimeOrder
     -- * Arithmetic functions
     , toPoint
     , scalarAdd
@@ -64,6 +69,7 @@ module Crypto.ECC.Edwards25519
     , pointAdd
     , pointDouble
     , pointMul
+    , pointMulByCofactor
     , pointsMulVarTime
     ) where
 
@@ -203,6 +209,19 @@ pointDecode bs
                     else return $ CryptoPassed (Point p)
 {-# NOINLINE pointDecode #-}
 
+-- | Test whether a point belongs to the prime-order subgroup
+-- generated by the base point.  Result is 'True' for the identity
+-- point.
+--
+-- @
+-- pointHasPrimeOrder p = 'pointNegate' p == 'pointMul' l_minus_one p
+-- @
+pointHasPrimeOrder :: Point -> Bool
+pointHasPrimeOrder (Point p) = unsafeDoIO $
+    withByteArray p $ \pp ->
+        fmap (/= 0) (ed25519_point_has_prime_order pp)
+{-# NOINLINE pointHasPrimeOrder #-}
+
 -- | Negate a point.
 pointNegate :: Point -> Point
 pointNegate (Point a) =
@@ -229,6 +248,17 @@ pointDouble (Point a) =
         withByteArray a $ \pa ->
              ed25519_point_double out pa
 
+-- | Multiply a point by h = 8.
+--
+-- @
+-- pointMulByCofactor p = 'pointMul' scalar_8 p
+-- @
+pointMulByCofactor :: Point -> Point
+pointMulByCofactor (Point a) =
+    Point $ B.allocAndFreeze pointArraySize $ \out ->
+        withByteArray a $ \pa ->
+             ed25519_point_mul_by_cofactor out pa
+
 -- | Scalar multiplication over curve edwards25519.
 --
 -- Note: when the scalar had reduction modulo L and the input point
@@ -299,6 +329,10 @@ foreign import ccall "cryptonite_ed25519_point_eq"
                      -> Ptr Point
                      -> IO CInt
 
+foreign import ccall "cryptonite_ed25519_point_has_prime_order"
+    ed25519_point_has_prime_order :: Ptr Point
+                                  -> IO CInt
+
 foreign import ccall "cryptonite_ed25519_point_negate"
     ed25519_point_negate :: Ptr Point -- minus_a
                          -> Ptr Point -- a
@@ -315,6 +349,11 @@ foreign import ccall "cryptonite_ed25519_point_double"
                          -> Ptr Point -- a
                          -> IO ()
 
+foreign import ccall "cryptonite_ed25519_point_mul_by_cofactor"
+    ed25519_point_mul_by_cofactor :: Ptr Point -- eight_a
+                                  -> Ptr Point -- a
+                                  -> IO ()
+
 foreign import ccall "cryptonite_ed25519_point_base_scalarmul"
     ed25519_point_base_scalarmul :: Ptr Point  -- scaled
                                  -> Ptr Scalar -- scalar
diff --git a/cbits/ed25519/ed25519-cryptonite-exts.h b/cbits/ed25519/ed25519-cryptonite-exts.h
index 530c8cf..8a74618 100644
--- a/cbits/ed25519/ed25519-cryptonite-exts.h
+++ b/cbits/ed25519/ed25519-cryptonite-exts.h
@@ -87,6 +87,25 @@ ED25519_FN(ed25519_point_eq) (const ge25519 *p, const ge25519 *q) {
     return eq;
 }
 
+static int
+ED25519_FN(ed25519_point_is_identity) (const ge25519 *p) {
+    static const unsigned char zero[32] = {0};
+    unsigned char check[32];
+    bignum25519 d;
+    int eq;
+
+    // pX = 0
+    curve25519_contract(check, p->x);
+    eq = ed25519_verify(check, zero, 32);
+
+    // pY - pZ = 0
+    curve25519_sub_reduce(d, p->y, p->z);
+    curve25519_contract(check, d);
+    eq &= ed25519_verify(check, zero, 32);
+
+    return eq;
+}
+
 void
 ED25519_FN(ed25519_point_negate) (ge25519 *r, const ge25519 *p) {
     curve25519_neg(r->x, p->x);
@@ -105,6 +124,13 @@ ED25519_FN(ed25519_point_double) (ge25519 *r, const ge25519 *p) {
     ge25519_double(r, p);
 }
 
+void
+ED25519_FN(ed25519_point_mul_by_cofactor) (ge25519 *r, const ge25519 *p) {
+    ge25519_double_partial(r, p);
+    ge25519_double_partial(r, r);
+    ge25519_double(r, r);
+}
+
 void
 ED25519_FN(ed25519_point_base_scalarmul) (ge25519 *r, const bignum256modm s) {
     ge25519_scalarmult_base_niels(r, ge25519_niels_base_multiples, s);
@@ -207,3 +233,14 @@ ED25519_FN(ed25519_base_double_scalarmul_vartime) (ge25519 *r, const bignum256mo
     // computes [s1]basepoint + [s2]p2
     ge25519_double_scalarmult_vartime(r, p2, s2, s1);
 }
+
+int
+ED25519_FN(ed25519_point_has_prime_order) (const ge25519 *p) {
+    static const bignum256modm sc_zero = {0};
+    ge25519 q;
+
+    // computes Q = m.P, vartime allowed because m is not secret
+    ED25519_FN(ed25519_base_double_scalarmul_vartime) (&q, sc_zero, p, modm_m);
+
+    return ED25519_FN(ed25519_point_is_identity) (&q);
+}
diff --git a/tests/ECC/Edwards25519.hs b/tests/ECC/Edwards25519.hs
index 5ad4d0f..602ae72 100644
--- a/tests/ECC/Edwards25519.hs
+++ b/tests/ECC/Edwards25519.hs
@@ -93,6 +93,10 @@ tests = testGroup "ECC.Edwards25519"
             p0 `propertyEq` pointAdd p (pointNegate p)
         , testProperty "doubling" $ \p ->
             pointAdd p p `propertyEq` pointDouble p
+        , testProperty "multiplication by cofactor" $ \p ->
+            pointMul s8 p `propertyEq` pointMulByCofactor p
+        , testProperty "prime order" $ \(PrimeOrder p) ->
+            True `propertyEq` pointHasPrimeOrder p
         , testCase "8-torsion point" $ do
             assertBool "mul by 4" $ p0 /= pointMul s4 torsion8
             assertBool "mul by 8" $ p0 == pointMul s8 torsion8
@@ -102,8 +106,8 @@ tests = testGroup "ECC.Edwards25519"
             p `propertyEq` pointMul s1 p
         , testProperty "scalarmult with two" $ \p ->
             pointDouble p `propertyEq` pointMul s2 p
-        , testProperty "scalarmult with curve order - 1" $ \(PrimeOrder p) ->
-            pointNegate p `propertyEq` pointMul sI p
+        , testProperty "scalarmult with curve order - 1" $ \p ->
+            pointHasPrimeOrder p === (pointNegate p == pointMul sI p)
         , testProperty "scalarmult commutative" $ \a b ->
             pointMul a (toPoint b) === pointMul b (toPoint a)
         , testProperty "scalarmult distributive" $ \x y (PrimeOrder p) ->