/**************************************************************************** PURPOSE : Try to find circles going through N lattice points. AUTHOR : Th. Wolf LANGUAGE: ANSI-C --------------------------------------------------------------------------- HISTORY : 23-May-1996 TW Created. 24-May-1996 TW Added SE triangle of NE quadrant. ****************************************************************************/ #include #include #include #include #include /****************************************************************************/ #define TRUE 1 #define FALSE 0 typedef struct POINT { double x, y; } POINT; typedef struct CIRCLE { POINT c; double r; int flag; /* Will be TRUE if we suspect that the solution we found might not be the smallest one. */ int *X, *Y; /* Dynamically allocated arrays holding the lattice points on the circle's perimeter. */ } CIRCLE; /****************************************************************************/ #define N_MAX 500 CIRCLE circles[N_MAX]; int suspect; /****************************************************************************/ int eq (double x, double y) { return (fabs (x - y) < DBL_EPSILON); /* Good enough for our purposes... */ } /****************************************************************************/ void get_circle (POINT p[], POINT *center, double *rr) /**** Note: 'rr' returns the square of the radius. */ { POINT sq[3]; int i; double den, nom; #define X(i) (p[i].x) #define Y(i) (p[i].y) #define SX(i) (sq[i].x) #define SY(i) (sq[i].y) for (i = 0; i < 3; i++) { SX(i) = X(i) * X(i); SY(i) = Y(i) * Y(i); } /* end for */ #if 0 /**** Note: these are the full formulas for calculating the center of a circle given three points. However, since p[0] always is the origin, we can simplify somewhat. */ den = Y(2) * (SX(0) - SX(1) + SY(0) - SY(1)) + Y(0) * (SX(1) - SX(2) + SY(1) - SY(2)) + Y(1) * (SX(2) - SX(0) + SY(2) - SY(0)); nom = 2 * ((X(0) - X(2)) * (Y(0) - Y(1)) - (X(1) - X(0)) * (Y(2) - Y(0))); center->x = den / nom; if (Y(0) == Y(1)) { den = 2 * center->x * (X(0) - X(2)) + SX(2) + SY(2) - SX(0) - SY(0); nom = 2 * (Y(2) - Y(0)); } else { den = 2 * center->x * (X(1) - X(0)) + SX(0) + SY(0) - SX(1) - SY(1); nom = 2 * (Y(0) - Y(1)); } /* end if */ center->y = den / nom; #else /**** These are the simplified formulas. Interestingly, there's no per- ceptible speed increase. */ den = Y(1) * (SX(2) + SY(2)) - Y(2) * (SX(1) + SY(1)); nom = 2 * (X(2) * Y(1) - X(1) * Y(2)); center->x = den / nom; if (Y(1) == 0.0) { den = 2 * center->x * X(2) + SX(2) + SY(2); nom = 2 * Y(2); } else { den = 2 * center->x * X(1) - SX(1) - SY(1); nom = -2 * Y(1); } /* end if */ center->y = den / nom; #endif *rr = center->x * center->x + center->y * center->y; /**** Now move the center as near as possible to the origin. */ center->x -= floor (center->x); center->y -= floor (center->y); #undef X #undef Y #undef SX #undef SY } /* end get_circle */ /****************************************************************************/ void check_points (POINT c, double rr) /**** For each integral x, compute the y coordinates and check whether they fall on a lattice point (i.e. the y coordinate also is integral). */ { double ry, /* The real y coordinate */ r = sqrt (rr); int left = floor (c.x - r); int right = ceil (c.x + r); int n, x, y; int X[N_MAX], Y[N_MAX]; /* Stores the coordinates of lattice points hit. */ n = 0; /* Number of lattice points hit so far */ for (x = left; x <= right; x++) { ry = rr - (x - c.x) * (x - c.x); if (eq (ry, 0.0)) ry = 0.0; /* Catch rounding errors */ if (ry >= 0) { ry = sqrt (ry) + c.y; y = ry + 0.5; /* Round to nearest integer */ if (eq (ry, y)) { if (n < N_MAX) { X[n] = x; Y[n] = y; } /* end if */ n++; } /* end if */ if (!eq (ry, c.y)) { /* Avoid counting extremal points twice */ ry = 2 * c.y - ry; y = ry - 0.5; if (eq (ry, y)) { if (n < N_MAX) { X[n] = x; Y[n] = y; } /* end if */ n++; } /* end if */ } /* end if */ } /* end if */ if (n >= N_MAX) return; /* We can't store it, so quit. */ } /* end for */ /**** Store the circle and its lattice points. */ if (r < circles[n].r || circles[n].r == 0.0) { if (circles[n].r == 0.0) { circles[n].X = malloc (n * sizeof (int)); circles[n].Y = malloc (n * sizeof (int)); circles[n].flag = suspect; } circles[n].r = r; circles[n].c = c; memcpy (circles[n].X, X, n * sizeof (int)); memcpy (circles[n].Y, Y, n * sizeof (int)); } /* end if */ } /* end check_points */ /****************************************************************************/ void generate_circles (int size) { /**** Idea: keep one point at the origin, the rotate two vectors around it. This gives us the three points we need for a circle. */ POINT a, b, c, p[3]; double r; p[0].x = p[0].y = 0.0; /**** First, try all combinations with both vectors in the lower half of the square. */ suspect = FALSE; for (a.y = 0; a.y >= -size; a.y--) { for (a.x = -size; a.x <= size; a.x++) { if (a.x == 0.0 && a.y == 0.0) continue; for (b.y = a.y; b.y >= -size; b.y--) { for (b.x = -size; b.x <= size; b.x++) { if (b.x == 0.0 && b.y == 0.0) continue; if (a.x == 0.0 || b.x == 0.0) { if (a.x == 0.0 && b.x == 0.0) continue; } else if (eq (a.y / a.x, b.y / b.x)) continue; /**** The three points are not all on one line. */ p[1] = a; p[2] = b; get_circle (p, &c, &r); check_points (c, r); } /* end for */ } /* end for */ } /* end for */ } /* end for */ /**** Now try all combinations with one vector in the lower half of the square and the other vector in the SE triangle of the NE quadrant. Note: this will sometimes find more solutions, but they're not necessarily minimal. They will be flagged in the output as such. */ suspect = TRUE; for (a.y = 0; a.y >= -size; a.y--) { for (a.x = -size; a.x <= size; a.x++) { if (a.x == 0.0 && a.y == 0.0) continue; for (b.y = 0; b.y <= size; b.y++) { for (b.x = b.y; b.x <= size; b.x++) { if (b.x == 0.0 && b.y == 0.0) continue; if (a.x == 0.0 || b.x == 0.0) { if (a.x == 0.0 && b.x == 0.0) continue; } else if (eq (a.y / a.x, b.y / b.x)) continue; /**** The three points are not all on one line. */ p[1] = a; p[2] = b; get_circle (p, &c, &r); check_points (c, r); } /* end for */ } /* end for */ } /* end for */ } /* end for */ } /* end generate_circles */ /****************************************************************************/ void print_points (int *x, int *y, int n) { int i = 0; while (i < n) { printf ("(%5d,%5d) ", *x++, *y++); if (++i % 5 == 0) putc ('\n', stdout); } if (n % 5 != 0) putc( '\n', stdout); } /* end print_points */ /****************************************************************************/ int main (int argc, char *argv[]) { int i, size, n; char *prog_name; clock_t start; if (argc > 1) { prog_name = (argv[0] ? argv[0] : "circle_lattice"); n = sscanf (argv[1], "%d", &size); if (!n) { fprintf (stderr, "%s: argument must be a number!\n", prog_name); return EXIT_FAILURE; } else if (size <= 0) { fprintf (stderr, "%s: argument must be positive!\n", prog_name); return EXIT_FAILURE; } else if (size > 30) { fprintf (stderr, "%s: argument too large!\n", prog_name); return EXIT_FAILURE; } else if (argc > 2) { fprintf (stderr, "%s: only one (integer) argument is allowed!\n", prog_name); fprintf (stderr, "%s: other arguments are ignored!\n", prog_name); } /* end if */ } else { size = 6; } /* end if */ start = clock (); generate_circles (size); start = clock () - start; printf ("The computation took %f seconds.\n", (double) start / CLOCKS_PER_SEC); for (i = 0; i < N_MAX; i++) { if (circles[i].r) { printf ("\nN=%3d: c = (%f, %f), r = %f%s\n", i, circles[i].c.x, circles[i].c.y, circles[i].r, circles[i].flag ? " ***" : ""); print_points (circles[i].X, circles[i].Y, i); free (circles[i].X); free (circles[i].Y); } /* end if */ } /* end for */ return EXIT_SUCCESS; } /* end main */ /****************************************************************************/ /* end circle_lattice.c */