Point group element generator and subgroup testing

Former-commit-id: 8d3e89facdac8e0744ada736d5b7453f23dc9ddc

qpms/pointgroups.c

@@ -1,5 +1,6 @@
 #include "pointgroups.h"
 #include <search.h>
+#include <stdlib.h>
 
 #define PAIRCMP(a, b) {\
   if ((a) < (b)) return -1;\
@@ -89,6 +90,37 @@ qpms_irot3_t *qpms_pg_canonical_elems(qpms_irot3_t *target,
   return target;
 }
 
+qpms_irot3_t *qpms_pg_elems(qpms_irot3_t *target, qpms_pointgroup_t g) {
+  QPMS_UNTESTED;
+  target = qpms_pg_canonical_elems(target, g.c, g.n);
+  qpms_gmi_t order = qpms_pg_order(g.c, g.n);
+  const qpms_irot3_t o = g.orientation, o_inv = qpms_irot3_inv(o);
+  for(qpms_gmi_t i = 0 ; i < order; ++i) 
+    target[i] = qpms_irot3_mult(o_inv, qpms_irot3_mult(target[i], o));
+  return target;
+}
+
+_Bool qpms_pg_is_subgroup(qpms_pointgroup_t small, qpms_pointgroup_t big) {
+  QPMS_UNTESTED;
+  qpms_gmi_t order_big = qpms_pg_order(big.c, big.n);
+  qpms_gmi_t order_small = qpms_pg_order(small.c, small.n);
+  if (!order_big || !order_small)
+    QPMS_NOT_IMPLEMENTED("Subgroup testing for infinite groups not implemented");
+  if (order_big < order_small) return false;
+  // TODO maybe some other fast-track negative decisions
+
+  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);
+
+  for(qpms_gmi_t smalli = 0; smalli < order_small; ++smalli) {
+    qpms_irot3_t *match = bsearch(&elems_small[smalli], elems_big, order_big,
+        sizeof(qpms_irot3_t), qpms_pg_irot3_approx_cmp_v);
+    if (!match) return false;
+  }
+
+  return true;
+}
 
 /// Returns the order of a given 3D point group type.
 /** For infinite groups returns 0. */

qpms/pointgroups.h

@@ -74,4 +74,16 @@ qpms_gmi_t qpms_pg_genset(qpms_pointgroup_class cls, ///< Point group class.
 		qpms_irot3_t gen[] ///< Target generator array
 		);
 
+/// Generates all elements of a given point group.
+/** The order of elements corresponds to the one obtained from qpms_pg_canonical_elems().
+ */
+qpms_irot3_t *qpms_pg_elems(qpms_irot3_t *target, ///< Target array (optional; if NULL, a new one is allocated)
+	       	qpms_pointgroup_t g ///< Specification of the point group.
+		);
+
+/// Checks whether \a a is subgroup of \a b (in a dirty general way).
+_Bool qpms_pg_is_subgroup(qpms_pointgroup_t a, qpms_pointgroup_t b);
+
+
+
 #endif //POINTGROUPS_H

qpms/qpms_types.h

@@ -366,7 +366,7 @@ typedef enum {
 /// Full characterisation of a 3D point group.
 typedef struct qpms_pointgroup_t {
 	qpms_pointgroup_class c; ///< Point group classification.
-	size_t n; ///< Order of the rotational subgroup \f$ \mathrm{C_n} \f$ of finite axial groups.
+	qpms_gmi_t n; ///< Order of the rotational subgroup \f$ \mathrm{C_n} \f$ of finite axial groups.
 	/// Transformation between this point group and the "canonical" point group of the same type.
 	/**
 	 * Each 3D point group of a given type (specified by the \a c 

qpms/quaternions.h

@@ -203,6 +203,15 @@ static inline qpms_irot3_t qpms_irot3_mult(const qpms_irot3_t p, const qpms_irot
 	return r;
 }
 
+/// Improper rotation inverse operation.
+static inline qpms_irot3_t qpms_irot3_inv(qpms_irot3_t p) {
+#ifndef NDEBUG
+	qpms_irot3_checkdet(p);
+#endif
+	p.rot = qpms_quat_inv(p.rot);
+	return p;
+}
+
 /// Improper rotation power \f$ p^n \f$.
 static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) {
 #ifndef NDEBUG