X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?a=blobdiff_plain;f=anneal.c;h=fd725f1844d577488a2e9d154fc03c155f9bf913;hb=8ff2e481a33d24276b85600c4f3fd7bc9cbf5a18;hp=4ddddb8934901dce708ad7270094365eaae2f169;hpb=259977336c8b6a3fbd6543927d18b840b8511bfc;p=moebius2.git

diff --git a/anneal.c b/anneal.c
index 4ddddb8..fd725f1 100644
--- a/anneal.c
+++ b/anneal.c
@@ -1,70 +1,7 @@
 /*
  * We try to find an optimal triangle grid
- *
- * Vertices in strip are numbered as follows:
- *
- *     ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3   ___ 4   __
- *         Y-1     Y-1      0       0       0       0       0
- *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
- *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
- *     X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___  3  ___ 4
- *     Y-2     Y-2     Y-2      1       1       1       1       1
- *       \    /  \    /  \    /  \    /  \    /  \    /  \    /
- *        \  /    \  /    \  /    \  /    \  /    \  /    \  /
- *     ___ X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3   __
- *         Y-3     Y-3     Y-3      2       2       2       2
- *
- *       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
- *
- *      X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   ___ 3
- *       1       1       1       1      Y-2     Y-2     Y-2     Y-2
- *        \    /  \    /  \    /  \    /  \    /  \    /  \    /
- *         \  /    \  /    \  /    \  /    \  /    \  /    \  /
- *      ___ X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0   ___ 1   ___ 2   __
- *           0       0       0       0      Y-1     Y-1     Y-1
- *
- * Node x,y for
- *   0 <= x < X     x = distance along
- *   0 <= y < Y     y = distance across
- *
- * Vertices are in reading order from diagram above ie x varies fastest.
- *
- * Y must be even.  The actual location opposite (0,0) is (X-(Y-1)/2,0),
- * and likewise opposite (0,Y-1) is ((Y-1)/2,0).
- *
- * To label edges, we counte anticlockwise[*] from to-the-right:
- *
- *			    \2   /1
- *			     \  /
- *			  ___ 0   __
- *			  3    1   0
- *			     /  \
- *			   4/   5\
- *
- *   [*] That is, in the direction represented as anticlockwise for
- *       the vertices (0,*)..(4,*) in the diagram above; and of course
- *       that is clockwise for the vertices (X-5,*)..(X-1,*).  The
- *       numbering has an actual discontinuity between (X-1,*) and (0,*).
- *
- * When we iterate over edges, we iterate first over vertices and then
- * over edges 0 to 2, disregarding edges 3 to 5.
  */
 
-#define XBITS 4
-#define X (1<<XBITS)
-#define YBITS 4
-#define Y (1<<YBITS)
-
-/* vertex number:   0000 | y     | x
- *                        YBITS   XBITS
- */
-
-#define N (X*Y)
-#define XMASK (X-1)
-#define YSHIFT XBITS
-#define Y1 (1 << YSHIFT)
-#define YMASK (Y-1 << YSHIFT)
-
   /* Energy has the following parts right now:
    */
   /*
@@ -100,9 +37,10 @@
    *     vd, edge PQ      vd      3
    *                             l
    *
-   *  By symmetry, this calculation gives the same answer
-   *  with R and S exchanged.  Looking at the projection in
-   *  the RMS plane:
+   *  (The dimensions of this are those of F_vd.)
+   *
+   *  By symmetry, this calculation gives the same answer with R and S
+   *  exchanged.  Looking at the projection in the RMS plane:
    *
    *
    *	       	       	       S'
@@ -125,66 +63,10 @@
    *
    *  In practice to avoid division by zero we'll add epsilon to l^3
    *  and the huge energy ought then to be sufficient for the model to
-   *  avoid being close to R=S.  */
-
-#define V6 6
-
-static int edge_end2(unsigned v1, int e) {
-  /* The topology is equivalent to that of a square lattice with only
-   * half of the diagonals.  Ie, the result of shearing the triangular
-   * lattice to make the lines of constant x vertical.  This gives
-   * these six directions:
-   *
-   *                2  1
-   *                | /
-   *                |/
-   *             3--*--0
-   * 		   /|
-   * 		  / |
-   * 		 4  5
-   *
-   * This also handily makes vertical the numbering discontinuity,
-   * where the join happens.
+   *  avoid being close to R=S.
    */
-  static const unsigned dx[V6]= {  +1,  +1,   0,  -1,  -1,   0  },
-                        dy[V6]= {   0, +Y1, +Y1,   0, -Y1, -Y1  };
-  unsigned x, y;
 
-  y= (v1 & YMASK) + dy[e];
-  if (y & ~YMASK) return -1;
-
-  x= (v1 & XMASK) + dx[e];
-  if (x & ~XMASK) {
-    y= (Y-1)*Y1 - y;
-    x &= XMASK;;
-  }
-  return x | y;
-}
-
-#define D3 3
-
-#define FOR_VERTEX(v) \
-  for ((v)=0; (v)<N; (v)++)
-
-#define FOR_VPEDGE(v,e) \
-  for ((e)=0; (e)<V6; (e)++)
-
-#define EDGE_END2 edge_end2
-
-#define FOR_VEDGE(v1,e,v2) \
-  FOR_VPEDGE((v1),(e))
-    if (((v2)= EDGE_END2((v1),(e))) < 0) ; else
-
-#define FOR_EDGE(v1,e,v2) \
-  FOR_VERTEX((v1)) \
-    FOR_VEDGE((v1),(e),(v2))
-
-#define FOR_COORD(k) \
-  for ((k)=0; (k)<D3; (k)++)
-
-#define K FOR_COORD(k)
-
-static double hypotD(const double p[], const double q[]) {
+double hypotD(const double p[D3], const double q[D3]) {
   int k;
   double pq[D3];
   gsl_vector v;
@@ -198,10 +80,16 @@ static double hypotD(const double p[], const double q[]) {
   return gsl_blas_snrm2(&v);
 }
 
+double hypotD2(const double p[D3], const double q[D3]) {
+  double d2= 0;
+  K d2= ffsqa(p[k] - q[k], d2);
+  return d2;
+}
+
 #ifdef FP_FAST_FMA
 # define ffsqa(factor,term) fma((factor),(factor),(term))
 #else
-# define ffsqu(factor,term) ((factor)*(factor)+(term))
+# define ffsqa(factor,term) ((factor)*(factor)+(term))
 #endif
 
 static const l3_epsison= 1e-6;
@@ -219,7 +107,7 @@ static double energy_function(const double vertices[N][D3]) {
     K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
     K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
     b= hypotD(vertices[pi], vertices[qi]);
-    d2= 0; K d2= ffsqa(m[k] - mprime[k], d2);
+    d2= hypotD2(m, mprime);
     l= hypotD(vertices[ri][k] - vertices[si][k]);
     l3 = l*l*l + l3_epsilon;