exploit symmetry

[moebius2.git] / half.c

diff --git a/half.c b/half.c
new file mode 100644 (file)
index 0000000..5ac69f8
--- /dev/null
+++ b/half.c
@@ -0,0 +1,157 @@
+/*
+ * actually do mapping from metaspace to realspace
+ */
+
+#include "half.h"
+
+/*
+ * Our metaspace dimensions are, in order:
+ *   - y realspace coordinate of symmetry axis vertices
+ *   - x,y,z realspace coordinates of even-y opposite-of-axis vertices
+ *   - x,y,z realspace coordinates of even-y non-axis vertices
+ *   - x,y,z realspace coordinates of odd-y vertices in canonical order
+ */
+
+#define FOR_MAP_DIMENSION__VERTEX_PAIR(body,mirror)    \
+  k=0; sign=-1; { body } { mirror } i++;               \
+  k=1; sign=+1; { body } { mirror } i++;               \
+  k=2; sign=-1; { body } { mirror } i++;
+
+/*#define FOR_MAP_DIMENSION_VPRIME             \
+  int vprime= N-1-v
+  */
+
+#define FOR_MAP_DIMENSION(body, mirror, degenerate)            \
+  do{                                                          \
+    int i,v,vprime,k,sign; /* public */                                \
+    int fmd_i, fmd_x, fmd_y;                                   \
+    i= 0;                                                      \
+    for (fmd_i=0; fmd_i<HALF_DEGRADED_NVERTICES; fmd_i++) {    \
+      v= fmd_i*X*2;                                            \
+      k=0;          { degenerate }                             \
+      k=1; sign=+1; { body } i++;                              \
+      k=2;          { degenerate }                             \
+    }                                                          \
+    for (fmd_y=0; fmd_y<Y; fmd_y+=2) {                         \
+      v= fmd_y << YSHIFT | (X/2);                              \
+      FOR_MAP_DIMENSION__VERTEX_PAIR(body,)                    \
+    }                                                          \
+    for (fmd_y=0; fmd_y<Y; fmd_y+=2)                           \
+      for (fmd_x=1; fmd_x<X/2; fmd_x++) {                      \
+       v= fmd_y << YSHIFT | fmd_x;                             \
+        vprime= ((Y-1) << YSHIFT | X) - v;                     \
+        FOR_MAP_DIMENSION__VERTEX_PAIR(body,mirror)            \
+      }                                                                \
+    for (fmd_y=1; fmd_y<Y; fmd_y+=2)                           \
+      for (fmd_x=0; fmd_x<X/2; fmd_x++)        {                       \
+       v= fmd_y << YSHIFT | fmd_x;                             \
+        vprime= ((Y-1) << YSHIFT | (X-1)) - v;                 \
+        FOR_MAP_DIMENSION__VERTEX_PAIR(body,mirror)            \
+      }                                                                \
+    assert(i == DIM);                                          \
+  }while(0)
+
+void pmap_dimensions(const struct Vertices *vs) {
+  int ixfrom[N][D3], count[N][D3], bad;
+  int vc,kc;
+
+  bad= 0;
+
+  FOR_VERTEX(vc)
+    FOR_COORD(kc) {
+      count[vc][kc]= 0;
+      ixfrom[vc][kc]= -1;
+    }
+  
+#define CHECK_V(v) assert(v >= 0); assert(v < N);
+#define CHECK_NEAR(a,b)                                        \
+  if (abs(a-b) >= 1e-6) { printf("    BAD!"); bad++; } \
+  putchar('\n');
+
+  FOR_MAP_DIMENSION(
+                   (printf("i=%2d  v=%03x k=%d  %# 20g\n",
+                           i,v,k,
+                           vs->a[v][k]
+                           ));
+                   CHECK_V(v)
+                   ixfrom[v][k]=i;
+                   count[v][k]++;,
+
+//                 FOR_MAP_DIMENSION_VPRIME;
+                   (printf("i=%2d v'=%03x k=%d  %# 20g := %# 20g",
+                           i,vprime,k,
+                           vs->a[vprime][k], sign*vs->a[v][k]
+                           ));
+                   CHECK_V(vprime)
+                   CHECK_NEAR(vs->a[vprime][k], sign*vs->a[v][k])
+                   ixfrom[vprime][k]=-i;
+                   count[vprime][k]++;,
+
+                   (printf("zero  v=%03x k=%d  %# 20g := 0",
+                           v,k,
+                           vs->a[v][k]
+                           ));
+                   CHECK_V(v)
+                   CHECK_NEAR(vs->a[v][k], 0)
+                   count[v][k]++;
+                   );
+
+  FOR_VERTEX(vc)
+    FOR_COORD(kc) {
+      if (count[vc][kc]==1) continue;
+      printf("BAD %03x %d count=%d ixfrom=%d\n",
+            vc,kc, count[vc][kc], ixfrom[vc][kc]);
+      bad++;
+    }
+  assert(!bad);
+}
+
+void map_dimensions(const double metavertex[DIM], Vertices out) {
+  FOR_MAP_DIMENSION(out[v     ][k]=      metavertex[i];,
+
+//                 FOR_MAP_DIMENSION_VPRIME;
+                   out[vprime][k]= sign*metavertex[i];,
+
+                   out[v     ][k]=      0;
+                   );
+}
+
+void unmap_dimensions(double metavertex[DIM], const struct Vertices *in) {
+  FOR_MAP_DIMENSION(
+                   metavertex[i]= in->a[v][k];,
+                   ,
+                   );
+}
+
+
+/*
+ *  OLD NUMBERING
+ *
+ * Our numbering is as follows:
+ *
+ *          :axis of symmetry
+ *               | :
+ *               | :
+ *     ___| 0  ___  1  ___  2  ___  3  ___  4  ___    up to X"-1
+ *        | :                                           where X"=(X/2)
+ *        / :\    /  \    /  \    /  \    /  \
+ *       /| : \  /    \  /    \  /    \  /    \
+ *        |___ 0   ___ 1   ___ 2   ___ 3   ___ 4
+ *        | :  +X"     +X"     +X"     +X"     +X"
+ *       \| : /  \    /  \    /  \    /  \    /
+ *        \ :/    \  /    \  /    \  /    \  /
+ *     ___|2*X"___ 1+  ___ 2+  ___ 3+  __  4+  ___
+ *        | :     2*X"   2*X"    2*X"    2*X"
+ *        / :\    /  \    /  \    /  \    /  \
+ *       /| : \  /    \  /    \  /    \  /    \
+ *        |___ 0+  ___ 1+  ___ 2+  ___ 3+  ___ 4+
+ *        | : 3*X"    3*X"    3*X"    3*X"    3*X"
+ *        | :
+ *
+ * X" is still a power of two so a vertex pair can be named by a
+ * number structured like this:
+ *
+ * vertex number":   0000 | y     | x
+ *                         YBITS   XBITS-1
+ *
+ */