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 reshape_apollonian 0
79 # define apollonian_handle_event 0
80 # include "xlockmore.h"		/* in xscreensaver distribution */
81 #else /* STANDALONE */
82 # include "xlock.h"		/* in xlockmore distribution */
83 #endif /* STANDALONE */
84 
85 #ifdef MODE_apollonian
86 
87 #define DEF_ALTGEOM  "True"
88 #define DEF_LABEL  "True"
89 
90 static Bool altgeom;
91 static Bool label;
92 
93 static XrmOptionDescRec opts[] =
94 {
95 	{(char *) "-altgeom", (char *) ".apollonian.altgeom", XrmoptionNoArg, (caddr_t) "on"},
96 	{(char *) "+altgeom", (char *) ".apollonian.altgeom", XrmoptionNoArg, (caddr_t) "off"},
97 	{(char *) "-label", (char *) ".apollonian.label", XrmoptionNoArg, (caddr_t) "on"},
98 	{(char *) "+label", (char *) ".apollonian.label", XrmoptionNoArg, (caddr_t) "off"},
99 };
100 static argtype vars[] =
101 {
102 	{(void *) & altgeom, (char *) "altgeom", (char *) "AltGeom", (char *) DEF_ALTGEOM, t_Bool},
103 	{(void *) & label, (char *) "label", (char *) "Label", (char *) DEF_LABEL, t_Bool},
104 };
105 static OptionStruct desc[] =
106 {
107         {(char *) "-/+altgeom", (char *) "turn on/off alternate geometries (off euclidean space, on includes spherical and hyperbolic)"},
108         {(char *) "-/+label", (char *) "turn on/off alternate space and number labeling"},
109 };
110 
111 ENTRYPOINT ModeSpecOpt apollonian_opts =
112 {sizeof opts / sizeof opts[0], opts, sizeof vars / sizeof vars[0], vars, desc};
113 
114 #ifdef DOFONT
115 extern XFontStruct *getFont(Display * display);
116 #endif
117 
118 #ifdef USE_MODULES
119 ModStruct   apollonian_description =
120 {"apollonian", "init_apollonian", "draw_apollonian", "release_apollonian",
121  "init_apollonian", "init_apollonian", "free_apollonian", &apollonian_opts,
122  1000000, 64, 20, 1, 64, 1.0, "",
123  "Shows Apollonian Circles", 0, NULL};
124 
125 #endif
126 
127 typedef struct {
128 	int a, b, c, d;
129 } apollonian_quadruple;
130 
131 typedef struct {
132         double e;       /* euclidean bend */
133         double s;       /* spherical bend */
134         double h;       /* hyperbolic bend */
135         double x, y;    /* euclidean bend times euclidean position */
136 } circle;
137 
138 enum space {
139   euclidean = 0, spherical = 1, hyperbolic = 2
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 calculated 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 #ifdef WIN32
307 #define FONT_HEIGHT 15
308 #define FONT_WIDTH 10
309 #define FONT_LENGTH 16
310 #else
311 #define FONT_HEIGHT 19
312 #define FONT_WIDTH 15
313 #define FONT_LENGTH 20
314 #endif
315 #define MAX_CHAR 10
316 #define K       2.15470053837925152902  /* 1+2/sqrt(3) */
317 #define MAXBEND 100 /* Do not want configurable by user since it will take too
318 	much time if increased. */
319 
320 static int
gcd(int a,int b)321 gcd(int a, int b)
322 {
323 	int r;
324 
325 	while (b) {
326 		r = a % b;
327 		a = b;
328 		b = r;
329        	}
330 	return a;
331 }
332 
333 static int
isqrt(int n)334 isqrt(int n)
335 {
336 	int y;
337 
338 	if (n < 0)
339 		return -1;
340 	y = (int) (sqrt((double) n) + 0.5);
341 	return ((n == y*y) ? y : -1);
342 }
343 
344 static void
dquad(int n,apollonian_quadruple * quad)345 dquad(int n, apollonian_quadruple *quad)
346 {
347 	int a, b, c, d;
348 	int counter = 0, B, C;
349 
350 	for (a = 0; a < MAXBEND; a++) {
351 		B = (int) (K * a);
352 		for (b = a + 1; b <= B; b++) {
353 			C = (int) (((a + b) * (a + b)) / (4.0 * (b - a)));
354 			for (c = b; c <= C; c++) {
355 				d = isqrt(b*c-a*(b+c));
356 				if (d >= 0 && (gcd(a,gcd(b,c)) <= 1)) {
357  					quad[counter].a = -a;
358  					quad[counter].b = b;
359  					quad[counter].c = c;
360 					quad[counter].d = -a+b+c-2*d;
361 					if (++counter >= n) {
362 						return;
363 					}
364 				}
365 			}
366 		}
367 	}
368 	(void) printf("found only %d below maximum bend of %d\n",
369 		counter, MAXBEND);
370 	for (; counter < n; counter++) {
371  		quad[counter].a = -1;
372  		quad[counter].b = 2;
373  		quad[counter].c = 2;
374 		quad[counter].d = 3;
375 	}
376 	return;
377 }
378 
379 /*
380  * Given a Descartes quadruple of bends (a,b,c,d), with a<0, find a
381  * quadruple of circles, represented by (bend,bend*x,bend*y), such
382  * that the circles have the given bends and the bends times the
383  * centers are integers.
384  *
385  * This just performs an exaustive search, assuming that the outer
386  * circle has center in the unit square.
387  *
388  * It is always sufficient to look in {(x,y):0<=y<=x<=1/2} for the
389  * center of the outer circle, but this may not lead to a packing
390  * that can be labelled with integer spherical and hyperbolic labels.
391  * To effect the smaller search, replace FOR(a) with
392  *
393  *      for (pa = ea/2; pa <= 0; pa++) for (qa = pa; qa <= 0; qa++)
394  */
395 
396 #define For(v,l,h)	for (v = l; v <= h; v++)
397 #define FOR(z)		For(p##z,lop##z,hip##z) For(q##z,loq##z,hiq##z)
398 #define H(z)		((e##z*e##z+p##z*p##z+q##z*q##z)%2)
399 #define UNIT(z)		((abs(e##z)-1)*(abs(e##z)-1) >= p##z*p##z+q##z*q##z)
400 #define T(z,w)		is_tangent(e##z,p##z,q##z,e##w,p##w,q##w)
401 #define LO(r,z)		lo##r##z = iceil(e##z*(r##a+1),ea)-1
402 #define HI(r,z)		hi##r##z = iflor(e##z*(r##a-1),ea)-1
403 #define B(z)		LO(p,z); HI(p,z); LO(q,z); HI(q,z)
404 
405 static int
is_quad(int a,int b,int c,int d)406 is_quad(int a, int b, int c, int d)
407 {
408 	int s;
409 
410 	s = a+b+c+d;
411 	return 2*(a*a+b*b+c*c+d*d) == s*s;
412 }
413 
414 static Bool
is_tangent(int e1,int p1,int q1,int e2,int p2,int q2)415 is_tangent(int e1, int p1, int q1, int e2, int p2, int q2)
416 {
417 	int dx, dy, s;
418 
419 	dx = p1*e2 - p2*e1;
420 	dy = q1*e2 - q2*e1;
421 	s = e1 + e2;
422 	return dx*dx + dy*dy == s*s;
423 }
424 
425 static int
iflor(int a,int b)426 iflor(int a, int b)
427 {
428 	int q;
429 
430 	if (b == 0) {
431 		(void) printf("iflor: b = 0\n");
432 		return 0;
433 	}
434 	if (a%b == 0)
435 		return a/b;
436 	q = abs(a)/abs(b);
437 	return ((a<0)^(b<0)) ? -q-1 : q;
438 }
439 
440 static int
iceil(int a,int b)441 iceil(int a, int b)
442 {
443 	int q;
444 
445 	if (b == 0) {
446 		(void) printf("iceil: b = 0\n");
447 		return 0;
448 	}
449 	if (a%b == 0)
450 		return a/b;
451 	q = abs(a)/abs(b);
452 	return ((a<0)^(b<0)) ? -q : 1+q;
453 }
454 
455 static double
geom(int geometry,int e,int p,int q)456 geom(int geometry, int e, int p, int q)
457 {
458 	int g = (geometry == spherical) ? -1 :
459 		(geometry == hyperbolic) ? 1 : 0;
460 
461 	if (g)
462 		return (e*e + (1.0 - p*p - q*q) * g) / (2.0*e);
463 	(void) printf("geom: g = 0\n");
464 	return ((double) e);
465 }
466 
467 static void
cquad(circle * c1,circle * c2,circle * c3,circle * c4)468 cquad(circle *c1, circle *c2, circle *c3, circle *c4)
469 {
470 	int ea, eb, ec, ed;
471 	int pa, pb, pc, pd;
472 	int qa, qb, qc, qd;
473 	int lopa, lopb, lopc, lopd;
474 	int hipa, hipb, hipc, hipd;
475 	int loqa, loqb, loqc, loqd;
476 	int hiqa, hiqb, hiqc, hiqd;
477 
478 	ea = (int) c1->e;
479 	eb = (int) c2->e;
480 	ec = (int) c3->e;
481 	ed = (int) c4->e;
482 	if (ea >= 0)
483 		(void) printf("ea = %d\n", ea);
484 	if (!is_quad(ea,eb,ec,ed))
485 		(void) printf("Error not quad %d %d %d %d\n", ea, eb, ec, ed);
486 	lopa = loqa = ea;
487 	hipa = hiqa = 0;
488 	FOR(a) {
489 		B(b); B(c); B(d);
490 		if (H(a) && UNIT(a)) FOR(b) {
491 			if (H(b) && T(a,b)) FOR(c) {
492 				if (H(c) && T(a,c) && T(b,c)) FOR(d) {
493 				  if (H(d) && T(a,d) && T(b,d) && T(c,d)) {
494 				    c1->s = geom(spherical, ea, pa, qa);
495 				    c1->h = geom(hyperbolic, ea, pa, qa);
496 				    c2->s = geom(spherical, eb, pb, qb);
497 				    c2->h = geom(hyperbolic, eb, pb, qb);
498 				    c3->s = geom(spherical, ec, pc, qc);
499 				    c3->h = geom(hyperbolic, ec, pc, qc);
500 				    c4->s = geom(spherical, ed, pd, qd);
501 				    c4->h = geom(hyperbolic, ed, pd, qd);
502 				  }
503 				}
504 			}
505 		}
506 	}
507 }
508 
509 static void
p(ModeInfo * mi,circle c)510 p(ModeInfo *mi, circle c)
511 {
512 	apollonianstruct *cp = &apollonians[MI_SCREEN(mi)];
513 	char string[15];
514 	double g, e;
515 	int g_width;
516 
517 #ifdef DEBUG
518 	(void) printf("c.e=%g c.s=%g c.h=%g  c.x=%g c.y=%g\n",
519 		c.e, c.s, c.h, c.x, c.y);
520 #endif
521 	g = (cp->geometry == spherical) ? c.s : (cp->geometry == hyperbolic) ?
522 		c.h : c.e;
523 	if (c.e < 0.0) {
524 		if (g < 0.0)
525 			g = -g;
526 		if (MI_NPIXELS(mi) <= 2)
527 			XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
528 				MI_WHITE_PIXEL(mi));
529 		else
530 			XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
531 				MI_PIXEL(mi, ((int) ((g + cp->color_offset) *
532 					g)) % MI_NPIXELS(mi)));
533 		XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
534 			((int) (cp->size * (-cp->c1.e) * (c.x - 1.0) /
535 				(-2.0 * c.e) + cp->size / 2.0 + cp->offset.x)),
536 			((int) (cp->size * (-cp->c1.e) * (c.y - 1.0) /
537 				(-2.0 * c.e) + cp->size / 2.0 + cp->offset.y)),
538 			(int) (cp->c1.e * cp->size / c.e),
539 			(int) (cp->c1.e * cp->size / c.e), 0, 23040);
540 		if (!cp->label) {
541 #ifdef DEBUG
542 			(void) printf("%g\n", -g);
543 #endif
544 			return;
545 		}
546 		(void) sprintf(string, "%g", (g == 0.0) ? 0 : -g);
547 		if (cp->size >= 10 * FONT_WIDTH) {
548 		  /* hard code these to corners */
549 		  XDrawString(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
550 			((int) (cp->size * c.x / (2.0 * c.e))) + cp->offset.x,
551 			((int) (cp->size * c.y / (2.0 * c.e))) + FONT_HEIGHT,
552 			string, (g == 0.0) ? 1 : ((g < 10.0) ? 2 :
553 				((g < 100.0) ? 3 : 4)));
554 		}
555 		if (cp->altgeom && MI_HEIGHT(mi) >= 30 * FONT_WIDTH) {
556 		  XDrawString(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
557 			((int) (cp->size * c.x / (2.0 * c.e) + cp->offset.x)),
558 			((int) (cp->size * c.y / (2.0 * c.e) + MI_HEIGHT(mi) -
559 			FONT_HEIGHT / 2)), (char *) space_string[(int) cp->geometry],
560 			strlen(space_string[(int) cp->geometry]));
561 		}
562 		return;
563 	}
564 	if (MI_NPIXELS(mi) <= 2)
565 		XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
566 	else
567 		XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
568 			MI_PIXEL(mi, ((int) ((g + cp->color_offset) * g)) %
569 				MI_NPIXELS(mi)));
570 	if (c.e == 0.0) {
571 		if (c.x == 0.0 && c.y != 0.0) {
572 			XDrawLine(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
573 				0, (int) ((c.y + 1.0) * cp->size / 2.0 + cp->offset.y),
574 				MI_WIDTH(mi),
575 				(int) ((c.y + 1.0) * cp->size / 2.0 + cp->offset.y));
576 		} else if (c.y == 0.0 && c.x != 0.0) {
577 			XDrawLine(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
578 				(int) ((c.x + 1.0) * cp->size / 2.0 + cp->offset.x), 0,
579 				(int) ((c.x + 1.0) * cp->size / 2.0 + cp->offset.x),
580 				MI_HEIGHT(mi));
581 		}
582 		return;
583 	}
584 	e = (cp->c1.e >= 0.0) ? 1.0 : -cp->c1.e;
585 	XFillArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
586 		((int) (cp->size * e * (c.x - 1.0) / (2.0 * c.e) +
587 			cp->size / 2.0 + cp->offset.x)),
588 		((int) (cp->size * e * (c.y - 1.0) / (2.0 * c.e) +
589 			cp->size / 2.0 + cp->offset.y)),
590 		(int) (e * cp->size / c.e), (int) (e * cp->size / c.e),
591 		0, 23040);
592 	if (!cp->label) {
593 #ifdef DEBUG
594 		(void) printf("%g\n", g);
595 #endif
596 		return;
597 	}
598 	if (MI_NPIXELS(mi) <= 2)
599 		XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_BLACK_PIXEL(mi));
600 	else
601 		XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
602 			MI_PIXEL(mi, ((int) ((g + cp->color_offset) * g) +
603 				MI_NPIXELS(mi) / 2) % MI_NPIXELS(mi)));
604 	g_width = (g < 10.0) ? 1: ((g < 100.0) ? 2 : 3);
605 	if (c.e < e * cp->size / (FONT_LENGTH + 5 * g_width) && g < 1000.0) {
606 		(void) sprintf(string, "%g", g);
607 		XDrawString(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
608 			((int) (cp->size * e * c.x / (2.0 * c.e) +
609 				cp->size / 2.0 + cp->offset.x)) -
610 				g_width * FONT_WIDTH / 2,
611 			((int) (cp->size * e * c.y / (2.0 * c.e) +
612 				cp->size / 2.0 + cp->offset.y)) +
613 				FONT_HEIGHT / 2,
614 			string, g_width);
615         }
616 }
617 
618 #define BIG 7
619 static void
f(ModeInfo * mi,circle c1,circle c2,circle c3,circle c4,int depth)620 f(ModeInfo *mi, circle c1, circle c2, circle c3, circle c4, int depth)
621 {
622 	apollonianstruct *cp = &apollonians[MI_SCREEN(mi)];
623 	int e = (int) ((cp->c1.e >= 0.0) ? 1.0 : -cp->c1.e);
624         circle c;
625 
626 	if (depth > MI_RECURSION_DEPTH(mi)) MI_RECURSION_DEPTH(mi) = depth;
627 
628         c.e = 2*(c1.e+c2.e+c3.e) - c4.e;
629         c.s = 2*(c1.s+c2.s+c3.s) - c4.s;
630         c.h = 2*(c1.h+c2.h+c3.h) - c4.h;
631         c.x = 2*(c1.x+c2.x+c3.x) - c4.x;
632         c.y = 2*(c1.y+c2.y+c3.y) - c4.y;
633         if (c.e == 0 ||
634 		c.e > cp->size * e || c.x / c.e > BIG || c.y / c.e > BIG ||
635             c.x / c.e < -BIG || c.y / c.e < -BIG)
636                 return;
637         p(mi, c);
638         f(mi, c2, c3, c, c1, depth+1);
639         f(mi, c1, c3, c, c2, depth+1);
640         f(mi, c1, c2, c, c3, depth+1);
641 }
642 
643 static void
free_apollonian_screen(Display * display,apollonianstruct * cp)644 free_apollonian_screen(
645 #ifdef DOFONT
646 		Display *display,
647 #endif
648 		apollonianstruct *cp)
649 {
650 	if (cp == NULL) {
651 		return;
652 	}
653 	if (cp->quad != NULL) {
654 		free(cp->quad);
655 		cp->quad = (apollonian_quadruple *) NULL;
656 	}
657 #ifdef DOFONT
658 	if (cp->gc != None) {
659 		XFreeGC(display, cp->gc);
660 		cp->gc = None;
661 	}
662 	if (cp->font != None) {
663 		XFreeFont(display, cp->font);
664 		cp->font = None;
665 	}
666 #endif
667 	cp = NULL;
668 }
669 
670 ENTRYPOINT void
free_apollonian(ModeInfo * mi)671 free_apollonian(ModeInfo * mi)
672 {
673 	apollonianstruct *ap = &apollonians[MI_SCREEN(mi)];
674 	free_apollonian_screen(
675 #ifdef DOFONT
676 		MI_DISPLAY(mi),
677 #endif
678 		ap);
679 }
680 
681 #ifndef DEBUG
682 static void
randomize_c(int randomize,circle * c)683 randomize_c(int randomize, circle * c)
684 {
685   if (randomize / 2) {
686     double temp;
687 
688     temp = c->x;
689     c->x = c->y;
690     c->y = temp;
691   }
692   if (randomize % 2) {
693     c->x = -c->x;
694     c->y = -c->y;
695   }
696 }
697 #endif
698 
699 ENTRYPOINT void
init_apollonian(ModeInfo * mi)700 init_apollonian(ModeInfo * mi)
701 {
702 	apollonianstruct *cp;
703 	int i;
704 
705 	MI_INIT (mi, apollonians);
706 	cp = &apollonians[MI_SCREEN(mi)];
707 
708 	cp->size = MAX(MIN(MI_WIDTH(mi), MI_HEIGHT(mi)) - 1, 1);
709 	cp->offset.x = (MI_WIDTH(mi) - cp->size) / 2;
710 	cp->offset.y = (MI_HEIGHT(mi) - cp->size) / 2;
711 	cp->color_offset = NRAND(MI_NPIXELS(mi));
712 
713 #ifdef DOFONT
714 	if (cp->font == None) {
715 		if ((cp->font = getFont(MI_DISPLAY(mi))) == None)
716 			return False;
717 	}
718 #endif
719 	cp->label = label;
720 	cp->altgeom = cp->label && altgeom;
721 
722 	if (cp->quad == NULL) {
723 		if (MI_COUNT(mi))
724 			cp->count = ABS(MI_COUNT(mi));
725 		else
726 			cp->count = 1;
727 		if ((cp->quad = (apollonian_quadruple *) malloc(cp->count *
728 			sizeof (apollonian_quadruple))) == NULL) {
729 			return;
730 		}
731 		dquad(cp->count, cp->quad);
732 	}
733 	cp->game = NRAND(PREDEF_CIRCLE_GAMES + cp->count);
734 	cp->geometry = (int) ((cp->game && cp->altgeom) ? NRAND(3) : 0);
735 
736 	if (cp->game < (int) PREDEF_CIRCLE_GAMES) {
737 		cp->c1 = examples[cp->game][0];
738 		cp->c2 = examples[cp->game][1];
739 		cp->c3 = examples[cp->game][2];
740 		cp->c4 = examples[cp->game][3];
741 		/* do not label non int */
742 		cp->label = cp->label && (cp->c4.e == (int) cp->c4.e);
743 	} else { /* uses results of dquad, all int */
744 		i = cp->game - PREDEF_CIRCLE_GAMES;
745 		cp->c1.e = cp->quad[i].a;
746 		cp->c2.e = cp->quad[i].b;
747 		cp->c3.e = cp->quad[i].c;
748 		cp->c4.e = cp->quad[i].d;
749 		if (cp->geometry != euclidean)
750 			cquad(&(cp->c1), &(cp->c2), &(cp->c3), &(cp->c4));
751 	}
752 	cp->time = 0;
753 	MI_CLEARWINDOW(mi);
754 	if (cp->game != 0) {
755 		double q123;
756 
757 		if (cp->c1.e == 0.0 || cp->c1.e == -cp->c2.e)
758 			return;
759 		cp->c1.x = 0.0;
760 		cp->c1.y = 0.0;
761 		cp->c2.x = -(cp->c1.e + cp->c2.e) / cp->c1.e;
762 		cp->c2.y = 0;
763 		q123 = sqrt(cp->c1.e * cp->c2.e + cp->c1.e * cp->c3.e +
764 			cp->c2.e * cp->c3.e);
765 #ifdef DEBUG
766 		(void) printf("q123 = %g, ", q123);
767 #endif
768 		cp->c3.x = (cp->c1.e * cp->c1.e - q123 * q123) / (cp->c1.e *
769 			(cp->c1.e + cp->c2.e));
770 		cp->c3.y = -2.0 * q123 / (cp->c1.e + cp->c2.e);
771 		q123 = -cp->c1.e - cp->c2.e + q123;
772 		cp->c4.x = (cp->c1.e * cp->c1.e - q123 * q123) / (cp->c1.e *
773 			(cp->c1.e + cp->c2.e));
774 		cp->c4.y = -2.0 * q123 / (cp->c1.e + cp->c2.e);
775 #ifdef DEBUG
776 		(void) printf("q124 = %g\n", q123);
777 		(void) printf("%g %g %g %g %g %g %g %g\n",
778 			cp->c1.x, cp->c1.y, cp->c2.x, cp->c2.y,
779 			cp->c3.x, cp->c3.y, cp->c4.x, cp->c4.y);
780 #endif
781 	}
782 #ifndef DEBUG
783 	if (LRAND() & 1) {
784 		cp->c3.y = -cp->c3.y;
785 		cp->c4.y = -cp->c4.y;
786 	}
787 	i = NRAND(4);
788 	randomize_c(i, &(cp->c1));
789 	randomize_c(i, &(cp->c2));
790 	randomize_c(i, &(cp->c3));
791 	randomize_c(i, &(cp->c4));
792 #endif
793 	MI_RECURSION_DEPTH(mi) = -1;
794 }
795 
796 ENTRYPOINT void
draw_apollonian(ModeInfo * mi)797 draw_apollonian(ModeInfo * mi)
798 {
799 	apollonianstruct *cp;
800 
801 	if (apollonians == NULL)
802 		return;
803 	cp = &apollonians[MI_SCREEN(mi)];
804 
805 
806 	MI_IS_DRAWN(mi) = True;
807 
808 	if (cp->time < 5) {
809 		switch (cp->time) {
810 		case 0:
811 			p(mi, cp->c1);
812 			p(mi, cp->c2);
813 			p(mi, cp->c3);
814 			p(mi, cp->c4);
815 			break;
816 		case 1:
817 			f(mi, cp->c1, cp->c2, cp->c3, cp->c4, 0);
818 			break;
819 		case 2:
820 			f(mi, cp->c1, cp->c2, cp->c4, cp->c3, 0);
821 			break;
822 		case 3:
823 			f(mi, cp->c1, cp->c3, cp->c4, cp->c2, 0);
824 			break;
825 		case 4:
826 			f(mi, cp->c2, cp->c3, cp->c4, cp->c1, 0);
827 		}
828 	}
829 	if (++cp->time > MI_CYCLES(mi))
830 		init_apollonian(mi);
831 }
832 
833 ENTRYPOINT void
release_apollonian(ModeInfo * mi)834 release_apollonian(ModeInfo * mi)
835 {
836 	if (apollonians != NULL) {
837 		int         screen;
838 
839 		for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
840 			free_apollonian_screen(
841 #ifdef DOFONT
842 				MI_DISPLAY(mi),
843 #endif
844 				&apollonians[screen]);
845 		free(apollonians);
846 		apollonians = (apollonianstruct *) NULL;
847 	}
848 }
849 
850 XSCREENSAVER_MODULE ("Apollonian", apollonian)
851 
852 #endif /* MODE_apollonian */
853