diff --git a/qpms/pointgroups.c b/qpms/pointgroups.c new file mode 100644 index 0000000..bec4334 --- /dev/null +++ b/qpms/pointgroups.c @@ -0,0 +1,95 @@ +#include "pointgroups.h" +#include + +#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; +} + + + + + diff --git a/qpms/pointgroups.h b/qpms/pointgroups.h index d4c7490..46e4a85 100644 --- a/qpms/pointgroups.h +++ b/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)) {