Untested generation of finite point groups

Former-commit-id: d61677f79685f69ba59f2f81ec7737b70b42d218

qpms/pointgroups.c

@@ -0,0 +1,95 @@
+#include "pointgroups.h"
+#include <search.h>
+
+#define PAIRCMP(a, b) {\
+  if ((a) < (b)) return -1;\
+  if ((a) > (b)) return 1;\
+}
+
+int qpms_pg_irot3_cmp(const qpms_irot3_t *a, const qpms_irot3_t *b) {
+  PAIRCMP(a->det, b->det);
+  PAIRCMP(creal(a->rot.a), creal(b->rot.a));
+  PAIRCMP(cimag(a->rot.a), cimag(b->rot.a));
+  PAIRCMP(creal(a->rot.b), creal(b->rot.b));
+  PAIRCMP(cimag(a->rot.b), cimag(b->rot.b));
+  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);
+}
+
+
+/// Generates the canonical elements of a given 3D point group type.
+qpms_irot3_t *qpms_pg_canonical_elems(qpms_irot3_t *target,
+    qpms_pointgroup_class cls, const qpms_gmi_t then) {
+  QPMS_UNTESTED;
+  qpms_gmi_t order = qpms_pg_order(cls, then);
+  QPMS_ENSURE(order, "Cannot generate an infinite group!");
+  if (!target) QPMS_CRASHING_MALLOC(target, order * sizeof(qpms_irot3_t));
+  target[0] = QPMS_IROT3_IDENTITY;
+  qpms_gmi_t ngen = qpms_pg_genset_size(cls, then);
+  qpms_irot3_t gens[ngen]; 
+  (void) qpms_pg_genset(cls, then, gens);
+
+  // Let's try it with a binary search tree, as an exercise :)
+  qpms_irot3_t tree[order];
+  void *root = NULL;
+  // put the starting element (identity) to the tree
+  (void) tsearch((void *) target, &root, qpms_pg_irot3_approx_cmp_v);
+  qpms_gmi_t n = 1; // No. of generated elements.
+
+  // And let's do the DFS without recursion; the "stack size" here (=order) might be excessive, but whatever
+  qpms_gmi_t gistack[order], //< generator indices (related to gens[])
+             srcstack[order], //< pre-image indices (related to target[])
+             si = 0; //< stack index
+  gistack[0] = 0;
+  srcstack[0] = 0;
+
+  while (si >= 0) { // DFS
+    if (gistack[si] < ngen) { // are there generators left at this level? If so, try another elem
+      if (n >= order) QPMS_WTF; // TODO some error message
+      target[n] = qpms_irot3_mult(gens[gistack[si]], target[srcstack[si]]);
+      if (tfind((void *) &(target[n]), &root, qpms_pg_irot3_approx_cmp_v))
+        // elem found, try it with another gen in the next iteration
+        gistack[si]++;
+      else {
+        // elem not found, add it to the tree and proceed to next level
+        (void) tsearch( &(target[n]), &root, qpms_pg_irot3_approx_cmp_v);
+        ++si;
+        gistack[si] = 0;
+        srcstack[si] = n;
+        ++n;
+      }
+    } else { // no generators left at this level, get to the previous level
+      --si; 
+      if (si >= 0) ++gistack[si];
+    }
+  }
+
+  QPMS_ENSURE(n == order, "Point group generation failure "
+      "(assumed group order = %d, got %d; qpms_pg_quat_cmp_atol = %g)",
+      order, n, qpms_pg_quat_cmp_atol);
+
+  while(root) tdelete(root, &root, qpms_pg_irot3_approx_cmp_v); // I hope this is the correct way.
+
+  return target;
+}
+
+
+
+
+

qpms/pointgroups.h

@@ -1,3 +1,7 @@
+/*! \file pointgroups.h
+ * \brief Quaternion-represented 3D point groups.
+ */
+
 #ifndef POINTGROUPS_H
 #define POINTGROUPS_H
 
@@ -5,6 +9,7 @@
 #include "quaternions.h"
 
 
+/// Returns true if the point group class belongs to one of the seven "axial" group series.
 static inline _Bool qpms_pg_is_finite_axial(qpms_pointgroup_class cls) {
 	switch(cls) {
         case QPMS_PGS_CN: 
@@ -20,10 +25,32 @@ static inline _Bool qpms_pg_is_finite_axial(qpms_pointgroup_class cls) {
 	}
 }
 
+/// Absolute tolerance threshold used internally to consider two different `qpms_irot3_t` instances equal.
+/**
+ * Used by @ref qpms_pg_irot3_approx_cmp_v.
+ * By default, set to @ref QPMS_QUAT_ATOL.
+ * It should work fine if the point group orders stay reasonable (thousands or less).
+ * Do not touch if unsure what you are doing.
+ */
+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.
+		);
+/// A `search.h` compatible ordering of qpms_irot3_t that considers close enough elements equal.
+/** The tolerance is determined by global variable @ref qpms_pg_quat_cmp_atol.
+ */
+int qpms_pg_irot3_approx_cmp_v(const void *, const void *);
+
 /// Returns the order of a given 3D point group type.
 /** For infinite groups returns 0. */
-static inline size_t qpms_pg_order(qpms_pointgroup_class cls, ///< Point group class.
-		size_t n ///< Number of rotations around main axis (only for finite axial groups).
+static inline qpms_gmi_t qpms_pg_order(qpms_pointgroup_class cls, ///< Point group class.
+		qpms_gmi_t n ///< Number of rotations around main axis (only for finite axial groups).
 		) {
 	if (qpms_pg_is_finite_axial(cls))
 		QPMS_ENSURE(n > 0, "n must be at least 1 for finite axial groups");
@@ -70,10 +97,18 @@ static inline size_t qpms_pg_order(qpms_pointgroup_class cls, ///< Point group c
 	}
 }
 
+/// Generates the canonical elements of a given 3D point group type.
+/** Uses the canonical generators and DPS. */
+qpms_irot3_t *qpms_pg_canonical_elems(
+		qpms_irot3_t *target, ///< Target array (optional; if NULL, a new one is allocated)
+		qpms_pointgroup_class cls, ///< Point group class.
+		qpms_gmi_t ///< Number of rotations around \a z axis (only for axial group classes).
+		);
+
 /// Returns the number of canonical generators of a given 3D point group type.
 /** TODO what does it do for infinite (Lie) groups? */
-static inline size_t qpms_pg_genset_size(qpms_pointgroup_class cls, ///< Point group class.
-		size_t n ///< Number of rotations around main axis (only for axial groups).
+static inline qpms_gmi_t qpms_pg_genset_size(qpms_pointgroup_class cls, ///< Point group class.
+		qpms_gmi_t n ///< Number of rotations around main axis (only for axial groups).
 		) {
 	if (qpms_pg_is_finite_axial(cls)) {
 		QPMS_ENSURE(n > 0, "n must be at least 1 for finite axial groups");
@@ -128,8 +163,8 @@ static inline size_t qpms_pg_genset_size(qpms_pointgroup_class cls, ///< Point g
 
 /// Fills an array of canonical generators of a given 3D point group type.
 /** TODO what does it do for infinite (Lie) groups? */
-static inline size_t qpms_pg_genset(qpms_pointgroup_class cls, ///< Point group class.
-		size_t n, ///< Number of rotations around main axis (only for axial groups).
+static inline qpms_gmi_t qpms_pg_genset(qpms_pointgroup_class cls, ///< Point group class.
+		qpms_gmi_t n, ///< Number of rotations around main axis (only for axial groups).
 		qpms_irot3_t gen[] ///< Target generator array
 		) {
 	if (qpms_pg_is_finite_axial(cls)) {