chiark / gitweb /
set threads_started
[moebius2.git] / half.c
1 /*
2  * actually do mapping from metaspace to realspace
3  */
4
5 #include "half.h"
6
7 /*
8  * Our metaspace dimensions are, in order:
9  *   - y realspace coordinate of symmetry axis vertices
10  *   - x,y,z realspace coordinates of even-y opposite-of-axis vertices
11  *   - x,y,z realspace coordinates of even-y non-axis vertices
12  *   - x,y,z realspace coordinates of odd-y vertices in canonical order
13  */
14
15 #define FOR_MAP_DIMENSION__VERTEX_PAIR(body,mirror)     \
16   k=0; sign=-1; { body } { mirror } i++;                \
17   k=1; sign=+1; { body } { mirror } i++;                \
18   k=2; sign=-1; { body } { mirror } i++;
19
20 /*#define FOR_MAP_DIMENSION_VPRIME              \
21   int vprime= N-1-v
22   */
23
24 #define FOR_MAP_DIMENSION(body, mirror, degenerate)             \
25   do{                                                           \
26     int i,v,vprime,k,sign; /* public */                         \
27     int fmd_i, fmd_x, fmd_y;                                    \
28     i= 0;                                                       \
29     for (fmd_i=0; fmd_i<HALF_DEGRADED_NVERTICES; fmd_i++) {     \
30       v= fmd_i*X*2;                                             \
31       k=0;          { degenerate }                              \
32       k=1; sign=+1; { body } i++;                               \
33       k=2;          { degenerate }                              \
34     }                                                           \
35     for (fmd_y=0; fmd_y<Y; fmd_y+=2) {                          \
36       v= fmd_y << YSHIFT | (X/2);                               \
37       FOR_MAP_DIMENSION__VERTEX_PAIR(body,)                     \
38     }                                                           \
39     for (fmd_y=0; fmd_y<Y; fmd_y+=2)                            \
40       for (fmd_x=1; fmd_x<X/2; fmd_x++) {                       \
41         v= fmd_y << YSHIFT | fmd_x;                             \
42         vprime= ((Y-1) << YSHIFT | X) - v;                      \
43         FOR_MAP_DIMENSION__VERTEX_PAIR(body,mirror)             \
44       }                                                         \
45     for (fmd_y=1; fmd_y<Y; fmd_y+=2)                            \
46       for (fmd_x=0; fmd_x<X/2; fmd_x++) {                       \
47         v= fmd_y << YSHIFT | fmd_x;                             \
48         vprime= ((Y-1) << YSHIFT | (X-1)) - v;                  \
49         FOR_MAP_DIMENSION__VERTEX_PAIR(body,mirror)             \
50       }                                                         \
51     assert(i == DIM);                                           \
52   }while(0)
53
54 void pmap_dimensions(const struct Vertices *vs) {
55   int ixfrom[N][D3], count[N][D3], bad;
56   int vc,kc;
57
58   bad= 0;
59
60   FOR_VERTEX(vc,INNER)
61     FOR_COORD(kc) {
62       count[vc][kc]= 0;
63       ixfrom[vc][kc]= -1;
64     }
65   
66 #define CHECK_V(v) assert(v >= 0); assert(v < N);
67 #define CHECK_NEAR(a,b)                                 \
68   if (abs(a-b) >= 1e-6) { printf("    BAD!"); bad++; }  \
69   putchar('\n');
70
71   FOR_MAP_DIMENSION(
72                     (printf("i=%2d  v=%03x k=%d  %# 20g\n",
73                             i,v,k,
74                             vs->a[v][k]
75                             ));
76                     CHECK_V(v)
77                     ixfrom[v][k]=i;
78                     count[v][k]++;,
79
80 //                  FOR_MAP_DIMENSION_VPRIME;
81                     (printf("i=%2d v'=%03x k=%d  %# 20g := %# 20g",
82                             i,vprime,k,
83                             vs->a[vprime][k], sign*vs->a[v][k]
84                             ));
85                     CHECK_V(vprime)
86                     CHECK_NEAR(vs->a[vprime][k], sign*vs->a[v][k])
87                     ixfrom[vprime][k]=-i;
88                     count[vprime][k]++;,
89
90                     (printf("zero  v=%03x k=%d  %# 20g := 0",
91                             v,k,
92                             vs->a[v][k]
93                             ));
94                     CHECK_V(v)
95                     CHECK_NEAR(vs->a[v][k], 0)
96                     count[v][k]++;
97                     );
98
99   FOR_VERTEX(vc,INNER)
100     FOR_COORD(kc) {
101       if (count[vc][kc]==1) continue;
102       printf("BAD %03x %d count=%d ixfrom=%d\n",
103              vc,kc, count[vc][kc], ixfrom[vc][kc]);
104       bad++;
105     }
106   assert(!bad);
107 }
108
109 void map_dimensions(const double metavertex[DIM], Vertices out) {
110   FOR_MAP_DIMENSION(out[v     ][k]=      metavertex[i];,
111
112 //                  FOR_MAP_DIMENSION_VPRIME;
113                     out[vprime][k]= sign*metavertex[i];,
114
115                     out[v     ][k]=      0;
116                     );
117 }
118
119 void unmap_dimensions(double metavertex[DIM], const struct Vertices *in) {
120   FOR_MAP_DIMENSION(
121                     metavertex[i]= in->a[v][k];,
122                     ,
123                     );
124 }
125
126
127 /*
128  *  OLD NUMBERING
129  *
130  * Our numbering is as follows:
131  *
132  *          :axis of symmetry
133  *        | :
134  *        | :
135  *     ___| 0  ___  1  ___  2  ___  3  ___  4  ___    up to X"-1
136  *        | :                                           where X"=(X/2)
137  *        / :\    /  \    /  \    /  \    /  \
138  *       /| : \  /    \  /    \  /    \  /    \
139  *        |___ 0   ___ 1   ___ 2   ___ 3   ___ 4
140  *        | :  +X"     +X"     +X"     +X"     +X"
141  *       \| : /  \    /  \    /  \    /  \    /
142  *        \ :/    \  /    \  /    \  /    \  /
143  *     ___|2*X"___ 1+  ___ 2+  ___ 3+  __  4+  ___
144  *        | :     2*X"    2*X"    2*X"    2*X"
145  *        / :\    /  \    /  \    /  \    /  \
146  *       /| : \  /    \  /    \  /    \  /    \
147  *        |___ 0+  ___ 1+  ___ 2+  ___ 3+  ___ 4+
148  *        | : 3*X"    3*X"    3*X"    3*X"    3*X"
149  *        | :
150  *
151  * X" is still a power of two so a vertex pair can be named by a
152  * number structured like this:
153  *
154  * vertex number":   0000 | y     | x
155  *                         YBITS   XBITS-1
156  *
157  */