qpms/pointgroups.c

   1 #include "pointgroups.h"
   2 #include <search.h>
   3 #include <stdlib.h>
   4 #include <assert.h>
   5
   6 double qpms_pg_quat_cmp_atol = QPMS_QUAT_ATOL;
   7
   8 #define PAIRCMP(a, b) {\
   9   if ((a) < (b)) return -1;\
  10   if ((a) > (b)) return 1;\
  11 }
  12
  13 #define PAIRCMP_ATOL(a, b, atol) {\
  14   if ((a) < (b) + (atol)) return -1;\
  15   if ((a) + (atol) > (b)) return 1;\
  16 }
  17
  18 int qpms_pg_irot3_approx_cmp(const qpms_irot3_t *p1,
  19     const qpms_irot3_t *p2, double atol) {
  20   assert(atol >= 0);
  21   PAIRCMP(p1->det, p2->det);
  22   const qpms_quat_t r1 = qpms_quat_standardise(p1->rot, atol),
  23         r2 = qpms_quat_standardise(p2->rot, atol);
  24   PAIRCMP_ATOL(creal(r1.a), creal(r2.a), atol);
  25   PAIRCMP_ATOL(cimag(r1.a), cimag(r2.a), atol);
  26   PAIRCMP_ATOL(creal(r1.b), creal(r2.b), atol);
  27   PAIRCMP_ATOL(cimag(r1.b), cimag(r2.b), atol);
  28   return 0;
  29 }
  30
  31 int qpms_pg_irot3_approx_cmp_v(const void *av, const void *bv) {
  32   const qpms_irot3_t *a = av, *b = bv;
  33   return qpms_pg_irot3_approx_cmp(a, b, qpms_pg_quat_cmp_atol);
  34 }
  35
  36
  37 /// Generates the canonical elements of a given 3D point group type.
  38 qpms_irot3_t *qpms_pg_canonical_elems(qpms_irot3_t *target,
  39     qpms_pointgroup_class cls, const qpms_gmi_t then) {
  40   QPMS_UNTESTED;
  41   const qpms_gmi_t order = qpms_pg_order(cls, then);
  42   QPMS_ENSURE(order, "Cannot generate an infinite group!");
  43   if (!target) QPMS_CRASHING_MALLOC(target, (1+order) * sizeof(qpms_irot3_t));
  44   target[0] = QPMS_IROT3_IDENTITY;
  45   qpms_gmi_t ngen = qpms_pg_genset_size(cls, then);
  46   qpms_irot3_t gens[ngen];
  47   (void) qpms_pg_genset(cls, then, gens);
  48
  49   // Let's try it with a binary search tree, as an exercise :)
  50   void *root = NULL;
  51   // put the starting element (identity) to the tree
  52   (void) tsearch((void *) target, &root, qpms_pg_irot3_approx_cmp_v);
  53   qpms_gmi_t n = 1; // No. of generated elements.
  54
  55   // And let's do the DFS without recursion; the "stack size" here (=order) might be excessive, but whatever
  56   qpms_gmi_t gistack[order], //< generator indices (related to gens[])
  57              srcstack[order], //< pre-image indices (related to target[])
  58              si = 0; //< stack index
  59   gistack[0] = 0;
  60   srcstack[0] = 0;
  61
  62   while (si >= 0) { // DFS
  63     if (gistack[si] < ngen) { // are there generators left at this level? If so, try another elem
  64       if (n > order) QPMS_WTF; // TODO some error message
  65       target[n] = qpms_irot3_mult(gens[gistack[si]], target[srcstack[si]]);
  66       if (tfind((void *) &(target[n]), &root, qpms_pg_irot3_approx_cmp_v))
  67         // elem found, try it with another gen in the next iteration
  68         gistack[si]++;
  69       else {
  70         // elem not found, add it to the tree and proceed to next level
  71         (void) tsearch( &(target[n]), &root, qpms_pg_irot3_approx_cmp_v);
  72         ++si;
  73         gistack[si] = 0;
  74         srcstack[si] = n;
  75         ++n;
  76       }
  77     } else { // no generators left at this level, get to the previous level
  78       --si;
  79       if (si >= 0) ++gistack[si];
  80     }
  81   }
  82
  83   QPMS_ENSURE(n == order, "Point group generation failure "
  84       "(assumed group order = %d, got %d; qpms_pg_quat_cmp_atol = %g)",
  85       order, n, qpms_pg_quat_cmp_atol);
  86
  87   while(root) tdelete(*(qpms_irot3_t **)root, &root, qpms_pg_irot3_approx_cmp_v); // I hope this is the correct way.
  88
  89   return target;
  90 }
  91
  92 qpms_irot3_t *qpms_pg_elems(qpms_irot3_t *target, qpms_pointgroup_t g) {
  93   QPMS_UNTESTED;
  94   target = qpms_pg_canonical_elems(target, g.c, g.n);
  95   qpms_gmi_t order = qpms_pg_order(g.c, g.n);
  96   const qpms_irot3_t o = g.orientation, o_inv = qpms_irot3_inv(o);
  97   for(qpms_gmi_t i = 0 ; i < order; ++i)
  98     target[i] = qpms_irot3_mult(o_inv, qpms_irot3_mult(target[i], o));
  99   return target;
 100 }
 101
 102 _Bool qpms_pg_is_subgroup(qpms_pointgroup_t small, qpms_pointgroup_t big) {
 103   QPMS_UNTESTED;
 104   qpms_gmi_t order_big = qpms_pg_order(big.c, big.n);
 105   qpms_gmi_t order_small = qpms_pg_order(small.c, small.n);
 106   if (!order_big || !order_small)
 107     QPMS_NOT_IMPLEMENTED("Subgroup testing for infinite groups not implemented");
 108   if (order_big < order_small) return false;
 109   // TODO maybe some other fast-track negative decisions
 110
 111   qpms_irot3_t *elems_small = qpms_pg_elems(NULL, small);
 112   qpms_irot3_t *elems_big = qpms_pg_elems(NULL, small);
 113   qsort(elems_big, order_big, sizeof(qpms_irot3_t), qpms_pg_irot3_approx_cmp_v);
 114
 115   for(qpms_gmi_t smalli = 0; smalli < order_small; ++smalli) {
 116     qpms_irot3_t *match = bsearch(&elems_small[smalli], elems_big, order_big,
 117         sizeof(qpms_irot3_t), qpms_pg_irot3_approx_cmp_v);
 118     if (!match) return false;
 119   }
 120
 121   return true;
 122 }
 123
 124 /// Returns the order of a given 3D point group type.
 125 /** For infinite groups returns 0. */
 126 qpms_gmi_t qpms_pg_order(qpms_pointgroup_class cls, ///< Point group class.
 127                 qpms_gmi_t n ///< Number of rotations around main axis (only for finite axial groups).
 128                 ) {
 129         if (qpms_pg_is_finite_axial(cls))
 130                 QPMS_ENSURE(n > 0, "n must be at least 1 for finite axial groups");
 131         switch(cls) {
 132         // Axial groups
 133         case QPMS_PGS_CN:
 134                 return n;
 135         case QPMS_PGS_S2N:
 136         case QPMS_PGS_CNH:
 137         case QPMS_PGS_CNV:
 138         case QPMS_PGS_DN:
 139                 return 2*n;
 140         case QPMS_PGS_DND:
 141         case QPMS_PGS_DNH:
 142                 return 4*n;
 143
 144         // Remaining polyhedral groups
 145         case QPMS_PGS_T:
 146                 return 12;
 147         case QPMS_PGS_TD:
 148         case QPMS_PGS_TH:
 149         case QPMS_PGS_O:
 150                 return 24;
 151         case QPMS_PGS_OH:
 152                 return 48;
 153         case QPMS_PGS_I:
 154                 return 60;
 155         case QPMS_PGS_IH:
 156                 return 120;
 157
 158         // Continuous axial groups
 159         case QPMS_PGS_CINF:
 160         case QPMS_PGS_CINFH:
 161         case QPMS_PGS_CINFV:
 162         case QPMS_PGS_DINF:
 163         case QPMS_PGS_DINFH:
 164
 165         // Remaining continuous groups
 166         case QPMS_PGS_SO3:
 167         case QPMS_PGS_O3:
 168                 return 0; // 0 is infinity :-)
 169         default:
 170                 QPMS_WTF;
 171         }
 172 }
 173
 174
 175 /// Returns the number of canonical generators of a given 3D point group type.
 176 /** TODO what does it do for infinite (Lie) groups? */
 177 qpms_gmi_t qpms_pg_genset_size(qpms_pointgroup_class cls, ///< Point group class.
 178                 qpms_gmi_t n ///< Number of rotations around main axis (only for axial groups).
 179                 ) {
 180         if (qpms_pg_is_finite_axial(cls)) {
 181                 QPMS_ENSURE(n > 0, "n must be at least 1 for finite axial groups");
 182                 // n = 1 needs special care:
 183                 if (n==1)
 184                         switch(cls) {
 185                                 case QPMS_PGS_CN:  return 0; // triv.
 186                                 case QPMS_PGS_S2N: return 1; // Z_2
 187                                 case QPMS_PGS_CNH: return 1; // Dih_1
 188                                 case QPMS_PGS_CNV: return 1; // Dih_1
 189                                 case QPMS_PGS_DN:  return 1; // Dih_1
 190                                 case QPMS_PGS_DND: return 2; // Dih_2
 191                                 case QPMS_PGS_DNH: return 2; // Dih_1 x Dih_1
 192                                 default: QPMS_WTF;
 193                         }
 194         }
 195
 196         switch(cls) {
 197         // Axial groups
 198         case QPMS_PGS_CN:  return 1; // Z_n
 199         case QPMS_PGS_S2N: return 1; // Z_{2n}
 200         case QPMS_PGS_CNH: return 2; // Z_n x Dih_1
 201         case QPMS_PGS_CNV: return 2; // Dih_n
 202         case QPMS_PGS_DN:  return 2; // Dih_n
 203         case QPMS_PGS_DND: return 2; // Dih_2n
 204         case QPMS_PGS_DNH: return 3; // Dih_n x Dih_1
 205
 206         // Remaining polyhedral groups
 207         case QPMS_PGS_T:   // return ???; // A_4
 208         case QPMS_PGS_TD:  // return 2; // S_4
 209         case QPMS_PGS_TH:  // A_4 x Z_2
 210         case QPMS_PGS_O:   // S_4
 211         case QPMS_PGS_OH:  // return 3; //  S_4 x Z_2
 212         case QPMS_PGS_I:   // A_5
 213         case QPMS_PGS_IH:  // A_5 x Z_2
 214
 215         // Continuous axial groups
 216         case QPMS_PGS_CINF:
 217         case QPMS_PGS_CINFH:
 218         case QPMS_PGS_CINFV:
 219         case QPMS_PGS_DINF:
 220         case QPMS_PGS_DINFH:
 221
 222         // Remaining continuous groups
 223         case QPMS_PGS_SO3:
 224         case QPMS_PGS_O3:
 225                 QPMS_NOT_IMPLEMENTED("Not yet implemented for this point group class.");
 226         default:
 227                 QPMS_WTF;
 228         }
 229 }
 230
 231 /// Fills an array of canonical generators of a given 3D point group type.
 232 /** TODO what does it do for infinite (Lie) groups? */
 233 qpms_gmi_t qpms_pg_genset(qpms_pointgroup_class cls, ///< Point group class.
 234                 qpms_gmi_t n, ///< Number of rotations around main axis (only for axial groups).
 235                 qpms_irot3_t gen[] ///< Target generator array
 236                 ) {
 237         if (qpms_pg_is_finite_axial(cls)) {
 238                 QPMS_ENSURE(n > 0, "n must be at least 1 for finite axial groups");
 239                 // n = 1 needs special care:
 240                 if (n==1)
 241                         switch(cls) {
 242                                 case QPMS_PGS_CN:
 243                                         return 0; // triv.
 244                                 case QPMS_PGS_S2N:
 245                                         gen[0] = QPMS_IROT3_INVERSION;
 246                                         return 1; // Z_2
 247                                 case QPMS_PGS_CNH:
 248                                         gen[0] = QPMS_IROT3_ZFLIP;
 249                                         return 1; // Dih_1
 250                                 case QPMS_PGS_CNV:
 251                                         gen[0] = QPMS_IROT3_XFLIP;
 252                                         return 1; // Dih_1
 253                                 case QPMS_PGS_DN:
 254                                         gen[0] = QPMS_IROT3_XROT_PI; // CHECKME
 255                                         return 1; // Dih_1
 256                                 case QPMS_PGS_DND:
 257                                         gen[0] = QPMS_IROT3_INVERSION;
 258                                         gen[1] = QPMS_IROT3_XFLIP;
 259                                         return 2; // Dih_2
 260                                 case QPMS_PGS_DNH: // CHECKME
 261                                         gen[0] = QPMS_IROT3_ZFLIP;
 262                                         gen[1] = QPMS_IROT3_XFLIP;
 263                                         return 2; // Dih_1 x Dih_1
 264                                 default: QPMS_WTF;
 265                         }
 266         }
 267
 268         switch(cls) {
 269                 // Axial groups
 270                 case QPMS_PGS_CN:
 271                         gen[0] = qpms_irot3_zrot_Nk(n, 1);
 272                         return 1; // Z_n
 273                 case QPMS_PGS_S2N:
 274                         gen[0].rot = qpms_quat_zrot_Nk(2*n, 1);
 275                         gen[0].det = -1;
 276                         return 1; // Z_{2n}
 277                 case QPMS_PGS_CNH:
 278                         gen[0] = qpms_irot3_zrot_Nk(n, 1);
 279                         gen[1] = QPMS_IROT3_ZFLIP;
 280                         return 2; // Z_n x Dih_1
 281                 case QPMS_PGS_CNV:
 282                         gen[0] = qpms_irot3_zrot_Nk(n, 1);
 283                         gen[1] = QPMS_IROT3_XFLIP;
 284                         return 2; // Dih_n
 285                 case QPMS_PGS_DN:
 286                         gen[0] = qpms_irot3_zrot_Nk(n, 1);
 287                         gen[1] = QPMS_IROT3_XROT_PI;
 288                         return 2; // Dih_n
 289                 case QPMS_PGS_DND:
 290                         gen[0].rot = qpms_quat_zrot_Nk(2*n, 1);
 291                         gen[0].det = -1;
 292                         gen[1] = QPMS_IROT3_XFLIP;
 293                         return 2; // Dih_2n
 294                 case QPMS_PGS_DNH:
 295                         gen[0] = qpms_irot3_zrot_Nk(n, 1);
 296                         gen[1] = QPMS_IROT3_ZFLIP;
 297                         gen[2] = QPMS_IROT3_XFLIP;
 298                         return 3; // Dih_n x Dih_1
 299
 300                                    // Remaining polyhedral groups
 301                 case QPMS_PGS_T:   // return ???; // A_4
 302                 case QPMS_PGS_TD:  // return 2; // S_4
 303                 case QPMS_PGS_TH:  // A_4 x Z_2
 304                 case QPMS_PGS_O:   // S_4
 305                 case QPMS_PGS_OH:  // return 3; //  S_4 x Z_2
 306                 case QPMS_PGS_I:   // A_5
 307                 case QPMS_PGS_IH:  // A_5 x Z_2
 308
 309                                    // Continuous axial groups
 310                 case QPMS_PGS_CINF:
 311                 case QPMS_PGS_CINFH:
 312                 case QPMS_PGS_CINFV:
 313                 case QPMS_PGS_DINF:
 314                 case QPMS_PGS_DINFH:
 315
 316                                    // Remaining continuous groups
 317                 case QPMS_PGS_SO3:
 318                 case QPMS_PGS_O3:
 319                                    QPMS_NOT_IMPLEMENTED("Not yet implemented for this point group class.");
 320                 default:
 321                 QPMS_WTF;
 322         }
 323 }