1 /* -*- Mode: C; tab-width: 4 -*- */
2 /* apollonian --- Apollonian Circles */
3
4 #if 0
5 static const char sccsid[] = "@(#)apollonian.c 5.02 2001/07/01 xlockmore";
6 #endif
7
8 /*-
9 * Copyright (c) 2000, 2001 by Allan R. Wilks <allan@research.att.com>.
10 *
11 * Permission to use, copy, modify, and distribute this software and its
12 * documentation for any purpose and without fee is hereby granted,
13 * provided that the above copyright notice appear in all copies and that
14 * both that copyright notice and this permission notice appear in
15 * supporting documentation.
16 *
17 * This file is provided AS IS with no warranties of any kind. The author
18 * shall have no liability with respect to the infringement of copyrights,
19 * trade secrets or any patents by this file or any part thereof. In no
20 * event will the author be liable for any lost revenue or profits or
21 * other special, indirect and consequential damages.
22 *
23 * radius r = 1 / c (curvature)
24 *
25 * Descartes Circle Theorem: (a, b, c, d are curvatures of tangential circles)
26 * Let a, b, c, d be the curvatures of for mutually (externally) tangent
27 * circles in the plane. Then
28 * a^2 + b^2 + c^2 + d^2 = (a + b + c + d)^2 / 2
29 *
30 * Complex Descartes Theorem: If the oriented curvatues and (complex) centers
31 * of an oriented Descrates configuration in the plane are a, b, c, d and
32 * w, x, y, z respectively, then
33 * a^2*w^2 + b^2*x^2 + c^2*y^2 + d^2*z^2 = (aw + bx + cy + dz)^2 / 2
34 * In addition these quantities satisfy
35 * a^2*w + b^2*x + c^2*y + d^2*z = (aw + bx + cy + dz)(a + b + c + d) / 2
36 *
37 * Enumerate root integer Descartes quadruples (a,b,c,d) satisfying the
38 * Descartes condition:
39 * 2(a^2+b^2+c^2+d^2) = (a+b+c+d)^2
40 * i.e., quadruples for which no application of the "pollinate" operator
41 * z <- 2(a+b+c+d) - 3*z,
42 * where z is in {a,b,c,d}, gives a quad of strictly smaller sum. This
43 * is equivalent to the condition:
44 * sum(a,b,c,d) >= 2*max(a,b,c,d)
45 * which, because of the Descartes condition, is equivalent to
46 * sum(a^2,b^2,c^2,d^2) >= 2*max(a,b,c,d)^2
47 *
48 *
49 * Todo:
50 * Add a small font
51 *
52 * Revision History:
53 * 25-Jun-2001: Converted from C and Postscript code by David Bagley
54 * Original code by Allan R. Wilks <allan@research.att.com>.
55 *
56 * From Circle Math Science News April 21, 2001 VOL. 254-255
57 * http://www.sciencenews.org/20010421/toc.asp
58 * Apollonian Circle Packings Assorted papers from Ronald L Graham,
59 * Jeffrey Lagarias, Colin Mallows, Allan Wilks, Catherine Yan
60 * http://front.math.ucdavis.edu/math.NT/0009113
61 * http://front.math.ucdavis.edu/math.MG/0101066
62 * http://front.math.ucdavis.edu/math.MG/0010298
63 * http://front.math.ucdavis.edu/math.MG/0010302
64 * http://front.math.ucdavis.edu/math.MG/0010324
65 */
66
67 #ifdef STANDALONE
68 # define MODE_apollonian
69 # define DEFAULTS "*delay: 1000000 \n" \
70 "*count: 64 \n" \
71 "*cycles: 20 \n" \
72 "*ncolors: 64 \n" \
73 "*font: fixed" "\n" \
74 "*fpsTop: true \n" \
75 "*fpsSolid: true \n" \
76 "*ignoreRotation: True" \
77
78 # define release_apollonian 0
79 # define reshape_apollonian 0
80 # define apollonian_handle_event 0
81 # include "xlockmore.h" /* in xscreensaver distribution */
82 #else /* STANDALONE */
83 # include "xlock.h" /* in xlockmore distribution */
84 #endif /* STANDALONE */
85
86 #ifdef MODE_apollonian
87
88 #define DEF_ALTGEOM "True"
89 #define DEF_LABEL "True"
90
91 static Bool altgeom;
92 static Bool label;
93
94 static XrmOptionDescRec opts[] =
95 {
96 {"-altgeom", ".apollonian.altgeom", XrmoptionNoArg, "on"},
97 {"+altgeom", ".apollonian.altgeom", XrmoptionNoArg, "off"},
98 {"-label", ".apollonian.label", XrmoptionNoArg, "on"},
99 {"+label", ".apollonian.label", XrmoptionNoArg, "off"},
100 };
101 static argtype vars[] =
102 {
103 {&altgeom, "altgeom", "AltGeom", DEF_ALTGEOM, t_Bool},
104 {&label, "label", "Label", DEF_LABEL, t_Bool},
105 };
106 static OptionStruct desc[] =
107 {
108 {"-/+altgeom", "turn on/off alternate geometries (off euclidean space, on includes spherical and hyperbolic)"},
109 {"-/+label", "turn on/off alternate space and number labeling"},
110 };
111
112 ENTRYPOINT ModeSpecOpt apollonian_opts =
113 {sizeof opts / sizeof opts[0], opts, sizeof vars / sizeof vars[0], vars, desc};
114
115 #ifdef DOFONT
116 extern XFontStruct *getFont(Display * display);
117 #endif
118
119 #ifdef USE_MODULES
120 ModStruct apollonian_description =
121 {"apollonian", "init_apollonian", "draw_apollonian", (char *) NULL,
122 "init_apollonian", "init_apollonian", "free_apollonian", &apollonian_opts,
123 1000000, 64, 20, 1, 64, 1.0, "",
124 "Shows Apollonian Circles", 0, NULL};
125
126 #endif
127
128 typedef struct {
129 int a, b, c, d;
130 } apollonian_quadruple;
131
132 typedef struct {
133 double e; /* euclidean bend */
134 double s; /* spherical bend */
135 double h; /* hyperbolic bend */
136 double x, y; /* euclidean bend times euclidean position */
137 } circle;
138 enum space {
139 euclidean = 0, spherical, hyperbolic
140 };
141
142 static const char * space_string[] = {
143 "euclidean",
144 "spherical",
145 "hyperbolic"
146 };
147
148 /*
149 Generate Apollonian packing starting with a quadruple of circles.
150 The four input lines each contain the 5-tuple (e,s,h,x,y) representing
151 the circle with radius 1/e and center (x/e,y/e). The s and h is propagated
152 like e, x and y, but can differ from e so as to represent different
153 geometries, spherical and hyperbolic, respectively. The "standard" picture,
154 for example (-1, 2, 2, 3), can be labeled for the three geometries.
155 Origins of circles z1, z2, z3, z4
156 a * z1 = 0
157 b * z2 = (a+b)/a
158 c * z3 = (q123 + a * i)^2/(a*(a+b)) where q123 = sqrt(a*b+a*c+b*c)
159 d * z4 = (q124 + a * i)^2/(a*(a+b)) where q124 = q123 - a - b
160 If (e,x,y) represents the Euclidean circle (1/e,x/e,y/e) (so that e is
161 the label in the standard picture) then the "spherical label" is
162 (e^2+x^2+y^2-1)/(2*e) (an integer!) and the "hyperbolic label", is
163 calulated by h + s = e.
164 */
165 static circle examples[][4] = {
166 { /* double semi-bounded */
167 { 0, 0, 0, 0, 1},
168 { 0, 0, 0, 0, -1},
169 { 1, 1, 1, -1, 0},
170 { 1, 1, 1, 1, 0}
171 },
172 #if 0
173 { /* standard */
174 {-1, 0, -1, 0, 0},
175 { 2, 1, 1, 1, 0},
176 { 2, 1, 1, -1, 0},
177 { 3, 2, 1, 0, 2}
178 },
179 { /* next simplest */
180 {-2, -1, -1, 0.0, 0},
181 { 3, 2, 1, 0.5, 0},
182 { 6, 3, 3, -2.0, 0},
183 { 7, 4, 3, -1.5, 2}
184 },
185 { /* */
186 {-3, -2, -1, 0.0, 0},
187 { 4, 3, 1, 1.0 / 3.0, 0},
188 {12, 7, 5, -3.0, 0},
189 {13, 8, 5, -8.0 / 3.0, 2}
190 },
191 { /* Mickey */
192 {-3, -2, -1, 0.0, 0},
193 { 5, 4, 1, 2.0 / 3.0, 0},
194 { 8, 5, 3, -4.0 / 3.0, -1},
195 { 8, 5, 3, -4.0 / 3.0, 1}
196 },
197 { /* */
198 {-4, -3, -1, 0.00, 0},
199 { 5, 4, 1, 0.25, 0},
200 {20, 13, 7, -4.00, 0},
201 {21, 14, 7, -3.75, 2}
202 },
203 { /* Mickey2 */
204 {-4, -2, -2, 0.0, 0},
205 { 8, 4, 4, 1.0, 0},
206 { 9, 5, 4, -0.75, -1},
207 { 9, 5, 4, -0.75, 1}
208 },
209 { /* Mickey3 */
210 {-5, -4, -1, 0.0, 0},
211 { 7, 6, 1, 0.4, 0},
212 {18, 13, 5, -2.4, -1},
213 {18, 13, 5, -2.4, 1}
214 },
215 { /* */
216 {-6, -5, -1, 0.0, 0},
217 { 7, 6, 1, 1.0 / 6.0, 0},
218 {42, 31, 11, -6.0, 0},
219 {43, 32, 11, -35.0 / 6.0, 2}
220 },
221 { /* */
222 {-6, -3, -3, 0.0, 0},
223 {10, 5, 5, 2.0 / 3.0, 0},
224 {15, 8, 7, -1.5, 0},
225 {19, 10, 9, -5.0 / 6.0, 2}
226 },
227 { /* asymmetric */
228 {-6, -5, -1, 0.0, 0.0},
229 {11, 10, 1, 5.0 / 6.0, 0.0},
230 {14, 11, 3, -16.0 / 15.0, -0.8},
231 {15, 12, 3, -0.9, 1.2}
232 },
233 #endif
234 /* Non integer stuff */
235 #define DELTA 2.154700538 /* ((3+2*sqrt(3))/3) */
236 { /* 3 fold symmetric bounded (x, y calculated later) */
237 { -1, -1, -1, 0.0, 0.0},
238 {DELTA, DELTA, DELTA, 1.0, 0.0},
239 {DELTA, DELTA, DELTA, 1.0, -1.0},
240 {DELTA, DELTA, DELTA, -1.0, 1.0}
241 },
242 { /* semi-bounded (x, y calculated later) */
243 #define ALPHA 2.618033989 /* ((3+sqrt(5))/2) */
244 { 1.0, 1.0, 1.0, 0, 0},
245 { 0.0, 0.0, 0.0, 0, -1},
246 {1.0/(ALPHA*ALPHA), 1.0/(ALPHA*ALPHA), 1.0/(ALPHA*ALPHA), -1, 0},
247 { 1.0/ALPHA, 1.0/ALPHA, 1.0/ALPHA, -1, 0}
248 },
249 { /* unbounded (x, y calculated later) */
250 /* #define PHI 1.618033989 *//* ((1+sqrt(5))/2) */
251 #define BETA 2.890053638 /* (PHI+sqrt(PHI)) */
252 { 1.0, 1.0, 1.0, 0, 0},
253 {1.0/(BETA*BETA*BETA), 1.0/(BETA*BETA*BETA), 1.0/(BETA*BETA*BETA), 1, 0},
254 { 1.0/(BETA*BETA), 1.0/(BETA*BETA), 1.0/(BETA*BETA), 1, 0},
255 { 1.0/BETA, 1.0/BETA, 1.0/BETA, 1, 0}
256 }
257 };
258
259 #define PREDEF_CIRCLE_GAMES (sizeof (examples) / (4 * sizeof (circle)))
260
261 #if 0
262 Euclidean
263 0, 0, 1, 1
264 -1, 2, 2, 3
265 -2, 3, 6, 7
266 -3, 5, 8, 8
267 -4, 8, 9, 9
268 -3, 4, 12, 13
269 -6, 11, 14, 15
270 -5, 7, 18, 18
271 -6, 10, 15, 19
272 -7, 12, 17, 20
273 -4, 5, 20, 21
274 -9, 18, 19, 22
275 -8, 13, 21, 24
276 Spherical
277 0, 1, 1, 2
278 -1, 2, 3, 4
279 -2, 4, 5, 5
280 -2, 3, 7, 8
281 Hyperbolic
282 -1, 1, 1, 1
283 0, 0, 1, 3
284 -2, 3, 5, 6
285 -3, 6, 6, 7
286 #endif
287
288 typedef struct {
289 int size;
290 XPoint offset;
291 int geometry;
292 circle c1, c2, c3, c4;
293 int color_offset;
294 int count;
295 Bool label, altgeom;
296 apollonian_quadruple *quad;
297 #ifdef DOFONT
298 XFontStruct *font;
299 #endif
300 int time;
301 int game;
302 } apollonianstruct;
303
304 static apollonianstruct *apollonians = (apollonianstruct *) NULL;
305
306 #define FONT_HEIGHT 19
307 #define FONT_WIDTH 15
308 #define FONT_LENGTH 20
309 #define MAX_CHAR 10
310 #define K 2.15470053837925152902 /* 1+2/sqrt(3) */
311 #define MAXBEND 100 /* Do not want configurable by user since it will take too
312 much time if increased. */
313
314 static int
gcd(int a,int b)315 gcd(int a, int b)
316 {
317 int r;
318
319 while (b) {
320 r = a % b;
321 a = b;
322 b = r;
323 }
324 return a;
325 }
326
327 static int
isqrt(int n)328 isqrt(int n)
329 {
330 int y;
331
332 if (n < 0)
333 return -1;
334 y = (int) (sqrt((double) n) + 0.5);
335 return ((n == y*y) ? y : -1);
336 }
337
338 static void
dquad(int n,apollonian_quadruple * quad)339 dquad(int n, apollonian_quadruple *quad)
340 {
341 int a, b, c, d;
342 int counter = 0, B, C;
343
344 for (a = 0; a < MAXBEND; a++) {
345 B = (int) (K * a);
346 for (b = a + 1; b <= B; b++) {
347 C = (int) (((a + b) * (a + b)) / (4.0 * (b - a)));
348 for (c = b; c <= C; c++) {
349 d = isqrt(b*c-a*(b+c));
350 if (d >= 0 && (gcd(a,gcd(b,c)) <= 1)) {
351 quad[counter].a = -a;
352 quad[counter].b = b;
353 quad[counter].c = c;
354 quad[counter].d = -a+b+c-2*d;
355 if (++counter >= n) {
356 return;
357 }
358 }
359 }
360 }
361 }
362 (void) printf("found only %d below maximum bend of %d\n",
363 counter, MAXBEND);
364 for (; counter < n; counter++) {
365 quad[counter].a = -1;
366 quad[counter].b = 2;
367 quad[counter].c = 2;
368 quad[counter].d = 3;
369 }
370 return;
371 }
372
373 /*
374 * Given a Descartes quadruple of bends (a,b,c,d), with a<0, find a
375 * quadruple of circles, represented by (bend,bend*x,bend*y), such
376 * that the circles have the given bends and the bends times the
377 * centers are integers.
378 *
379 * This just performs an exaustive search, assuming that the outer
380 * circle has center in the unit square.
381 *
382 * It is always sufficient to look in {(x,y):0<=y<=x<=1/2} for the
383 * center of the outer circle, but this may not lead to a packing
384 * that can be labelled with integer spherical and hyperbolic labels.
385 * To effect the smaller search, replace FOR(a) with
386 *
387 * for (pa = ea/2; pa <= 0; pa++) for (qa = pa; qa <= 0; qa++)
388 */
389
390 #define For(v,l,h) for (v = l; v <= h; v++)
391 #define FOR(z) For(p##z,lop##z,hip##z) For(q##z,loq##z,hiq##z)
392 #define H(z) ((e##z*e##z+p##z*p##z+q##z*q##z)%2)
393 #define UNIT(z) ((abs(e##z)-1)*(abs(e##z)-1) >= p##z*p##z+q##z*q##z)
394 #define T(z,w) is_tangent(e##z,p##z,q##z,e##w,p##w,q##w)
395 #define LO(r,z) lo##r##z = iceil(e##z*(r##a+1),ea)-1
396 #define HI(r,z) hi##r##z = iflor(e##z*(r##a-1),ea)-1
397 #define B(z) LO(p,z); HI(p,z); LO(q,z); HI(q,z)
398
399 static int
is_quad(int a,int b,int c,int d)400 is_quad(int a, int b, int c, int d)
401 {
402 int s;
403
404 s = a+b+c+d;
405 return 2*(a*a+b*b+c*c+d*d) == s*s;
406 }
407
408 static Bool
is_tangent(int e1,int p1,int q1,int e2,int p2,int q2)409 is_tangent(int e1, int p1, int q1, int e2, int p2, int q2)
410 {
411 int dx, dy, s;
412
413 dx = p1*e2 - p2*e1;
414 dy = q1*e2 - q2*e1;
415 s = e1 + e2;
416 return dx*dx + dy*dy == s*s;
417 }
418
419 static int
iflor(int a,int b)420 iflor(int a, int b)
421 {
422 int q;
423
424 if (b == 0) {
425 (void) printf("iflor: b = 0\n");
426 return 0;
427 }
428 if (a%b == 0)
429 return a/b;
430 q = abs(a)/abs(b);
431 return ((a<0)^(b<0)) ? -q-1 : q;
432 }
433
434 static int
iceil(int a,int b)435 iceil(int a, int b)
436 {
437 int q;
438
439 if (b == 0) {
440 (void) printf("iceil: b = 0\n");
441 return 0;
442 }
443 if (a%b == 0)
444 return a/b;
445 q = abs(a)/abs(b);
446 return ((a<0)^(b<0)) ? -q : 1+q;
447 }
448
449 static double
geom(int geometry,int e,int p,int q)450 geom(int geometry, int e, int p, int q)
451 {
452 int g = (geometry == spherical) ? -1 :
453 (geometry == hyperbolic) ? 1 : 0;
454
455 if (g)
456 return (e*e + (1.0 - p*p - q*q) * g) / (2.0*e);
457 (void) printf("geom: g = 0\n");
458 return e;
459 }
460
461 static void
cquad(circle * c1,circle * c2,circle * c3,circle * c4)462 cquad(circle *c1, circle *c2, circle *c3, circle *c4)
463 {
464 int ea, eb, ec, ed;
465 int pa, pb, pc, pd;
466 int qa, qb, qc, qd;
467 int lopa, lopb, lopc, lopd;
468 int hipa, hipb, hipc, hipd;
469 int loqa, loqb, loqc, loqd;
470 int hiqa, hiqb, hiqc, hiqd;
471
472 ea = (int) c1->e;
473 eb = (int) c2->e;
474 ec = (int) c3->e;
475 ed = (int) c4->e;
476 if (ea >= 0)
477 (void) printf("ea = %d\n", ea);
478 if (!is_quad(ea,eb,ec,ed))
479 (void) printf("Error not quad %d %d %d %d\n", ea, eb, ec, ed);
480 lopa = loqa = ea;
481 hipa = hiqa = 0;
482 FOR(a) {
483 B(b); B(c); B(d);
484 if (H(a) && UNIT(a)) FOR(b) {
485 if (H(b) && T(a,b)) FOR(c) {
486 if (H(c) && T(a,c) && T(b,c)) FOR(d) {
487 if (H(d) && T(a,d) && T(b,d) && T(c,d)) {
488 c1->s = geom(spherical, ea, pa, qa);
489 c1->h = geom(hyperbolic, ea, pa, qa);
490 c2->s = geom(spherical, eb, pb, qb);
491 c2->h = geom(hyperbolic, eb, pb, qb);
492 c3->s = geom(spherical, ec, pc, qc);
493 c3->h = geom(hyperbolic, ec, pc, qc);
494 c4->s = geom(spherical, ed, pd, qd);
495 c4->h = geom(hyperbolic, ed, pd, qd);
496 }
497 }
498 }
499 }
500 }
501 }
502
503 static void
p(ModeInfo * mi,circle c)504 p(ModeInfo *mi, circle c)
505 {
506 apollonianstruct *cp = &apollonians[MI_SCREEN(mi)];
507 char string[15];
508 double g, e;
509 int g_width;
510
511 #ifdef DEBUG
512 (void) printf("c.e=%g c.s=%g c.h=%g c.x=%g c.y=%g\n",
513 c.e, c.s, c.h, c.x, c.y);
514 #endif
515 g = (cp->geometry == spherical) ? c.s : (cp->geometry == hyperbolic) ?
516 c.h : c.e;
517 if (c.e < 0.0) {
518 if (g < 0.0)
519 g = -g;
520 if (MI_NPIXELS(mi) <= 2)
521 XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
522 MI_WHITE_PIXEL(mi));
523 else
524 XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
525 MI_PIXEL(mi, ((int) ((g + cp->color_offset) *
526 g)) % MI_NPIXELS(mi)));
527 XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
528 ((int) (cp->size * (-cp->c1.e) * (c.x - 1.0) /
529 (-2.0 * c.e) + cp->size / 2.0 + cp->offset.x)),
530 ((int) (cp->size * (-cp->c1.e) * (c.y - 1.0) /
531 (-2.0 * c.e) + cp->size / 2.0 + cp->offset.y)),
532 (int) (cp->c1.e * cp->size / c.e),
533 (int) (cp->c1.e * cp->size / c.e), 0, 23040);
534 if (!cp->label) {
535 #ifdef DEBUG
536 (void) printf("%g\n", -g);
537 #endif
538 return;
539 }
540 (void) sprintf(string, "%g", (g == 0.0) ? 0 : -g);
541 if (cp->size >= 10 * FONT_WIDTH) {
542 /* hard code these to corners */
543 XDrawString(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
544 ((int) (cp->size * c.x / (2.0 * c.e))) + cp->offset.x,
545 ((int) (cp->size * c.y / (2.0 * c.e))) + FONT_HEIGHT,
546 string, (g == 0.0) ? 1 : ((g < 10.0) ? 2 :
547 ((g < 100.0) ? 3 : 4)));
548 }
549 if (cp->altgeom && MI_HEIGHT(mi) >= 30 * FONT_WIDTH) {
550 XDrawString(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
551 ((int) (cp->size * c.x / (2.0 * c.e) + cp->offset.x)),
552 ((int) (cp->size * c.y / (2.0 * c.e) + MI_HEIGHT(mi) -
553 FONT_HEIGHT / 2)), (char *) space_string[cp->geometry],
554 strlen(space_string[cp->geometry]));
555 }
556 return;
557 }
558 if (MI_NPIXELS(mi) <= 2)
559 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
560 else
561 XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
562 MI_PIXEL(mi, ((int) ((g + cp->color_offset) * g)) %
563 MI_NPIXELS(mi)));
564 if (c.e == 0.0) {
565 if (c.x == 0.0 && c.y != 0.0) {
566 XDrawLine(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
567 0, (int) ((c.y + 1.0) * cp->size / 2.0 + cp->offset.y),
568 MI_WIDTH(mi),
569 (int) ((c.y + 1.0) * cp->size / 2.0 + cp->offset.y));
570 } else if (c.y == 0.0 && c.x != 0.0) {
571 XDrawLine(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
572 (int) ((c.x + 1.0) * cp->size / 2.0 + cp->offset.x), 0,
573 (int) ((c.x + 1.0) * cp->size / 2.0 + cp->offset.x),
574 MI_HEIGHT(mi));
575 }
576 return;
577 }
578 e = (cp->c1.e >= 0.0) ? 1.0 : -cp->c1.e;
579 XFillArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
580 ((int) (cp->size * e * (c.x - 1.0) / (2.0 * c.e) +
581 cp->size / 2.0 + cp->offset.x)),
582 ((int) (cp->size * e * (c.y - 1.0) / (2.0 * c.e) +
583 cp->size / 2.0 + cp->offset.y)),
584 (int) (e * cp->size / c.e), (int) (e * cp->size / c.e),
585 0, 23040);
586 if (!cp->label) {
587 #ifdef DEBUG
588 (void) printf("%g\n", g);
589 #endif
590 return;
591 }
592 if (MI_NPIXELS(mi) <= 2)
593 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_BLACK_PIXEL(mi));
594 else
595 XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
596 MI_PIXEL(mi, ((int) ((g + cp->color_offset) * g) +
597 MI_NPIXELS(mi) / 2) % MI_NPIXELS(mi)));
598 g_width = (g < 10.0) ? 1: ((g < 100.0) ? 2 : 3);
599 if (c.e < e * cp->size / (FONT_LENGTH + 5 * g_width) && g < 1000.0) {
600 (void) sprintf(string, "%g", g);
601 XDrawString(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
602 ((int) (cp->size * e * c.x / (2.0 * c.e) +
603 cp->size / 2.0 + cp->offset.x)) -
604 g_width * FONT_WIDTH / 2,
605 ((int) (cp->size * e * c.y / (2.0 * c.e) +
606 cp->size / 2.0 + cp->offset.y)) +
607 FONT_HEIGHT / 2,
608 string, g_width);
609 }
610 }
611
612 #define BIG 7
613 static void
f(ModeInfo * mi,circle c1,circle c2,circle c3,circle c4,int depth)614 f(ModeInfo *mi, circle c1, circle c2, circle c3, circle c4, int depth)
615 {
616 apollonianstruct *cp = &apollonians[MI_SCREEN(mi)];
617 int e = (int) ((cp->c1.e >= 0.0) ? 1.0 : -cp->c1.e);
618 circle c;
619
620 if (depth > mi->recursion_depth) mi->recursion_depth = depth;
621
622 c.e = 2*(c1.e+c2.e+c3.e) - c4.e;
623 c.s = 2*(c1.s+c2.s+c3.s) - c4.s;
624 c.h = 2*(c1.h+c2.h+c3.h) - c4.h;
625 c.x = 2*(c1.x+c2.x+c3.x) - c4.x;
626 c.y = 2*(c1.y+c2.y+c3.y) - c4.y;
627 if (c.e == 0 ||
628 c.e > cp->size * e || c.x / c.e > BIG || c.y / c.e > BIG ||
629 c.x / c.e < -BIG || c.y / c.e < -BIG)
630 return;
631 p(mi, c);
632 f(mi, c2, c3, c, c1, depth+1);
633 f(mi, c1, c3, c, c2, depth+1);
634 f(mi, c1, c2, c, c3, depth+1);
635 }
636
637 ENTRYPOINT void
free_apollonian(ModeInfo * mi)638 free_apollonian (ModeInfo * mi)
639 {
640 apollonianstruct *cp = &apollonians[MI_SCREEN(mi)];
641
642 if (cp->quad != NULL) {
643 (void) free((void *) cp->quad);
644 cp->quad = (apollonian_quadruple *) NULL;
645 }
646 #ifdef DOFONT
647 if (cp->gc != None) {
648 XFreeGC(display, cp->gc);
649 cp->gc = None;
650 }
651 if (cp->font != None) {
652 XFreeFont(display, cp->font);
653 cp->font = None;
654 }
655 #endif
656 }
657
658 #ifndef DEBUG
659 static void
randomize_c(int randomize,circle * c)660 randomize_c(int randomize, circle * c)
661 {
662 if (randomize / 2) {
663 double temp;
664
665 temp = c->x;
666 c->x = c->y;
667 c->y = temp;
668 }
669 if (randomize % 2) {
670 c->x = -c->x;
671 c->y = -c->y;
672 }
673 }
674 #endif
675
676 ENTRYPOINT void
init_apollonian(ModeInfo * mi)677 init_apollonian (ModeInfo * mi)
678 {
679 apollonianstruct *cp;
680 int i;
681
682 MI_INIT (mi, apollonians);
683 cp = &apollonians[MI_SCREEN(mi)];
684
685 cp->size = MAX(MIN(MI_WIDTH(mi), MI_HEIGHT(mi)) - 1, 1);
686 cp->offset.x = (MI_WIDTH(mi) - cp->size) / 2;
687 cp->offset.y = (MI_HEIGHT(mi) - cp->size) / 2;
688 cp->color_offset = NRAND(MI_NPIXELS(mi));
689
690 #ifdef DOFONT
691 if (cp->font == None) {
692 if ((cp->font = getFont(MI_DISPLAY(mi))) == None)
693 return False;
694 }
695 #endif
696 cp->label = label;
697 cp->altgeom = cp->label && altgeom;
698
699 if (cp->quad == NULL) {
700 cp->count = ABS(MI_COUNT(mi));
701 if ((cp->quad = (apollonian_quadruple *) malloc(cp->count *
702 sizeof (apollonian_quadruple))) == NULL) {
703 return;
704 }
705 dquad(cp->count, cp->quad);
706 }
707 cp->game = NRAND(PREDEF_CIRCLE_GAMES + cp->count);
708 cp->geometry = (cp->game && cp->altgeom) ? NRAND(3) : 0;
709
710 if (cp->game < PREDEF_CIRCLE_GAMES) {
711 cp->c1 = examples[cp->game][0];
712 cp->c2 = examples[cp->game][1];
713 cp->c3 = examples[cp->game][2];
714 cp->c4 = examples[cp->game][3];
715 /* do not label non int */
716 cp->label = cp->label && (cp->c4.e == (int) cp->c4.e);
717 } else { /* uses results of dquad, all int */
718 i = cp->game - PREDEF_CIRCLE_GAMES;
719 cp->c1.e = cp->quad[i].a;
720 cp->c2.e = cp->quad[i].b;
721 cp->c3.e = cp->quad[i].c;
722 cp->c4.e = cp->quad[i].d;
723 if (cp->geometry)
724 cquad(&(cp->c1), &(cp->c2), &(cp->c3), &(cp->c4));
725 }
726 cp->time = 0;
727 MI_CLEARWINDOW(mi);
728 if (cp->game != 0) {
729 double q123;
730
731 if (cp->c1.e == 0.0 || cp->c1.e == -cp->c2.e)
732 return;
733 cp->c1.x = 0.0;
734 cp->c1.y = 0.0;
735 cp->c2.x = -(cp->c1.e + cp->c2.e) / cp->c1.e;
736 cp->c2.y = 0;
737 q123 = sqrt(cp->c1.e * cp->c2.e + cp->c1.e * cp->c3.e +
738 cp->c2.e * cp->c3.e);
739 #ifdef DEBUG
740 (void) printf("q123 = %g, ", q123);
741 #endif
742 cp->c3.x = (cp->c1.e * cp->c1.e - q123 * q123) / (cp->c1.e *
743 (cp->c1.e + cp->c2.e));
744 cp->c3.y = -2.0 * q123 / (cp->c1.e + cp->c2.e);
745 q123 = -cp->c1.e - cp->c2.e + q123;
746 cp->c4.x = (cp->c1.e * cp->c1.e - q123 * q123) / (cp->c1.e *
747 (cp->c1.e + cp->c2.e));
748 cp->c4.y = -2.0 * q123 / (cp->c1.e + cp->c2.e);
749 #ifdef DEBUG
750 (void) printf("q124 = %g\n", q123);
751 (void) printf("%g %g %g %g %g %g %g %g\n",
752 cp->c1.x, cp->c1.y, cp->c2.x, cp->c2.y,
753 cp->c3.x, cp->c3.y, cp->c4.x, cp->c4.y);
754 #endif
755 }
756 #ifndef DEBUG
757 if (LRAND() & 1) {
758 cp->c3.y = -cp->c3.y;
759 cp->c4.y = -cp->c4.y;
760 }
761 i = NRAND(4);
762 randomize_c(i, &(cp->c1));
763 randomize_c(i, &(cp->c2));
764 randomize_c(i, &(cp->c3));
765 randomize_c(i, &(cp->c4));
766 #endif
767
768 mi->recursion_depth = -1;
769 }
770
771 ENTRYPOINT void
draw_apollonian(ModeInfo * mi)772 draw_apollonian (ModeInfo * mi)
773 {
774 apollonianstruct *cp;
775
776 if (apollonians == NULL)
777 return;
778 cp = &apollonians[MI_SCREEN(mi)];
779
780
781 MI_IS_DRAWN(mi) = True;
782
783 if (cp->time < 5) {
784 switch (cp->time) {
785 case 0:
786 p(mi, cp->c1);
787 p(mi, cp->c2);
788 p(mi, cp->c3);
789 p(mi, cp->c4);
790 break;
791 case 1:
792 f(mi, cp->c1, cp->c2, cp->c3, cp->c4, 0);
793 break;
794 case 2:
795 f(mi, cp->c1, cp->c2, cp->c4, cp->c3, 0);
796 break;
797 case 3:
798 f(mi, cp->c1, cp->c3, cp->c4, cp->c2, 0);
799 break;
800 case 4:
801 f(mi, cp->c2, cp->c3, cp->c4, cp->c1, 0);
802 }
803 }
804 if (++cp->time > MI_CYCLES(mi))
805 init_apollonian(mi);
806 }
807
808 XSCREENSAVER_MODULE ("Apollonian", apollonian)
809
810 #endif /* MODE_apollonian */
811