+ /* Find a linear representation of r using eqs[0]. */
+ rt[0] = -(eqs[0][0]*xt[0] + eqs[0][1]*yt[0])/eqs[0][2];
+ rt[1] = (eqs[0][3] - eqs[0][0]*xt[1] -
+ eqs[0][1]*yt[1])/eqs[0][2];
+
+ /* Construct the quadratic equation. */
+ q[0] = -rt[0]*rt[0];
+ q[1] = -2*rt[0]*rt[1];
+ q[2] = -rt[1]*rt[1];
+ q[0] += xt[0]*xt[0];
+ q[1] += 2*xt[0]*(xt[1]-dots[0]->x);
+ q[2] += (xt[1]-dots[0]->x)*(xt[1]-dots[0]->x);
+ q[0] += yt[0]*yt[0];
+ q[1] += 2*yt[0]*(yt[1]-dots[0]->y);
+ q[2] += (yt[1]-dots[0]->y)*(yt[1]-dots[0]->y);
+
+ /* And solve it. */
+ disc = q[1]*q[1] - 4*q[0]*q[2];
+ if (disc >= 0) {
+ double t;
+
+ disc = sqrt(disc);
+
+ t = (-q[1] + disc) / (2*q[0]);
+ cx[cn] = xt[0]*t + xt[1];
+ cy[cn] = yt[0]*t + yt[1];
+ cn++;
+
+ t = (-q[1] - disc) / (2*q[0]);
+ cx[cn] = xt[0]*t + xt[1];
+ cy[cn] = yt[0]*t + yt[1];
+ cn++;
+ }
+ } else if (nedges == 1) {
+ /*
+ * Two dots and an edge. This one's another
+ * quadratic equation.
+ *
+ * The point we want must lie on the perpendicular
+ * bisector of the two dots; that much is obvious.
+ * So we can construct a parametrisation of that
+ * bisecting line, giving linear formulae for x,y
+ * in terms of t. We can also express the distance
+ * from the edge as such a linear formula.
+ *
+ * Then we set that equal to the radius of the
+ * circle passing through the two points, which is
+ * a Pythagoras exercise; that gives rise to a
+ * quadratic in t, which we solve.
+ */
+ double xt[2], yt[2], rt[2]; /* a,b: {x,y,r}=at+b */
+ double q[3]; /* a,b,c: at^2+bt+c=0 */
+ double disc;
+ double halfsep;
+
+ /* Find parametric formulae for x,y. */
+ {
+ int x1 = dots[0]->x, x2 = dots[1]->x;
+ int y1 = dots[0]->y, y2 = dots[1]->y;
+ int dx = x2-x1, dy = y2-y1;
+ double d = sqrt((double)dx*dx + (double)dy*dy);
+
+ xt[1] = (x1+x2)/2.0;
+ yt[1] = (y1+y2)/2.0;
+ /* It's convenient if we have t at standard scale. */
+ xt[0] = -dy/d;
+ yt[0] = dx/d;
+
+ /* Also note down half the separation between
+ * the dots, for use in computing the circle radius. */
+ halfsep = 0.5*d;
+ }
+
+ /* Find a parametric formula for r. */
+ {
+ int x1 = edgedot1[0]->x, x2 = edgedot2[0]->x;
+ int y1 = edgedot1[0]->y, y2 = edgedot2[0]->y;
+ int dx = x2-x1, dy = y2-y1;
+ double d = sqrt((double)dx*dx + (double)dy*dy);
+ rt[0] = (xt[0]*dy - yt[0]*dx) / d;
+ rt[1] = ((xt[1]-x1)*dy - (yt[1]-y1)*dx) / d;
+ }
+
+ /* Construct the quadratic equation. */
+ q[0] = rt[0]*rt[0];
+ q[1] = 2*rt[0]*rt[1];
+ q[2] = rt[1]*rt[1];
+ q[0] -= 1;
+ q[2] -= halfsep*halfsep;
+
+ /* And solve it. */
+ disc = q[1]*q[1] - 4*q[0]*q[2];
+ if (disc >= 0) {
+ double t;
+
+ disc = sqrt(disc);
+
+ t = (-q[1] + disc) / (2*q[0]);
+ cx[cn] = xt[0]*t + xt[1];
+ cy[cn] = yt[0]*t + yt[1];
+ cn++;
+
+ t = (-q[1] - disc) / (2*q[0]);
+ cx[cn] = xt[0]*t + xt[1];
+ cy[cn] = yt[0]*t + yt[1];
+ cn++;
+ }
+ } else if (nedges == 0) {
+ /*
+ * Three dots. This is another linear matrix
+ * equation, this time with each row of the matrix
+ * representing the perpendicular bisector between
+ * two of the points. Of course we only need two
+ * such lines to find their intersection, so we
+ * need only solve a 2x2 matrix equation.
+ */
+
+ double matrix[4], vector[2], vector2[2];
+ int m;
+
+ for (m = 0; m < 2; m++) {
+ int x1 = dots[m]->x, x2 = dots[m+1]->x;
+ int y1 = dots[m]->y, y2 = dots[m+1]->y;
+ int dx = x2-x1, dy = y2-y1;
+
+ /*
+ * ((x,y) - (x1,y1)) . (dx,dy) = 1/2 |(dx,dy)|^2
+ *
+ * => 2x dx + 2y dy = dx^2+dy^2 + (2 x1 dx + 2 y1 dy)
+ */
+ matrix[2*m+0] = 2*dx;
+ matrix[2*m+1] = 2*dy;
+ vector[m] = ((double)dx*dx + (double)dy*dy +
+ 2.0*x1*dx + 2.0*y1*dy);
+ }
+
+ if (solve_2x2_matrix(matrix, vector, vector2)) {
+ cx[cn] = vector2[0];
+ cy[cn] = vector2[1];
+ cn++;
+ }