Some 3D lattice functionality, little refactoring.

Former-commit-id: 8e6572c7352e3b9a75bc13585d688f1c72d378e1

qpms/latticegens.c

@@ -439,7 +439,7 @@ typedef struct PGen_xyWeb_StateData {
   cart2_t b1, b2; // lattice vectors
   cart2_t offset; // offset of the zeroth lattice point from origin (TODO will be normalised to the WS cell)
   // TODO type rectangular vs. triangular
-  LatticeType2 lt;
+  LatticeFlags lf;
 } PGen_xyWeb_StateData;
 
 // Constructor
@@ -450,7 +450,7 @@ PGen PGen_xyWeb_new(cart2_t b1, cart2_t b2, double rtol, cart2_t offset, double
   s->inc_maxR = inc_maxR;
   l2d_reduceBasis(b1, b2, &(s->b1), &(s->b2));
   s->offset = offset; // TODO shorten into the WS cell ?
-  s->lt = l2d_classifyLattice(s->b1, s->b2, rtol);
+  s->lf = l2d_detectRightAngles(s->b1, s->b2, rtol);
   s->layer_min_height = l2d_hexWebInCircleRadius(s->b1, s->b2);
   s->layer = ceil(s->minR/s->layer_min_height);
   if(!inc_minR && (s->layer * s->layer_min_height)  <= minR)
@@ -485,7 +485,7 @@ PGenCart2ReturnData PGen_xyWeb_next_cart2(PGen *g) {
         s->i = s->layer;
         s->j = 0;
       }
-      else if(s->lt & (SQUARE | RECTANGULAR)) {
+      else if(s->lf & ORTHOGONAL_01) {
         // rectangular or square lattice, four perimeters
         switch(s->phase) {
           case 0: // initial i = l, j = 0

qpms/lattices.h

@@ -73,6 +73,9 @@ static inline point2d point2d_fromxy(const double x, const double y) {
 int qpms_reduce_lattice_basis(double *b, ///< Array of dimension [bsize][ndim].
 		const size_t bsize, ///< Number of the basis vectors (dimensionality of the lattice).
 		const size_t ndim, ///< Dimension of the space into which the lattice is embedded.
+		/// Lovász condition parameter \f$ \delta \f$.
+		/** Polynomial time complexity guaranteed for \f$\delta \in (1/4,1)\f$.
+		 */
 		double delta
 		);
 
@@ -655,12 +658,25 @@ typedef enum {
 } SpaceGroup2;
 #endif
 
+// Just for detecting the right angles (needed for generators).
+typedef enum {
+	NOT_ORTHOGONAL = 0,
+	ORTHOGONAL_01 = 1,
+	ORTHOGONAL_12 = 2,
+	ORTHOGONAL_02 = 4
+} LatticeFlags;
+
+
+
 /*
  * Lagrange-Gauss reduction of a 2D basis.
  * The output shall satisfy |out1| <= |out2| <= |out2 - out1|
  */
 void l2d_reduceBasis(cart2_t in1, cart2_t in2, cart2_t *out1, cart2_t *out2);
 
+// This one uses LLL reduction.
+void l3d_reduceBasis(const cart3_t in[3], cart3_t out[3]);
+
 /* 
  * This gives the "ordered shortest triple" of base vectors (each pair from the triple
  * is a base) and there may not be obtuse angle between o1, o2 and between o2, o3
@@ -684,6 +700,10 @@ bool l2d_is_obtuse(cart2_t i1, cart2_t i2);
  */
 LatticeType2 l2d_classifyLattice(cart2_t b1, cart2_t b2, double rtol);
 
+// Detects right angles.
+LatticeFlags l2d_detectRightAngles(cart2_t b1, cart2_t b2, double rtol);
+LatticeFlags l3d_detectRightAngles(const cart3_t basis[3], double rtol);
+
 // Other functions in lattices2d.py: TODO?
 // range2D()
 // generateLattice()

qpms/lattices2d.c

@@ -677,6 +677,11 @@ void l2d_reduceBasis(cart2_t b1, cart2_t b2, cart2_t *out1, cart2_t *out2){
   *out2 = b2;
 }
 
+void l3d_reduceBasis(const cart3_t in[3], cart3_t out[3]) {
+  memcpy(out, in, 3*sizeof(cart3_t));
+  QPMS_ENSURE_SUCCESS(qpms_reduce_lattice_basis((double *)out, 3, 3, 1.));
+}
+
 /*
  * This gives the "ordered shortest triple" of base vectors (each pair from the triple
  * is a base) and there may not be obtuse angle between o1, o2 and between o2, o3
@@ -765,6 +770,33 @@ LatticeType2 l2d_classifyLattice(cart2_t b1, cart2_t b2, double rtol)
     return OBLIQUE;
 }
 
+LatticeFlags l2d_detectRightAngles(cart2_t b1, cart2_t b2, double rtol)
+{
+  l2d_reduceBasis(b1, b2, &b1, &b2);
+  cart2_t ht = cart2_substract(b2, b1);
+  double b1s = cart2_normsq(b1), b2s = cart2_normsq(b2), hts = cart2_normsq(ht);
+  double eps = rtol * (b2s + b1s); //FIXME what should eps be?
+  if (hts - b2s - b1s <= eps)
+    return ORTHOGONAL_01;
+  else 
+    return NOT_ORTHOGONAL;
+}
+
+LatticeFlags l3d_detectRightAngles(const cart3_t basis_nr[3], double rtol) 
+{
+  cart3_t b[3];
+  l3d_reduceBasis(basis_nr, b);
+  LatticeFlags res = NOT_ORTHOGONAL;
+  for (int i = 0; i < 3; ++i) {
+    cart3_t ba = b[i], bb = b[(i+1) % 3];
+    cart3_t ht = cart3_substract(ba, bb);
+    double bas = cart3_normsq(ba), bbs = cart3_normsq(ba), hts = cart3_normsq(ht);
+    double eps = rtol * (bas + bbs);
+    if (hts - bbs - bas <= eps)
+      res |= ((LatticeFlags[]){ORTHOGONAL_01, ORTHOGONAL_12, ORTHOGONAL_02})[i];
+  }
+  return res;
+}
 
 # if 0
 // variant

qpms/vectors.h

@@ -99,11 +99,20 @@ static inline cart2_t cart3xy2cart2(const cart3_t a) {
 	return c;
 }
 
-/// 3D vector euclidian norm.
+/// 3D vector dot product.
-static inline double cart3norm(const cart3_t v) {
+static inline double cart3_dot(const cart3_t a, const cart3_t b) {
-	return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
+	return a.x * b.x + a.y * b.y + a.z * b.z;
 }
 
+/// 3D vector euclidian norm squared.
+static inline double cart3_normsq(const cart3_t a) {
+	return cart3_dot(a, a);
+}
+
+/// 3D vector euclidian norm.
+static inline double cart3norm(const cart3_t v) {
+	return sqrt(cart3_normsq(v));
+}
 
 /// 3D cartesian to spherical coordinates conversion. See @ref coord_conversions.
 static inline sph_t cart2sph(const cart3_t cart) {
@@ -224,11 +233,6 @@ static inline complex double csphvec_dotnc(const csphvec_t a, const csphvec_t b)
 	return a.rc * b.rc + a.thetac * b.thetac + a.phic * b.phic;
 }
 
-/// 3D vector dot product.
-static inline double cart3_dot(const cart3_t a, const cart3_t b) {
-	return a.x * b.x + a.y * b.y + a.z * b.z;
-}
-
 /// Spherical coordinate system scaling.
 static inline sph_t sph_scale(double c, const sph_t s) {
 	sph_t res = {c * s.r, s.theta, s.phi};