qpms/latticegens.c

   1 #include "lattices.h"
   2 #include <limits.h>
   3 #include <math.h>
   4 #include <string.h>
   5
   6 #define MIN(x,y) ((x)<(y)?(x):(y))
   7
   8 // generic converting extractors
   9
  10 PGenPolReturnData PGen_next_pol_from_cart2(PGen *g) {
  11   const PGenCart2ReturnData c = PGen_next_cart2(g);
  12   if (c.flags & PGEN_DONE)
  13     return PGenPolDoneVal;
  14   else {
  15     PGenPolReturnData p;
  16     p.flags = (c.flags & ~QPMS_COORDS_BITRANGE) | QPMS_COORDS_POL;
  17     p.point_pol = cart2pol(c.point_cart2);
  18     return p;
  19   }
  20 }
  21
  22 PGenCart2ReturnData PGen_next_cart2_from_pol(PGen *g) {
  23   const PGenPolReturnData p = PGen_next_pol(g);
  24   if (p.flags & PGEN_DONE)
  25     return PGenCart2DoneVal;
  26   else {
  27     PGenCart2ReturnData c;
  28     c.flags = (p.flags & ~QPMS_COORDS_BITRANGE) | QPMS_COORDS_CART2;
  29     c.point_cart2 = pol2cart(p.point_pol);
  30     return c;
  31   }
  32 }
  33
  34 PGenSphReturnData PGen_next_sph_from_cart3(PGen *g) {
  35   const PGenCart3ReturnData c = PGen_next_cart3(g);
  36   if (c.flags & PGEN_DONE)
  37     return PGenSphDoneVal;
  38   else {
  39     PGenSphReturnData s;
  40     s.flags = (c.flags & ~QPMS_COORDS_BITRANGE) | QPMS_COORDS_SPH;
  41     s.point_sph = cart2sph(c.point_cart3);
  42     return s;
  43   }
  44 }
  45
  46 PGenCart3ReturnData PGen_next_cart3_from_cart2xy(PGen *g) {
  47   const PGenCart2ReturnData c2 = PGen_next_cart2(g);
  48   if (c2.flags & PGEN_DONE)
  49     return PGenCart3DoneVal;
  50   else {
  51     PGenCart3ReturnData c3;
  52     c3.flags = (c2.flags & ~QPMS_COORDS_BITRANGE) | QPMS_COORDS_CART3;
  53     c3.point_cart3 = cart22cart3xy(c2.point_cart2);
  54     return c3;
  55   }
  56 }
  57
  58 PGenSphReturnData PGen_next_sph_from_cart2(PGen *g) {
  59   const PGenCart2ReturnData c = PGen_next_cart2(g);
  60   if (c.flags & PGEN_DONE)
  61     return PGenSphDoneVal;
  62   else {
  63     PGenSphReturnData s;
  64     s.flags = (c.flags & ~QPMS_COORDS_BITRANGE) | QPMS_COORDS_SPH;
  65     s.point_sph = cart22sph(c.point_cart2);
  66     return s;
  67   }
  68 }
  69
  70 PGenCart3ReturnData PGen_next_cart3_from_sph(PGen *g) {
  71   const PGenSphReturnData s = PGen_next_sph(g);
  72   if (s.flags & PGEN_DONE)
  73     return PGenCart3DoneVal;
  74   else {
  75     PGenCart3ReturnData c;
  76     c.flags = (s.flags & ~QPMS_COORDS_BITRANGE) | QPMS_COORDS_CART3;
  77     c.point_cart3 = sph2cart(s.point_sph);
  78     return c;
  79   }
  80 }
  81
  82 // here, various "classes" of the PGenSph point generators are implemented.
  83
  84
  85 // const PGenSphReturnData PGenSphDoneVal = {PGEN_DONE, {0,0,0}}; // defined already in lattices.h
  86 // const PGenCart3ReturnData PGenCart3DoneVal = {PGEN_DONE, {0,0,0}}; // defined already in lattices.h
  87
  88 // General structure of a generator implementation looks like this:
  89
  90 #if 0
  91 //==== PGen_NAME ====
  92
  93 extern const PGenClassInfo PGen_NAME; // forward declaration needed by constructor (may be placed in header file instead)
  94
  95 // Internal state structure
  96 typedef struct PGen_NAME_StateData {
  97   ...
  98 } PGen_NAME_StateData;
  99
 100 // Constructor
 101 PGenSph PGen_NAME_new(...) {
 102   g->stateData = malloc(sizeof(PGen_NAME_StateData));
 103   ...
 104   PGenSph g = {&PGen_NAME, (void *) stateData};
 105   return g;
 106 }
 107
 108 // Dectructor
 109 void PGen_NAME_destructor(PGen *g) {
 110   ...
 111   free(g->stateData);
 112   g->stateData = NULL;
 113 }
 114
 115 // Extractor, spherical coordinate output
 116 PGenSphReturnData PGen_NAME_next_sph(PGen *g) {
 117   if (g->stateData == NULL) // already destroyed
 118     return PGenSphDoneVal;
 119   else {
 120     PGen_NAME_StateData *s = (PGen_NAME_StateData *) g->stateData;
 121     if (... /* there are still points to be generated */) {
 122       ...
 123       PGenSphReturnData retval = {.../*flags*/, .../*thePoint*/};
 124       return retval;
 125     } else {
 126       PGen_destroy(g);
 127       return PGenSphDoneVal;
 128     }
 129   }
 130 }
 131
 132 // Extractor, 3D cartesian coordinate output
 133 PGenCart3ReturnData PGen_NAME_next_cart3(PGen *g) {
 134   if (g->stateData == NULL) // already destroyed
 135     return PGenCart3DoneVal;
 136   else {
 137     PGen_NAME_StateData *s = (PGen_NAME_StateData *) g->stateData;
 138     if (... /* there are still points to be generated */) {
 139       ...
 140       PGenCart3ReturnData retval = {.../*flags*/, .../*thePoint*/};
 141       return retval;
 142     } else {
 143       PGen_destroy(g);
 144       return PGenCart3DoneVal;
 145     }
 146   }
 147 }
 148
 149 // Class metadata structure; TODO maybe this can rather be done by macro.
 150 const PGenClassInfo PGen_NAME = {
 151   "PGen_NAME",
 152   ?, //dimensionality
 153   PGEN_COORDS_????, // native coordinate system
 154   // some of the _next_... fun pointers can be NULL
 155   PGen_NAME_next,
 156   PGen_NAME_next_z,
 157   PGen_NAME_next_pol,
 158   PGen_NAME_next_sph,
 159   PGen_NAME_next_cart2,
 160   PGen_NAME_next_cart3,
 161   PGen_NAME_fetch,
 162   PGen_NAME_fetch_z,
 163   PGen_NAME_fetch_pol,
 164   PGen_NAME_fetch_sph,
 165   PGen_NAME_fetch_cart2,
 166   PGen_NAME_fetch_cart3,
 167   PGen_NAME_destructor
 168 };
 169
 170 #endif // 0
 171
 172
 173 //==== PGenSph_FromPoint2DArray ====
 174
 175 // Internal state structure
 176 typedef struct PGen_FromPoint2DArray_StateData {
 177   const point2d *base;
 178   size_t len;
 179   size_t currentIndex;
 180 }PGen_FromPoint2DArray_StateData;
 181
 182 // Constructor
 183 PGen PGen_FromPoint2DArray_new(const point2d *points, size_t len) {
 184   PGen_FromPoint2DArray_StateData *stateData = malloc(sizeof(PGen_FromPoint2DArray_StateData));
 185   stateData->base = points;
 186   stateData->len = len;
 187   stateData->currentIndex = 0;
 188   PGen g = {&PGen_FromPoint2DArray, (void *) stateData};
 189   return g;
 190 }
 191
 192 // Destructor
 193 void PGen_FromPoint2DArray_destructor(PGen *g) {
 194   free(g->stateData);
 195   g->stateData = NULL;
 196 }
 197
 198 // Extractor, 2D cartesian (native)
 199 PGenCart2ReturnData PGen_FromPoint2DArray_next_cart2(PGen *g) {
 200   if (g->stateData == NULL) // already destroyed
 201     return PGenCart2DoneVal;
 202   else {
 203     PGen_FromPoint2DArray_StateData *s = (PGen_FromPoint2DArray_StateData *) g->stateData;
 204     if (s->currentIndex < s->len) {
 205       cart2_t thePoint = s->base[s->currentIndex];
 206       ++(s->currentIndex);
 207       PGenCart2ReturnData retval = {(PGEN_NOTDONE | PGEN_AT_XY | PGEN_COORDS_CART2), thePoint};
 208       return retval;
 209     } else {
 210       PGen_destroy(g);
 211       return PGenCart2DoneVal;
 212     }
 213   }
 214 }
 215
 216 // Extractor, spherical
 217 PGenSphReturnData PGen_FromPoint2DArray_next_sph(PGen *g) {
 218   if (g->stateData == NULL) // already destroyed
 219     return PGenSphDoneVal;
 220   else {
 221     PGen_FromPoint2DArray_StateData *s = (PGen_FromPoint2DArray_StateData *) g->stateData;
 222     if (s->currentIndex < s->len) {
 223       sph_t thePoint = cart22sph(s->base[s->currentIndex]);
 224       ++(s->currentIndex);
 225       PGenSphReturnData retval = {(PGEN_AT_XY | PGEN_COORDS_SPH), thePoint};
 226       return retval;
 227     } else {
 228       PGen_destroy(g);
 229       return PGenSphDoneVal;
 230     }
 231   }
 232 }
 233
 234 const PGenClassInfo PGen_FromPoint2DArray = {
 235   "PGen_FromPoint2DArray",
 236   2, // dimensionality
 237   PGEN_COORDS_CART2,
 238   NULL,//PGen_FromPoint2DArray_next,
 239   NULL,
 240   NULL,//PGen_FromPoint2DArray_next_pol,
 241   PGen_FromPoint2DArray_next_sph,
 242   PGen_FromPoint2DArray_next_cart2,
 243   NULL,//PGen_FromPoint2DArray_next_cart3,
 244   NULL,//PGen_FromPoint2DArray_fetch,
 245   NULL,//PGen_FromPoint2DArray_fetch_z,
 246   NULL,//PGen_FromPoint2DArray_fetch_pol,
 247   NULL,//PGen_FromPoint2DArray_fetch_sph,
 248   NULL,//PGen_FromPoint2DArray_fetch_cart2,
 249   NULL,//PGen_FromPoint2DArray_fetch_cart3,
 250   PGen_FromPoint2DArray_destructor,
 251 };
 252
 253
 254
 255 //==== PGen_1D ====
 256 //equidistant points along the z-axis;
 257
 258 extern const PGenClassInfo PGen_1D; // forward declaration needed by constructor (may be placed in header file instead)
 259
 260 /* // This had to go to the header file:
 261 enum PGen_1D_incrementDirection{
 262     //PGEN1D_POSITIVE_INC, // not implemented
 263     //PGEN1D_NEGATIVE_INC, // not implemented
 264     PGEN1D_INC_FROM_ORIGIN,
 265     PGEN1D_INC_TOWARDS_ORIGIN
 266 };
 267 */
 268
 269 // Internal state structure
 270 typedef struct PGen_1D_StateData {
 271   long ptindex;
 272   //long stopindex;
 273   double minR, maxR;
 274   bool inc_minR, inc_maxR;
 275   double a; // lattice period
 276   double offset; // offset of the zeroth lattice point from origin (will be normalised to interval [-a/2,a/2]
 277   enum PGen_1D_incrementDirection incdir;
 278   //bool skip_origin;
 279 } PGen_1D_StateData;
 280
 281 static inline long ptindex_inc(long i) {
 282   if (i > 0)
 283     return -i;
 284   else
 285     return -i + 1;
 286 }
 287
 288 static inline long ptindex_dec(long i) {
 289   if (i > 0)
 290     return -i + 1;
 291   else
 292     return -i;
 293 }
 294
 295 // Constructor, specified by maximum and maximum absolute value
 296 PGen PGen_1D_new_minMaxR(double period, double offset, double minR, bool inc_minR, double maxR, bool inc_maxR,
 297     PGen_1D_incrementDirection incdir) {
 298   PGen_1D_StateData *s = malloc(sizeof(PGen_1D_StateData));
 299   s->minR = minR;
 300   s->maxR = maxR;
 301   s->inc_minR = inc_minR;
 302   s->inc_maxR = inc_maxR;
 303   s->incdir = incdir;
 304   period = fabs(period);
 305   double offset_normalised = offset - period * floor(offset / period); // shift to interval [0, period]
 306   if (offset_normalised > period / 2) offset_normalised -= period; // and to interval [-period/2, period/2]
 307   s->offset = offset_normalised;
 308   if (offset_normalised > 0) // reverse the direction so that the conditions in _next() are hit in correct order
 309         period *= -1;
 310   switch(s->incdir) {
 311     double curR;
 312     case PGEN_1D_INC_FROM_ORIGIN:
 313       s->ptindex = floor(minR / fabs(period));
 314       while ( (curR = fabs(s->offset + s->ptindex * period)) < minR || (!inc_minR && curR <= minR))
 315         s->ptindex = ptindex_inc(s->ptindex);
 316       break;
 317     case PGEN_1D_INC_TOWARDS_ORIGIN:
 318       s->ptindex = - ceil(maxR / fabs(period));
 319       while ( (curR = fabs(s->offset + s->ptindex * period)) > maxR || (!inc_minR && curR >= maxR))
 320         s->ptindex = ptindex_dec(s->ptindex);
 321       break;
 322     default:
 323       QPMS_WTF;
 324   }
 325   s->a = period;
 326
 327   PGen g = {&PGen_1D, (void *) s};
 328   return g;
 329 }
 330
 331 // Dectructor
 332 void PGen_1D_destructor(PGen *g) {
 333   free(g->stateData);
 334   g->stateData = NULL;
 335 }
 336
 337 // Extractor 1D number
 338 PGenZReturnData PGen_1D_next_z(PGen *g) {
 339   if (g->stateData == NULL) // already destroyed
 340     return PGenZDoneVal;
 341   PGen_1D_StateData *s = (PGen_1D_StateData *) g->stateData;
 342   const double zval = s->ptindex * s->a + s->offset;
 343   const double r = fabs(zval);
 344   bool theEnd = false;
 345   switch (s->incdir) {
 346     case PGEN_1D_INC_FROM_ORIGIN:
 347       if (r < s->maxR || (s->inc_maxR && r == s->maxR))
 348         s->ptindex = ptindex_inc(s->ptindex);
 349       else theEnd = true;
 350       break;
 351     case PGEN_1D_INC_TOWARDS_ORIGIN:
 352       if (r > s->minR || (s->inc_minR && r == s->minR)) {
 353         if (s->ptindex == 0) // handle "underflow"
 354           s->minR = INFINITY;
 355         else
 356           s->ptindex = ptindex_dec(s->ptindex);
 357       } else theEnd = true;
 358       break;
 359     default:
 360       QPMS_INVALID_ENUM(s->incdir);
 361   }
 362   if (!theEnd) {
 363       const PGenZReturnData retval = {PGEN_NOTDONE | PGEN_AT_Z, zval};
 364       return retval;
 365   } else {
 366       PGen_destroy(g);
 367       return PGenZDoneVal;
 368   }
 369 }
 370
 371 // Extractor spherical coordinates // TODO remove/simplify
 372 PGenSphReturnData PGen_1D_next_sph(PGen *g) {
 373   if (g->stateData == NULL) // already destroyed
 374     return PGenSphDoneVal;
 375   PGen_1D_StateData *s = (PGen_1D_StateData *) g->stateData;
 376   const double zval = s->ptindex * s->a + s->offset;
 377   const double r = fabs(zval);
 378   bool theEnd = false;
 379   switch (s->incdir) {
 380     case PGEN_1D_INC_FROM_ORIGIN:
 381       if (r < s->maxR || (s->inc_maxR && r == s->maxR))
 382         s->ptindex = ptindex_inc(s->ptindex);
 383       else theEnd = true;
 384       break;
 385     case PGEN_1D_INC_TOWARDS_ORIGIN:
 386       if (r > s->minR || (s->inc_minR && r == s->minR)) {
 387         if (s->ptindex == 0) // handle "underflow"
 388           s->minR = INFINITY;
 389         else
 390           s->ptindex = ptindex_dec(s->ptindex);
 391       } else theEnd = true;
 392       break;
 393     default:
 394         QPMS_INVALID_ENUM(s->incdir);
 395   }
 396   if (!theEnd) {
 397       const PGenSphReturnData retval = {PGEN_NOTDONE | PGEN_AT_Z | PGEN_COORDS_SPH,
 398         {r, zval >= 0 ? 0 : M_PI, 0}};
 399       return retval;
 400   } else {
 401       PGen_destroy(g);
 402       return PGenSphDoneVal;
 403   }
 404 }
 405
 406 // Class metadata structure; TODO maybe this can rather be done by macro.
 407 const PGenClassInfo PGen_1D = {
 408   "PGen_1D",
 409   1, // dimensionality
 410   PGEN_COORDS_CART1,
 411   NULL, //PGen_1D_next,
 412   PGen_1D_next_z,
 413   NULL,//PGen_1D_next_pol,
 414   PGen_1D_next_sph,
 415   NULL,//PGen_1D_next_cart2,
 416   NULL,//PGen_1D_next_cart3,
 417   NULL,//PGen_1D_fetch,
 418   NULL,//PGen_1D_fetch_z,
 419   NULL,//PGen_1D_fetch_pol,
 420   NULL,//PGen_1D_fetch_sph,
 421   NULL,//PGen_1D_fetch_cart2,
 422   NULL,//PGen_1D_fetch_cart3,
 423   PGen_1D_destructor
 424 };
 425
 426 //==== PGen_xyWeb ====
 427 // 2D lattice generator in the "spiderweb" style, generated in the "perimetre" order,
 428 // not strictly ordered (or limited) by distance from origin.
 429 // The minR and maxR here refer to the TODO WWHAT
 430
 431 extern const PGenClassInfo PGen_xyWeb; // forward declaration needed by constructor (may be placed in header file instead)
 432
 433 // Internal state structure
 434 typedef struct PGen_xyWeb_StateData {
 435   long i, j;
 436   unsigned short phase; // 0 to 5
 437   long layer;
 438   long last_layer; // generation stops when layer > last_layer
 439   double layer_min_height; // this * layer is what minR and maxR are compared to
 440   double minR, maxR;
 441   bool inc_minR, inc_maxR;
 442   cart2_t b1, b2; // lattice vectors
 443   cart2_t offset; // offset of the zeroth lattice point from origin (TODO will be normalised to the WS cell)
 444   // TODO type rectangular vs. triangular
 445   LatticeFlags lf;
 446 } PGen_xyWeb_StateData;
 447
 448 // Constructor
 449 PGen PGen_xyWeb_new(cart2_t b1, cart2_t b2, double rtol, cart2_t offset, double minR, bool inc_minR, double maxR, bool inc_maxR) {
 450   PGen_xyWeb_StateData *s =  malloc(sizeof(PGen_xyWeb_StateData));
 451   s->minR = minR; s->maxR = maxR;
 452   s->inc_minR = inc_minR;
 453   s->inc_maxR = inc_maxR;
 454   l2d_reduceBasis(b1, b2, &(s->b1), &(s->b2));
 455   s->offset = offset; // TODO shorten into the WS cell ?
 456   s->lf = l2d_detectRightAngles(s->b1, s->b2, rtol);
 457   s->layer_min_height = l2d_hexWebInCircleRadius(s->b1, s->b2);
 458   s->layer = ceil(s->minR/s->layer_min_height);
 459   if(!inc_minR && (s->layer * s->layer_min_height)  <= minR)
 460     ++(s->layer);
 461   s->i = s->layer; s->j = 0; s->phase = 0; // init indices
 462   s->last_layer = floor(s->maxR/s->layer_min_height);
 463   if(!inc_maxR && (s->last_layer * s->layer_min_height) >= maxR)
 464     --(s->last_layer);
 465   PGen g = {&PGen_xyWeb, (void *) s};
 466   return g;
 467 }
 468
 469 // Destructor
 470 void PGen_xyWeb_destructor(PGen *g) {
 471   free(g->stateData);
 472   g->stateData = NULL;
 473 }
 474
 475 // Extractor (2D cartesian, native)
 476 PGenCart2ReturnData PGen_xyWeb_next_cart2(PGen *g) {
 477   if (g->stateData == NULL) // already destroyed
 478     return PGenCart2DoneVal;
 479   else {
 480     PGen_xyWeb_StateData * const s = (PGen_xyWeb_StateData *) g->stateData;
 481     assert(s->layer >= 0);
 482     if (s->layer <= s->last_layer) {
 483       const cart2_t thePoint = cart2_add(s->offset,
 484           cart2_add(cart2_scale(s->i, s->b1), cart2_scale(s->j, s->b2)));
 485       if(s->layer == 0) { // origin is unique, proceed with next layer
 486         ++s->layer;
 487         s->phase = 0;
 488         s->i = s->layer;
 489         s->j = 0;
 490       }
 491       else if(s->lf & ORTHOGONAL_01) {
 492         // rectangular or square lattice, four perimeters
 493         switch(s->phase) {
 494           case 0: // initial i = l, j = 0
 495             --s->i;
 496             ++s->j;
 497             if(s->i <= 0) ++s->phase;
 498             break;
 499           case 1: // initial i = 0, j = l
 500             --s->i;
 501             --s->j;
 502             if(s->j <= 0) ++s->phase;
 503             break;
 504           case 2: // initial i = -l, j = 0
 505             ++s->i;
 506             --s->j;
 507             if(s->i >= 0) ++s->phase;
 508             break;
 509           case 3: // initial i = 0, j = -l
 510             ++s->i;
 511             ++s->j;
 512             if(s->j >= 0) ++s->phase;
 513             break;
 514           default:
 515             QPMS_WTF;
 516         }
 517         if(s->phase == 4) { // phase overflow, start new layer
 518           ++s->layer;
 519           s->phase = 0;
 520           s->i = s->layer;
 521           s->j = 0;
 522         }
 523       } else { // non-rectangular lattice, six perimeters
 524         switch(s->phase) {
 525           case 0:
 526             --s->i;
 527             ++s->j;
 528             if(s->i <= 0) ++s->phase;
 529             break;
 530           case 1:
 531             --s->i;
 532             if(s->i + s->j <= 0) ++s->phase;
 533             break;
 534           case 2:
 535             --s->j;
 536             if(s->j <= 0) ++s->phase;
 537             break;
 538           case 3:
 539             ++s->i;
 540             --s->j;
 541             if(s->i >= 0) ++s->phase;
 542             break;
 543           case 4:
 544             ++s->i;
 545             if(s->i + s->j >= 0) ++s->phase;
 546             break;
 547           case 5:
 548             ++s->j;
 549             if(s->j >= 0) ++s->phase;
 550             break;
 551           default:
 552             QPMS_WTF;
 553         }
 554         if(s->phase == 6) { // phase overflow, start next layer
 555           ++s->layer;
 556           s->phase = 0;
 557           s->i = s->layer;
 558           s->j = 0;
 559         }
 560       }
 561       PGenCart2ReturnData retval = {(PGEN_NOTDONE | PGEN_AT_XY | PGEN_COORDS_CART2), thePoint};
 562       return retval;
 563     } else {
 564       PGen_destroy(g);
 565       return PGenCart2DoneVal;
 566     }
 567   }
 568 }
 569
 570 // Class metadata structure; TODO maybe this can rather be done by macro.
 571 const PGenClassInfo PGen_xyWeb = {
 572   "PGen_xyWeb",
 573   2,
 574   PGEN_COORDS_CART2,
 575   NULL,//PGen_xyWeb_next, // FIXME I should really implement this.
 576   NULL,//PGen_xyWeb_next_z,
 577   PGen_next_pol_from_cart2, //NULL,//PGen_xyWeb_next_pol,
 578   PGen_next_sph_from_cart2, //NULL,//PGen_xyWeb_next_sph,
 579   PGen_xyWeb_next_cart2, // native
 580   PGen_next_cart3_from_cart2xy, //NULL,//PGen_xyWeb_next_cart3,
 581   NULL,//PGen_xyWeb_fetch, // FIXME I should really implement this
 582   NULL,//PGen_xyWeb_fetch_z,
 583   NULL,//PGen_xyWeb_fetch_pol,
 584   NULL,//PGen_xyWeb_fetch_sph,
 585   NULL,//PGen_xyWeb_fetch_cart2,
 586   NULL,//PGen_xyWeb_fetch_cart3,
 587   PGen_xyWeb_destructor
 588 };
 589
 590
 591 size_t PGen_xyWeb_sizecap(cart2_t b1, cart2_t b2, double rtol, cart2_t offset,
 592     double minR, bool inc_minR, double maxR, bool inc_maxR)
 593 {
 594   l2d_reduceBasis(b1, b2, &b1, &b2);
 595   LatticeFlags lf = l2d_detectRightAngles(b1, b2, rtol);
 596   double layer_min_height = l2d_hexWebInCircleRadius(b1, b2);
 597   long layer = ceil(minR / layer_min_height);
 598   if(!inc_minR && (layer * layer_min_height)  <= minR)
 599     ++layer;
 600   long last_layer = floor(maxR / layer_min_height);
 601   if(!inc_maxR && (last_layer * layer_min_height) >= maxR)
 602     --(last_layer);
 603   // TODO less crude estimate (this one should be safe, however)
 604   return ((lf & ORTHOGONAL_01) ? 4 : 6) * (last_layer - layer + 1);
 605 }
 606
 607
 608 //==== PGen_LatticeRadialHeap ====
 609
 610
 611 extern const PGenClassInfo PGen_LatticeRadialHeap; // forward declaration needed by constructor (may be placed in header file instead)
 612
 613 // Internal state structure
 614 typedef struct PGen_LatticeRadialHeap_StateData {
 615   int ldim; // 2 or 3, must
 616   int sdim;
 617   int layer;
 618   // Minimal distance of a point from the origin in the last layer (if the top of the heap exceeds this, a new layer must be evaluated
 619   double layer_min_r;
 620   size_t heap_len;
 621   size_t heap_capacity;
 622   double *r_heap;
 623   int *coord_heap;
 624   double b[0]; // basis vectors and offset
 625   double offset_r;
 626   double minR, maxR;
 627   bool inc_minR, inc_maxR;
 628 } PGen_LatticeRadialHeap_StateData;
 629
 630 static inline double nd2norm(const double a[], int d) {
 631   double n = 0;
 632   for(int i = 0; i < d; ++i)
 633     n += a[i]*a[i];
 634   return sqrt(n);
 635 }
 636
 637 // General Constructor
 638 PGen_LatticeRadialHeap_StateData *PGen_LatticeRadialHeap_new(int ldim, int sdim, double bvectors[],
 639     double offset[], double minR, double maxR,
 640     bool inc_minR, bool inc_maxR) {
 641   QPMS_UNTESTED;
 642   PGen_LatticeRadialHeap_StateData *s =
 643        malloc(sizeof(PGen_LatticeRadialHeap_StateData) + (ldim + 1) * sdim * sizeof(double));
 644   s->ldim = ldim;
 645   s->sdim = sdim;
 646   memcpy(s->b, bvectors, ldim * sdim * sizeof(double));
 647   if (offset) {
 648     memcpy(s->b + ldim * sdim, offset, sdim * sizeof(double));
 649     s->offset_r = nd2norm(s->b + ldim*sdim, sdim);
 650   } else {
 651     for (size_t i = 0; i < sdim; ++i)
 652       s->b[ldim*sdim + i] = 0;
 653     s->offset_r = 0;
 654   }
 655   s->heap_len = 0;
 656   s->heap_capacity = 1024;
 657   QPMS_CRASHING_MALLOC(s->r_heap, sizeof(*s->r_heap) * s->heap_capacity);
 658   QPMS_CRASHING_MALLOC(s->coord_heap, sizeof(*s->coord_heap) * ldim * s->heap_capacity);
 659   s->layer = -1;
 660   s->layer_min_r = -INFINITY;
 661   s->inc_minR = inc_minR; s->inc_maxR = inc_maxR;
 662
 663   return s;
 664 }
 665
 666 // 2D constructor
 667 PGen PGen_LatticeRadialHeap2D_new(cart2_t b1, cart2_t b2, cart2_t offset,
 668     double minR, bool inc_minR, double maxR, bool inc_maxR) {
 669   double bvectors[4];
 670   double offset_a[2];
 671   cart2_to_double_array(&bvectors[0], b1);
 672   cart2_to_double_array(&bvectors[2], b2);
 673   cart2_to_double_array(offset_a, offset);
 674   PGen g = {.c = &PGen_LatticeRadialHeap2D,
 675     .stateData = PGen_LatticeRadialHeap_new(2, 2, bvectors, offset_a,
 676         minR, maxR, inc_minR, inc_maxR)};
 677   return g;
 678 }
 679
 680 // 3D constructor
 681 PGen PGen_LatticeRadialHeap3D_new(const cart3_t *b1, const cart3_t *b2, const cart3_t *b3,
 682                 const cart3_t *offset, double minR, bool inc_minR, double maxR, bool inc_maxR) {
 683   double bvectors[9];
 684   double offset_a[3];
 685   int ldim = 0;
 686   if (b1) {cart3_to_double_array(&bvectors[3*ldim], *b1); ++ldim;}
 687   if (b2) {cart3_to_double_array(&bvectors[3*ldim], *b2); ++ldim;}
 688   if (b3) {cart3_to_double_array(&bvectors[3*ldim], *b3); ++ldim;}
 689   QPMS_ENSURE(ldim > 0, "At least one basis vector must be specified (non-NULL)");
 690   PGen g = {.c = &PGen_LatticeRadialHeap3D,
 691     .stateData = PGen_LatticeRadialHeap_new(ldim, 3, bvectors,
 692        offset ? (cart3_to_double_array(offset_a, *offset), offset_a) : NULL,
 693        minR, maxR, inc_minR, inc_maxR)};
 694   return g;
 695 }
 696
 697 // Increment a counter array with constant sum; in the beginning, both counter[] and counter_cumsum[]
 698 // are expected to be init'd to {thesum, 0, 0, ..., 0}
 699 // Everything is expected to be positive (although the types are signed)
 700 static inline _Bool counter_increment(const int ldim, int counter[],
 701     int counter_cumsum[], const int thesum) {
 702   for(int i = 1; i < ldim; ++i) {
 703     if(counter_cumsum[i] < thesum) { // Can increment this, do it.
 704       ++counter[i];
 705       ++counter_cumsum[i];
 706       for(int j = i - 1; j > 0; --j) { // Zero the preceding ones
 707         counter[j] = 0;
 708         counter_cumsum[j] = counter_cumsum[i];
 709       }
 710       // Determine the last digit
 711       counter[0] = thesum - counter_cumsum[1];
 712       return 1; // Incremented successfully
 713     }
 714   }
 715   return 0; // Could not increment (counter exhausted)
 716 }
 717
 718 // Assuming the counter is initialised to all non-negative values, step through
 719 // all the possible sign combinations. In the end, return them back to non-negative
 720 // and return false.
 721 static inline _Bool counter_signcycle(const int ldim, int counter[]) {
 722   for(int i = 0; i < ldim; ++i) { // Flip signs until a sign flipped to negative (cool, right?)
 723     counter[i] = -counter[i];
 724     if (counter[i] < 0)
 725       return 1;
 726   }
 727   return 0;
 728 }
 729
 730 static inline double PGen_LatticeRadialHeap_nextlayer(PGen_LatticeRadialHeap_StateData *s) {
 731   double mindist = 0;
 732   s->layer++;
 733   double minr = +INFINITY;
 734   const int ldim = s->ldim;
 735   const int sdim = s->sdim;
 736   int *counter, *counter_cumsum, *tmp;
 737   QPMS_CRASHING_CALLOC(counter, s->ldim * 2,  sizeof(*counter));
 738   counter_cumsum = counter + s->ldim;
 739   counter[0] = s->layer;
 740   counter_cumsum[0] = s->layer;
 741
 742   do {
 743     if (s->heap_len >= s->heap_capacity - 1) { // Check heap capacity
 744       s->heap_capacity = s->heap_capacity >= 1024 ? s->heap_capacity * 2 : 1024;
 745       QPMS_CRASHING_REALLOC(s->r_heap, sizeof(*s->r_heap) * s->heap_capacity);
 746       QPMS_CRASHING_REALLOC(s->coord_heap, sizeof(*s->coord_heap) * ldim * s->heap_capacity);
 747     }
 748     double r = 0;
 749     for(int i = 0; i < s->sdim; ++i) { // calculate r
 750       double component = s->b[ldim * sdim + i]; // offset
 751       for (int j = 0; j < s->ldim; ++j)
 752         component += s->b[j * sdim + i] * counter[j];
 753       r += component * component;
 754     }
 755     r = sqrt(r);
 756     minr = MIN(r, minr);
 757     // Add to the heaps
 758     int position = s->heap_len++;
 759     s->r_heap[position] = r;
 760     memcpy(&s->coord_heap[ldim * position], counter, sizeof(*s->coord_heap) * ldim);
 761     // bubble-up
 762     while(position > 0) {
 763       int parent = (position - 1) / 2;
 764       if  (s->r_heap[parent] > r)  { // swap
 765         s->r_heap[position] = s->r_heap[parent];
 766         memcpy(&s->coord_heap[ldim * position], &s->coord_heap[ldim * parent], sizeof(*s->coord_heap) * ldim);
 767         s->r_heap[parent] = r;
 768         memcpy(&s->coord_heap[ldim * parent], counter, sizeof(*s->coord_heap) * ldim);
 769         position = parent;
 770       }
 771       else break;
 772     }
 773   } while (counter_signcycle(s->ldim, counter)
 774       || counter_increment(s->ldim, counter, counter_cumsum, s->layer));
 775
 776   free(counter);
 777   return minr;
 778 }
 779
 780
 781
 782 // sdim-independent generator method, returns the integer lattice coordinates
 783 // N.B. This ALWAYS produces, not checking against maxR or destructing the generator itself
 784 // (although it does discard the points with distance smaller (or equal) than minR)
 785 // TODO maybe I want to return (double) r instead of (int) 0
 786 int PGen_LatticeRadialHeap_fillNext_intcoords(PGen_LatticeRadialHeap_StateData *s, int target[]) {
 787   bool hit = false;
 788   while(!hit) {
 789     // Ensure that we have sufficiently filled heap
 790     while (s->heap_len < 1 || s->r_heap[0] + s->offset_r > s->layer_min_r)
 791       s->layer_min_r = PGen_LatticeRadialHeap_nextlayer(s);
 792     double r = s->r_heap[0];
 793     hit = (r > s->minR || (s->inc_minR && r == s->minR));
 794     if (hit) memcpy(target, s->coord_heap, s->ldim * sizeof(*target));
 795     // Heap extract anyway
 796     // Move last element to root
 797     --(s->heap_len);
 798     s->r_heap[0] = s->r_heap[s->heap_len];
 799     memmove(s->coord_heap, &s->coord_heap[s->ldim * s->heap_len], s->ldim * sizeof(*s->coord_heap));
 800     // Bubble down
 801     int pos = 0;
 802     while(1) {
 803       int largest = pos, kidL = 2*pos+1, kidR = 2*pos+2;
 804       if (kidL < s->heap_len && s->r_heap[kidL] > s->r_heap[largest])
 805         largest = kidL;
 806       if (kidR < s->heap_len && s->r_heap[kidR] > s->r_heap[largest])
 807         largest = kidR;
 808       if (largest == pos)
 809         break;
 810       else { // swap
 811         s->r_heap[pos] = s->r_heap[largest];
 812         s->r_heap[largest] = r;
 813         memcpy(&s->coord_heap[s->ldim * pos], &s->coord_heap[s->ldim * largest], s->ldim * sizeof(*s->coord_heap));
 814         memcpy(&s->coord_heap[s->ldim * largest], &s->coord_heap[s->ldim * s->heap_len /*it's still there*/],
 815             s->ldim * sizeof(*s->coord_heap));
 816       }
 817     }
 818   }
 819   return 0;
 820 }
 821
 822 // sdim-independent generator method, fills the point's cartesian coordinates, including the offset
 823 // N.B. This ALWAYS produces, returning 0 if maxR not exceeded, else -2
 824 // (it does discard the points with distance smaller (or equal) than minR)
 825 #define LRH_STACKBUFSIZ 3 // Avoid heap allocation for typical dimensions
 826 int PGen_LatticeRadialHeap_fillNext_cart(PGen *g, double target[]) {
 827   if (g->stateData == NULL) // already destroyed
 828     return -1; // LPTODO some other return value?
 829   else {
 830     PGen_LatticeRadialHeap_StateData * const s = g->stateData;
 831     int stackbuf[LRH_STACKBUFSIZ];
 832     int *buf;
 833     if (s->sdim > LRH_STACKBUFSIZ) {
 834       QPMS_CRASHING_MALLOC(buf, s->sdim * sizeof(*buf));
 835     } else
 836       buf = stackbuf;
 837     QPMS_ENSURE_SUCCESS(PGen_LatticeRadialHeap_fillNext_intcoords(s, buf));
 838     // Calculate the actual point
 839     double r = 0;
 840     for(int i = 0; i < s->sdim; ++i) { // calculate r
 841       double component = s->b[s->ldim * s->sdim + i]; // offset
 842       for (int j = 0; j < s->ldim; ++j)
 843         component += s->b[j * s->sdim + i] * buf[j];
 844       r += component * component;
 845       target[i] = component;
 846     }
 847     r = sqrt(r);
 848     if (s->sdim > LRH_STACKBUFSIZ) free(buf);
 849
 850     if (r < s->maxR || (r == s->maxR && s->inc_maxR))
 851       return 0;
 852     else
 853       return -2;
 854   }
 855 }
 856
 857 // Destructor
 858 void PGen_LatticeRadialHeap_destructor(PGen *g) {
 859   PGen_LatticeRadialHeap_StateData *s = g->stateData;
 860   free(s->r_heap);
 861   free(s->coord_heap);
 862   free(g->stateData);
 863   g->stateData = NULL;
 864 }
 865
 866 // Extractor, 2D cartesian coordinate output
 867 PGenCart2ReturnData PGen_LatticeRadialHeap2D_next_cart2(PGen *g) {
 868   if (g->stateData == NULL) // already destroyed
 869     return PGenCart2DoneVal;
 870   else {
 871     PGen_LatticeRadialHeap_StateData *s = (PGen_LatticeRadialHeap_StateData *) g->stateData;
 872     QPMS_ENSURE(s->sdim <= 2, "Attemted to get a 2D point from %nD generator.", (int)(s->sdim));
 873     double target[2] = {0,0};
 874     if (PGen_LatticeRadialHeap_fillNext_cart(g, target) == 0 /* there are still points to be generated */) {
 875       PGenCart2ReturnData retval = {PGEN_COORDS_CART2 | PGEN_NOTDONE, cart2_from_double_array(target)};
 876       return retval;
 877     } else { // Maybe more checking for errors?
 878       PGen_destroy(g);
 879       return PGenCart2DoneVal;
 880     }
 881   }
 882 }
 883
 884 // Extractor, 3D cartesian coordinate output
 885 PGenCart3ReturnData PGen_LatticeRadialHeap3D_next_cart3(PGen *g) {
 886   if (g->stateData == NULL) // already destroyed
 887     return PGenCart3DoneVal;
 888   else {
 889     PGen_LatticeRadialHeap_StateData *s = (PGen_LatticeRadialHeap_StateData *) g->stateData;
 890     QPMS_ENSURE(s->sdim <= 3, "Attemted to get a 3D point from %nD generator.", (int)(s->sdim));
 891     double target[3] = {0,0,0};
 892     if (PGen_LatticeRadialHeap_fillNext_cart(g, target) == 0 /* there are still points to be generated */) {
 893       PGenCart3ReturnData retval = {PGEN_COORDS_CART3 | PGEN_NOTDONE, cart3_from_double_array(target)};
 894       return retval;
 895     } else { // Maybe more checking for errors?
 896       PGen_destroy(g);
 897       return PGenCart3DoneVal;
 898     }
 899   }
 900 }
 901
 902 // Class metadata structure; TODO maybe this can rather be done by macro.
 903 const PGenClassInfo PGen_LatticeRadialHeap2D = {
 904   "PGen_LatticeRadialHeap2D", // 2D labels the space, not lattice dimension
 905   2, //dimensionality
 906   PGEN_COORDS_CART2, // native coordinate system
 907   // some of the _next_... fun pointers can be NULL
 908   NULL, //PGen_LatticeRadialHeap2D_next,
 909   NULL, //PGen_LatticeRadialHeap2D_next_z,
 910   PGen_next_pol_from_cart2, // NULL, //PGen_LatticeRadialHeap2D_next_pol,
 911   PGen_next_sph_from_cart2, // NULL, //PGen_LatticeRadialHeap2D_next_sph,
 912   PGen_LatticeRadialHeap2D_next_cart2,
 913   PGen_next_cart3_from_cart2xy, // NULL, //PGen_LatticeRadialHeap2D_next_cart3,
 914   NULL, //PGen_LatticeRadialHeap2D_fetch,
 915   NULL, //PGen_LatticeRadialHeap2D_fetch_z,
 916   NULL, //PGen_LatticeRadialHeap2D_fetch_pol,
 917   NULL, //PGen_LatticeRadialHeap2D_fetch_sph,
 918   NULL, //PGen_LatticeRadialHeap2D_fetch_cart2,
 919   NULL, //PGen_LatticeRadialHeap2D_fetch_cart3,
 920   PGen_LatticeRadialHeap_destructor
 921 };
 922
 923 // Class metadata structure; TODO maybe this can rather be done by macro.
 924 const PGenClassInfo PGen_LatticeRadialHeap3D = {
 925   "PGen_LatticeRadialHeap3D", // 3D labels the space, not lattice dimension
 926   3, //dimensionality
 927   PGEN_COORDS_CART3, // native coordinate system
 928   // some of the _next_... fun pointers can be NULL
 929   NULL, //PGen_LatticeRadialHeap3D_next,
 930   NULL, //PGen_LatticeRadialHeap3D_next_z,
 931   NULL, //PGen_LatticeRadialHeap3D_next_pol,
 932   PGen_next_sph_from_cart3, // NULL, //PGen_LatticeRadialHeap3D_next_sph,
 933   NULL, //PGen_LatticeRadialHeap3D_next_cart2,
 934   PGen_LatticeRadialHeap3D_next_cart3,
 935   NULL, //PGen_LatticeRadialHeap3D_fetch,
 936   NULL, //PGen_LatticeRadialHeap3D_fetch_z,
 937   NULL, //PGen_LatticeRadialHeap3D_fetch_pol,
 938   NULL, //PGen_LatticeRadialHeap3D_fetch_sph,
 939   NULL, //PGen_LatticeRadialHeap3D_fetch_cart2,
 940   NULL, //PGen_LatticeRadialHeap3D_fetch_cart3,
 941   PGen_LatticeRadialHeap_destructor
 942 };
 943
 944