diff --git a/qpms/pointgroups.c b/qpms/pointgroups.c
index 4daed5e..0125655 100644
--- a/qpms/pointgroups.c
+++ b/qpms/pointgroups.c
@@ -1,36 +1,33 @@
 #include "pointgroups.h"
 #include <search.h>
 #include <stdlib.h>
+#include <assert.h>
+
+double qpms_pg_quat_cmp_atol = QPMS_QUAT_ATOL;
 
 #define PAIRCMP(a, b) {\
   if ((a) < (b)) return -1;\
   if ((a) > (b)) return 1;\
 }
 
-// THIS IS NUMERICALLY UNSTABLE! FIXME!!!!!!
-int qpms_pg_irot3_cmp(const qpms_irot3_t *p1, const qpms_irot3_t *p2) {
-  QPMS_WARN("NUMERICALLY UNSTABLE! DON'T USE ME!")
+#define PAIRCMP_ATOL(a, b, atol) {\
+  if ((a) < (b) + (atol)) return -1;\
+  if ((a) + (atol) > (b)) return 1;\
+}
+
+int qpms_pg_irot3_approx_cmp(const qpms_irot3_t *p1, 
+    const qpms_irot3_t *p2, double atol) {
+  assert(atol >= 0);
   PAIRCMP(p1->det, p2->det);
-  const qpms_quat_t r1 = qpms_quat_standardise(p1->rot), r2 = qpms_quat_standardise(p2->rot);
-  PAIRCMP(creal(r1.a), creal(r2.a));
-  PAIRCMP(cimag(r1.a), cimag(r2.a));
-  PAIRCMP(creal(r1.b), creal(r2.b));
-  PAIRCMP(cimag(r1.b), cimag(r2.b));
+  const qpms_quat_t r1 = qpms_quat_standardise(p1->rot, atol),
+        r2 = qpms_quat_standardise(p2->rot, atol);
+  PAIRCMP_ATOL(creal(r1.a), creal(r2.a), atol);
+  PAIRCMP_ATOL(cimag(r1.a), cimag(r2.a), atol);
+  PAIRCMP_ATOL(creal(r1.b), creal(r2.b), atol);
+  PAIRCMP_ATOL(cimag(r1.b), cimag(r2.b), atol);
   return 0;
 }
 
-int qpms_pg_irot3_cmp_v(const void *av, const void *bv) {
-  const qpms_irot3_t *a = av, *b = bv;
-  return qpms_pg_irot3_cmp(a, b);
-}
-
-int qpms_pg_irot3_approx_cmp(const qpms_irot3_t *a, const qpms_irot3_t *b, double atol) {
-  if (qpms_irot3_isclose(*a, *b, atol)) return 0;
-  else return qpms_pg_irot3_cmp(a, b);
-}
-
-double qpms_pg_quat_cmp_atol = QPMS_QUAT_ATOL;
-
 int qpms_pg_irot3_approx_cmp_v(const void *av, const void *bv) {
   const qpms_irot3_t *a = av, *b = bv;
   return qpms_pg_irot3_approx_cmp(a, b, qpms_pg_quat_cmp_atol);
@@ -113,7 +110,7 @@ _Bool qpms_pg_is_subgroup(qpms_pointgroup_t small, qpms_pointgroup_t big) {
 
   qpms_irot3_t *elems_small = qpms_pg_elems(NULL, small);
   qpms_irot3_t *elems_big = qpms_pg_elems(NULL, small);
-  qsort(elems_big, order_big, sizeof(qpms_irot3_t), qpms_pg_irot3_cmp_v);
+  qsort(elems_big, order_big, sizeof(qpms_irot3_t), qpms_pg_irot3_approx_cmp_v);
 
   for(qpms_gmi_t smalli = 0; smalli < order_small; ++smalli) {
     qpms_irot3_t *match = bsearch(&elems_small[smalli], elems_big, order_big,
diff --git a/qpms/pointgroups.h b/qpms/pointgroups.h
index bc83a81..0fc92fe 100644
--- a/qpms/pointgroups.h
+++ b/qpms/pointgroups.h
@@ -1,7 +1,6 @@
 /*! \file pointgroups.h
  * \brief Quaternion-represented 3D point groups.
  */
-
 #ifndef POINTGROUPS_H
 #define POINTGROUPS_H
 
@@ -34,10 +33,6 @@ static inline _Bool qpms_pg_is_finite_axial(qpms_pointgroup_class cls) {
  */
 extern double qpms_pg_quat_cmp_atol;
 
-/// An ordering of qpms_irot3_t.
-int qpms_pg_irot3_cmp(const qpms_irot3_t *, const qpms_irot3_t *);
-/// A `search.h` and `qsort()` compatible ordering of qpms_irot3_t.
-int qpms_pg_irot3_cmp_v(const void *, const void *);
 /// An ordering of qpms_irot3_t that considers close enough elements equal.
 int qpms_pg_irot3_approx_cmp(const qpms_irot3_t *, const qpms_irot3_t *, 
 		double atol ///< Absolute tolerance for the quaternion part difference.
diff --git a/qpms/quaternions.h b/qpms/quaternions.h
index a1e93dd..eed1ce8 100644
--- a/qpms/quaternions.h
+++ b/qpms/quaternions.h
@@ -99,16 +99,15 @@ static inline bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, d
 
 /// "Standardises" a quaternion to have the largest component "positive".
 /**
- * FIXME
- * NUMERICALLY UNSTABLE. DON'T USE
  * This is to remove the ambiguity stemming from the double cover of SO(3).
  */
-static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p) {
+static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p, double atol) {
+	//assert(atol >= 0);
 	double maxabs = 0;
 	int maxi = 0;
 	const double *arr = (double *) &(p.a);
 	for(int i = 0; i < 4; ++i)
-		if (fabs(arr[i]) > maxabs) {
+		if (fabs(arr[i]) > maxabs + atol) {
 			maxi = i;
 			maxabs = fabs(arr[i]);
 		}
@@ -122,8 +121,8 @@ static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p) {
 /// Test approximate equality of "standardised" quaternions, so that \f$-q\f$ is considered equal to \f$q\f$.
 static inline bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
 	return qpms_quat_norm(qpms_quat_sub(
-				qpms_quat_standardise(p),
-				qpms_quat_standardise(q))) <= atol;
+				qpms_quat_standardise(p, atol),
+				qpms_quat_standardise(q, atol))) <= atol;
 }
 
 /// Norm of the quaternion imaginary (vector) part.