QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

This might somehow more stable 

Former-commit-id: bc836c111d6a0bd14aa539c7d35e8987ff96d3da
This commit is contained in:
Marek Nečada 2019-07-25 07:59:43 +03:00
parent b939b1c177
commit 22d9f05ef8
3 changed files with 23 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
qpms/pointgroups.c
View File

@ -1,36 +1,33 @@
#include "pointgroups.h" #include "pointgroups.h"
#include <search.h> #include <search.h>
#include <stdlib.h> #include <stdlib.h>
#include <assert.h>
double qpms_pg_quat_cmp_atol = QPMS_QUAT_ATOL;
#define PAIRCMP(a, b) {\ #define PAIRCMP(a, b) {\
if ((a) < (b)) return -1;\ if ((a) < (b)) return -1;\
if ((a) > (b)) return 1;\ if ((a) > (b)) return 1;\
} }
// THIS IS NUMERICALLY UNSTABLE! FIXME!!!!!! #define PAIRCMP_ATOL(a, b, atol) {\
int qpms_pg_irot3_cmp(const qpms_irot3_t *p1, const qpms_irot3_t *p2) { if ((a) < (b) + (atol)) return -1;\
QPMS_WARN("NUMERICALLY UNSTABLE! DON'T USE ME!") 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); PAIRCMP(p1->det, p2->det);
const qpms_quat_t r1 = qpms_quat_standardise(p1->rot), r2 = qpms_quat_standardise(p2->rot); const qpms_quat_t r1 = qpms_quat_standardise(p1->rot, atol),
PAIRCMP(creal(r1.a), creal(r2.a)); r2 = qpms_quat_standardise(p2->rot, atol);
PAIRCMP(cimag(r1.a), cimag(r2.a)); PAIRCMP_ATOL(creal(r1.a), creal(r2.a), atol);
PAIRCMP(creal(r1.b), creal(r2.b)); PAIRCMP_ATOL(cimag(r1.a), cimag(r2.a), atol);
PAIRCMP(cimag(r1.b), cimag(r2.b)); PAIRCMP_ATOL(creal(r1.b), creal(r2.b), atol);
PAIRCMP_ATOL(cimag(r1.b), cimag(r2.b), atol);
return 0; 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) { int qpms_pg_irot3_approx_cmp_v(const void *av, const void *bv) {
const qpms_irot3_t *a = av, *b = bv; const qpms_irot3_t *a = av, *b = bv;
return qpms_pg_irot3_approx_cmp(a, b, qpms_pg_quat_cmp_atol); 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_small = qpms_pg_elems(NULL, small);
qpms_irot3_t *elems_big = 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) { for(qpms_gmi_t smalli = 0; smalli < order_small; ++smalli) {
qpms_irot3_t *match = bsearch(&elems_small[smalli], elems_big, order_big, qpms_irot3_t *match = bsearch(&elems_small[smalli], elems_big, order_big,

5
qpms/pointgroups.h
View File

@ -1,7 +1,6 @@
/*! \file pointgroups.h /*! \file pointgroups.h
* \brief Quaternion-represented 3D point groups. * \brief Quaternion-represented 3D point groups.
*/ */
#ifndef POINTGROUPS_H #ifndef POINTGROUPS_H
#define 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; 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. /// 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 *, int qpms_pg_irot3_approx_cmp(const qpms_irot3_t *, const qpms_irot3_t *,
double atol ///< Absolute tolerance for the quaternion part difference. double atol ///< Absolute tolerance for the quaternion part difference.

11
qpms/quaternions.h
View File

@ -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". /// "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). * 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; double maxabs = 0;
int maxi = 0; int maxi = 0;
const double *arr = (double *) &(p.a); const double *arr = (double *) &(p.a);
for(int i = 0; i < 4; ++i) for(int i = 0; i < 4; ++i)
if (fabs(arr[i]) > maxabs) { if (fabs(arr[i]) > maxabs + atol) {
maxi = i; maxi = i;
maxabs = fabs(arr[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$. /// 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) { 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( return qpms_quat_norm(qpms_quat_sub(
qpms_quat_standardise(p), qpms_quat_standardise(p, atol),
qpms_quat_standardise(q))) <= atol; qpms_quat_standardise(q, atol))) <= atol;
} }
/// Norm of the quaternion imaginary (vector) part. /// Norm of the quaternion imaginary (vector) part.