From http://www.jwz.org/xscreensaver/xscreensaver-5.30.tar.gz
[xscreensaver] / hacks / apollonian.c
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 refresh_apollonian 0
79 # include "xlockmore.h"         /* in xscreensaver distribution */
80 # include "erase.h"
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         {"-altgeom", ".apollonian.altgeom", XrmoptionNoArg, "on"},
96         {"+altgeom", ".apollonian.altgeom", XrmoptionNoArg, "off"},
97         {"-label", ".apollonian.label", XrmoptionNoArg, "on"},
98         {"+label", ".apollonian.label", XrmoptionNoArg, "off"},
99 };
100 static argtype vars[] =
101 {
102         {&altgeom, "altgeom", "AltGeom", DEF_ALTGEOM, t_Bool},
103         {&label,   "label",   "Label",   DEF_LABEL,   t_Bool},
104 };
105 static OptionStruct desc[] =
106 {
107         {"-/+altgeom", "turn on/off alternate geometries (off euclidean space, on includes spherical and hyperbolic)"},
108         {"-/+label", "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", (char *) NULL, &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 enum space {
138   euclidean = 0, spherical, hyperbolic
139 };
140
141 static const char * space_string[] = {
142   "euclidean",
143   "spherical",
144   "hyperbolic"
145 };
146
147 /*
148 Generate Apollonian packing starting with a quadruple of circles.
149 The four input lines each contain the 5-tuple (e,s,h,x,y) representing
150 the circle with radius 1/e and center (x/e,y/e).  The s and h is propagated
151 like e, x and y, but can differ from e so as to represent different
152 geometries, spherical and hyperbolic, respectively.  The "standard" picture,
153 for example (-1, 2, 2, 3), can be labeled for the three geometries.
154 Origins of circles z1, z2, z3, z4
155 a * z1 = 0
156 b * z2 = (a+b)/a
157 c * z3 = (q123 + a * i)^2/(a*(a+b)) where q123 = sqrt(a*b+a*c+b*c)
158 d * z4 = (q124 + a * i)^2/(a*(a+b)) where q124 = q123 - a - b
159 If (e,x,y) represents the Euclidean circle (1/e,x/e,y/e) (so that e is
160 the label in the standard picture) then the "spherical label" is
161 (e^2+x^2+y^2-1)/(2*e) (an integer!)  and the "hyperbolic label", is 
162 calulated by h + s = e.
163 */
164 static circle examples[][4] = {
165 { /* double semi-bounded */
166         { 0, 0, 0,   0,  1},
167         { 0, 0, 0,   0, -1},
168         { 1, 1, 1,  -1,  0},
169         { 1, 1, 1,   1,  0}
170 },
171 #if 0
172 { /* standard */
173         {-1, 0, -1,   0,  0},
174         { 2, 1,  1,   1,  0},
175         { 2, 1,  1,  -1,  0},
176         { 3, 2,  1,   0,  2}
177 },
178 { /* next simplest */
179         {-2, -1, -1,   0.0,  0},
180         { 3,  2,  1,   0.5,  0},
181         { 6,  3,  3,  -2.0,  0},
182         { 7,  4,  3,  -1.5,  2}
183 },
184 { /*  */
185         {-3, -2, -1,         0.0,  0},
186         { 4,  3,  1,   1.0 / 3.0,  0},
187         {12,  7,  5,        -3.0,  0},
188         {13,  8,  5,  -8.0 / 3.0,  2}
189 },
190 { /* Mickey */
191         {-3, -2, -1,         0.0,  0},
192         { 5,  4,  1,   2.0 / 3.0,  0},
193         { 8,  5,  3,  -4.0 / 3.0, -1},
194         { 8,  5,  3,  -4.0 / 3.0,  1}
195 },
196 { /*  */
197         {-4, -3, -1,   0.00,  0},
198         { 5,  4,  1,   0.25,  0},
199         {20, 13,  7,  -4.00,  0},
200         {21, 14,  7,  -3.75,  2}
201 },
202 { /* Mickey2 */
203         {-4, -2, -2,    0.0,  0},
204         { 8,  4,  4,    1.0,  0},
205         { 9,  5,  4,  -0.75, -1},
206         { 9,  5,  4,  -0.75,  1}
207 },
208 { /* Mickey3 */
209         {-5,  -4, -1,   0.0,  0},
210         { 7,   6,  1,   0.4,  0},
211         {18,  13,  5,  -2.4, -1},
212         {18,  13,  5,  -2.4,  1}
213 },
214 { /*  */
215         {-6, -5, -1,          0.0,  0},
216         { 7,  6,  1,    1.0 / 6.0,  0},
217         {42, 31, 11,         -6.0,  0},
218         {43, 32, 11,  -35.0 / 6.0,  2}
219 },
220 { /*  */
221         {-6, -3, -3,         0.0,  0},
222         {10,  5,  5,   2.0 / 3.0,  0},
223         {15,  8,  7,        -1.5,  0},
224         {19, 10,  9,  -5.0 / 6.0,  2}
225 },
226 { /* asymmetric */
227         {-6, -5, -1,           0.0,  0.0},
228         {11, 10,  1,     5.0 / 6.0,  0.0},
229         {14, 11,  3,  -16.0 / 15.0, -0.8},
230         {15, 12,  3,          -0.9,  1.2} 
231 },
232 #endif
233 /* Non integer stuff */
234 #define DELTA 2.154700538 /* ((3+2*sqrt(3))/3) */
235 { /* 3 fold symmetric bounded (x, y calculated later) */
236         {   -1,    -1,    -1,   0.0,  0.0},
237         {DELTA, DELTA, DELTA,   1.0,  0.0},
238         {DELTA, DELTA, DELTA,   1.0, -1.0},
239         {DELTA, DELTA, DELTA,  -1.0,  1.0} 
240 },
241 { /* semi-bounded (x, y calculated later) */
242 #define ALPHA 2.618033989 /* ((3+sqrt(5))/2) */
243         {              1.0,               1.0,               1.0,   0,  0},
244         {              0.0,               0.0,               0.0,   0, -1},
245         {1.0/(ALPHA*ALPHA), 1.0/(ALPHA*ALPHA), 1.0/(ALPHA*ALPHA),  -1,  0},
246         {        1.0/ALPHA,         1.0/ALPHA,         1.0/ALPHA,  -1,  0}
247 },
248 { /* unbounded (x, y calculated later) */
249 /* #define PHI 1.618033989 *//* ((1+sqrt(5))/2) */
250 #define BETA 2.890053638 /* (PHI+sqrt(PHI)) */
251         {                 1.0,                  1.0,                  1.0,  0,  0},
252         {1.0/(BETA*BETA*BETA), 1.0/(BETA*BETA*BETA), 1.0/(BETA*BETA*BETA),  1,  0},
253         {     1.0/(BETA*BETA),      1.0/(BETA*BETA),      1.0/(BETA*BETA),  1,  0},
254         {            1.0/BETA,             1.0/BETA,             1.0/BETA,  1,  0}
255 }
256 };
257
258 #define PREDEF_CIRCLE_GAMES  (sizeof (examples) / (4 * sizeof (circle)))
259
260 #if 0
261 Euclidean
262 0, 0, 1, 1
263 -1, 2, 2, 3
264 -2, 3, 6, 7
265 -3, 5, 8, 8
266 -4, 8, 9, 9
267 -3, 4, 12, 13
268 -6, 11, 14, 15
269  -5, 7, 18, 18
270  -6, 10, 15, 19
271  -7, 12, 17, 20
272  -4, 5, 20, 21
273  -9, 18, 19, 22
274  -8, 13, 21, 24
275 Spherical
276 0, 1, 1, 2
277  -1, 2, 3, 4
278  -2, 4, 5, 5
279  -2, 3, 7, 8
280 Hyperbolic
281 -1, 1, 1, 1
282  0, 0, 1, 3
283  -2, 3, 5, 6
284  -3, 6, 6, 7
285 #endif
286
287 typedef struct {
288         int         size;
289         XPoint      offset;
290         int         geometry;
291         circle      c1, c2, c3, c4;
292         int         color_offset;
293         int         count;
294         Bool        label, altgeom;
295         apollonian_quadruple  *quad;
296 #ifdef DOFONT
297         XFontStruct *font;
298 #endif
299         int         time;
300         int         game;
301 #ifdef STANDALONE
302   eraser_state *eraser;
303 #endif
304 } apollonianstruct;
305
306 static apollonianstruct *apollonians = (apollonianstruct *) NULL;
307
308 #define FONT_HEIGHT 19
309 #define FONT_WIDTH 15
310 #define FONT_LENGTH 20
311 #define MAX_CHAR 10
312 #define K       2.15470053837925152902  /* 1+2/sqrt(3) */
313 #define MAXBEND 100 /* Do not want configurable by user since it will take too
314         much time if increased. */
315
316 static int
317 gcd(int a, int b)
318 {
319         int r;
320
321         while (b) {
322                 r = a % b;
323                 a = b;
324                 b = r;
325         }
326         return a;
327 }
328
329 static int
330 isqrt(int n)
331 {
332         int y;
333
334         if (n < 0)
335                 return -1;
336         y = (int) (sqrt((double) n) + 0.5);
337         return ((n == y*y) ? y : -1);
338 }
339
340 static void
341 dquad(int n, apollonian_quadruple *quad)
342 {
343         int a, b, c, d;
344         int counter = 0, B, C;
345         
346         for (a = 0; a < MAXBEND; a++) {
347                 B = (int) (K * a);
348                 for (b = a + 1; b <= B; b++) {
349                         C = (int) (((a + b) * (a + b)) / (4.0 * (b - a)));
350                         for (c = b; c <= C; c++) {
351                                 d = isqrt(b*c-a*(b+c));
352                                 if (d >= 0 && (gcd(a,gcd(b,c)) <= 1)) {
353                                         quad[counter].a = -a;
354                                         quad[counter].b = b;
355                                         quad[counter].c = c;
356                                         quad[counter].d = -a+b+c-2*d;
357                                         if (++counter >= n) {
358                                                 return;
359                                         }
360                                 }
361                         }
362                 }
363         }
364         (void) printf("found only %d below maximum bend of %d\n",
365                 counter, MAXBEND);
366         for (; counter < n; counter++) {
367                 quad[counter].a = -1; 
368                 quad[counter].b = 2;
369                 quad[counter].c = 2;
370                 quad[counter].d = 3;
371         }
372         return;
373 }
374
375 /*
376  * Given a Descartes quadruple of bends (a,b,c,d), with a<0, find a
377  * quadruple of circles, represented by (bend,bend*x,bend*y), such
378  * that the circles have the given bends and the bends times the
379  * centers are integers.
380  *
381  * This just performs an exaustive search, assuming that the outer
382  * circle has center in the unit square.
383  *
384  * It is always sufficient to look in {(x,y):0<=y<=x<=1/2} for the
385  * center of the outer circle, but this may not lead to a packing
386  * that can be labelled with integer spherical and hyperbolic labels.
387  * To effect the smaller search, replace FOR(a) with
388  *
389  *      for (pa = ea/2; pa <= 0; pa++) for (qa = pa; qa <= 0; qa++)
390  */
391
392 #define For(v,l,h)      for (v = l; v <= h; v++)
393 #define FOR(z)          For(p##z,lop##z,hip##z) For(q##z,loq##z,hiq##z)
394 #define H(z)            ((e##z*e##z+p##z*p##z+q##z*q##z)%2)
395 #define UNIT(z)         ((abs(e##z)-1)*(abs(e##z)-1) >= p##z*p##z+q##z*q##z)
396 #define T(z,w)          is_tangent(e##z,p##z,q##z,e##w,p##w,q##w)
397 #define LO(r,z)         lo##r##z = iceil(e##z*(r##a+1),ea)-1
398 #define HI(r,z)         hi##r##z = iflor(e##z*(r##a-1),ea)-1
399 #define B(z)            LO(p,z); HI(p,z); LO(q,z); HI(q,z)
400
401 static int
402 is_quad(int a, int b, int c, int d)
403 {
404         int s;
405
406         s = a+b+c+d;
407         return 2*(a*a+b*b+c*c+d*d) == s*s;
408 }
409
410 static Bool
411 is_tangent(int e1, int p1, int q1, int e2, int p2, int q2)
412 {
413         int dx, dy, s;
414
415         dx = p1*e2 - p2*e1;
416         dy = q1*e2 - q2*e1;
417         s = e1 + e2;
418         return dx*dx + dy*dy == s*s;
419 }
420
421 static int
422 iflor(int a, int b)
423 {
424         int q;
425
426         if (b == 0) {
427                 (void) printf("iflor: b = 0\n");
428                 return 0;
429         }
430         if (a%b == 0)
431                 return a/b;
432         q = abs(a)/abs(b);
433         return ((a<0)^(b<0)) ? -q-1 : q;
434 }
435
436 static int
437 iceil(int a, int b)
438 {
439         int q;
440
441         if (b == 0) {
442                 (void) printf("iceil: b = 0\n");
443                 return 0;
444         }
445         if (a%b == 0)
446                 return a/b;
447         q = abs(a)/abs(b);
448         return ((a<0)^(b<0)) ? -q : 1+q;
449 }
450
451 static double
452 geom(int geometry, int e, int p, int q)
453 {
454         int g = (geometry == spherical) ? -1 :
455                 (geometry == hyperbolic) ? 1 : 0;
456
457         if (g)
458                 return (e*e + (1.0 - p*p - q*q) * g) / (2.0*e);
459         (void) printf("geom: g = 0\n");
460         return e;
461 }
462
463 static void
464 cquad(circle *c1, circle *c2, circle *c3, circle *c4)
465 {
466         int ea, eb, ec, ed;
467         int pa, pb, pc, pd;
468         int qa, qb, qc, qd;
469         int lopa, lopb, lopc, lopd;
470         int hipa, hipb, hipc, hipd;
471         int loqa, loqb, loqc, loqd;
472         int hiqa, hiqb, hiqc, hiqd;
473
474         ea = (int) c1->e;
475         eb = (int) c2->e;
476         ec = (int) c3->e;
477         ed = (int) c4->e;
478         if (ea >= 0)
479                 (void) printf("ea = %d\n", ea);
480         if (!is_quad(ea,eb,ec,ed))
481                 (void) printf("Error not quad %d %d %d %d\n", ea, eb, ec, ed);
482         lopa = loqa = ea;
483         hipa = hiqa = 0;
484         FOR(a) {
485                 B(b); B(c); B(d);
486                 if (H(a) && UNIT(a)) FOR(b) {
487                         if (H(b) && T(a,b)) FOR(c) {
488                                 if (H(c) && T(a,c) && T(b,c)) FOR(d) {
489                                   if (H(d) && T(a,d) && T(b,d) && T(c,d)) {
490                                     c1->s = geom(spherical, ea, pa, qa);
491                                     c1->h = geom(hyperbolic, ea, pa, qa);
492                                     c2->s = geom(spherical, eb, pb, qb);
493                                     c2->h = geom(hyperbolic, eb, pb, qb);
494                                     c3->s = geom(spherical, ec, pc, qc);
495                                     c3->h = geom(hyperbolic, ec, pc, qc);
496                                     c4->s = geom(spherical, ed, pd, qd);
497                                     c4->h = geom(hyperbolic, ed, pd, qd);
498                                   }
499                                 }
500                         }
501                 }
502         }
503 }
504
505 static void
506 p(ModeInfo *mi, circle c)
507 {
508         apollonianstruct *cp = &apollonians[MI_SCREEN(mi)];
509         char string[10];
510         double g, e;
511         int g_width;
512
513 #ifdef DEBUG
514         (void) printf("c.e=%g c.s=%g c.h=%g  c.x=%g c.y=%g\n",
515                 c.e, c.s, c.h, c.x, c.y);
516 #endif
517         g = (cp->geometry == spherical) ? c.s : (cp->geometry == hyperbolic) ?
518                 c.h : c.e;
519         if (c.e < 0.0) {
520                 if (g < 0.0)
521                         g = -g;
522                 if (MI_NPIXELS(mi) <= 2)
523                         XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
524                                 MI_WHITE_PIXEL(mi));
525                 else
526                         XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
527                                 MI_PIXEL(mi, ((int) ((g + cp->color_offset) *
528                                         g)) % MI_NPIXELS(mi)));
529                 XDrawArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
530                         ((int) (cp->size * (-cp->c1.e) * (c.x - 1.0) /
531                                 (-2.0 * c.e) + cp->size / 2.0 + cp->offset.x)),
532                         ((int) (cp->size * (-cp->c1.e) * (c.y - 1.0) /
533                                 (-2.0 * c.e) + cp->size / 2.0 + cp->offset.y)),
534                         (int) (cp->c1.e * cp->size / c.e),
535                         (int) (cp->c1.e * cp->size / c.e), 0, 23040);
536                 if (!cp->label) {
537 #ifdef DEBUG
538                         (void) printf("%g\n", -g);
539 #endif
540                         return;
541                 }
542                 (void) sprintf(string, "%g", (g == 0.0) ? 0 : -g);
543                 if (cp->size >= 10 * FONT_WIDTH) {
544                   /* hard code these to corners */
545                   XDrawString(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
546                         ((int) (cp->size * c.x / (2.0 * c.e))) + cp->offset.x,
547                         ((int) (cp->size * c.y / (2.0 * c.e))) + FONT_HEIGHT,
548                         string, (g == 0.0) ? 1 : ((g < 10.0) ? 2 :
549                                 ((g < 100.0) ? 3 : 4)));
550                 }
551                 if (cp->altgeom && MI_HEIGHT(mi) >= 30 * FONT_WIDTH) {
552                   XDrawString(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
553                         ((int) (cp->size * c.x / (2.0 * c.e) + cp->offset.x)),
554                         ((int) (cp->size * c.y / (2.0 * c.e) + MI_HEIGHT(mi) -
555                         FONT_HEIGHT / 2)), (char *) space_string[cp->geometry],
556                         strlen(space_string[cp->geometry]));
557                 }
558                 return;
559         }
560         if (MI_NPIXELS(mi) <= 2)
561                 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_WHITE_PIXEL(mi));
562         else
563                 XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
564                         MI_PIXEL(mi, ((int) ((g + cp->color_offset) * g)) %
565                                 MI_NPIXELS(mi)));
566         if (c.e == 0.0) {
567                 if (c.x == 0.0 && c.y != 0.0) {
568                         XDrawLine(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
569                                 0, (int) ((c.y + 1.0) * cp->size / 2.0 + cp->offset.y),
570                                 MI_WIDTH(mi),
571                                 (int) ((c.y + 1.0) * cp->size / 2.0 + cp->offset.y));
572                 } else if (c.y == 0.0 && c.x != 0.0) {
573                         XDrawLine(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
574                                 (int) ((c.x + 1.0) * cp->size / 2.0 + cp->offset.x), 0,
575                                 (int) ((c.x + 1.0) * cp->size / 2.0 + cp->offset.x),
576                                 MI_HEIGHT(mi));
577                 }
578                 return;
579         }
580         e = (cp->c1.e >= 0.0) ? 1.0 : -cp->c1.e;
581         XFillArc(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
582                 ((int) (cp->size * e * (c.x - 1.0) / (2.0 * c.e) +
583                         cp->size / 2.0 + cp->offset.x)),
584                 ((int) (cp->size * e * (c.y - 1.0) / (2.0 * c.e) +
585                         cp->size / 2.0 + cp->offset.y)),
586                 (int) (e * cp->size / c.e), (int) (e * cp->size / c.e),
587                 0, 23040);
588         if (!cp->label) {
589 #ifdef DEBUG
590                 (void) printf("%g\n", g);
591 #endif
592                 return;
593         }
594         if (MI_NPIXELS(mi) <= 2)
595                 XSetForeground(MI_DISPLAY(mi), MI_GC(mi), MI_BLACK_PIXEL(mi));
596         else
597                 XSetForeground(MI_DISPLAY(mi), MI_GC(mi),
598                         MI_PIXEL(mi, ((int) ((g + cp->color_offset) * g) +
599                                 MI_NPIXELS(mi) / 2) % MI_NPIXELS(mi)));
600         g_width = (g < 10.0) ? 1: ((g < 100.0) ? 2 : 3);
601         if (c.e < e * cp->size / (FONT_LENGTH + 5 * g_width) && g < 1000.0) {
602                 (void) sprintf(string, "%g", g);
603                 XDrawString(MI_DISPLAY(mi), MI_WINDOW(mi), MI_GC(mi),
604                         ((int) (cp->size * e * c.x / (2.0 * c.e) +
605                                 cp->size / 2.0 + cp->offset.x)) -
606                                 g_width * FONT_WIDTH / 2,
607                         ((int) (cp->size * e * c.y / (2.0 * c.e) +
608                                 cp->size / 2.0 + cp->offset.y)) +
609                                 FONT_HEIGHT / 2,
610                         string, g_width);
611         }
612 }
613
614 #define BIG 7
615 static void
616 f(ModeInfo *mi, circle c1, circle c2, circle c3, circle c4, int depth)
617 {
618         apollonianstruct *cp = &apollonians[MI_SCREEN(mi)];
619         int e = (int) ((cp->c1.e >= 0.0) ? 1.0 : -cp->c1.e);
620         circle c;
621
622         if (depth > mi->recursion_depth) mi->recursion_depth = depth;
623
624         c.e = 2*(c1.e+c2.e+c3.e) - c4.e;
625         c.s = 2*(c1.s+c2.s+c3.s) - c4.s;
626         c.h = 2*(c1.h+c2.h+c3.h) - c4.h;
627         c.x = 2*(c1.x+c2.x+c3.x) - c4.x;
628         c.y = 2*(c1.y+c2.y+c3.y) - c4.y;
629         if (c.e == 0 ||
630             c.e > cp->size * e || c.x / c.e > BIG || c.y / c.e > BIG ||
631             c.x / c.e < -BIG || c.y / c.e < -BIG)
632                 return;
633         p(mi, c);
634         f(mi, c2, c3, c, c1, depth+1);
635         f(mi, c1, c3, c, c2, depth+1);
636         f(mi, c1, c2, c, c3, depth+1);
637 }
638
639 ENTRYPOINT void
640 free_apollonian (Display *display, apollonianstruct *cp)
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
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
677 init_apollonian (ModeInfo * mi)
678 {
679         apollonianstruct *cp;
680         int i;
681
682         if (apollonians == NULL) {
683                 if ((apollonians = (apollonianstruct *) calloc(MI_NUM_SCREENS(mi),
684                                              sizeof (apollonianstruct))) == NULL)
685                         return;
686         }
687         cp = &apollonians[MI_SCREEN(mi)];
688
689         cp->size = MAX(MIN(MI_WIDTH(mi), MI_HEIGHT(mi)) - 1, 1);
690         cp->offset.x = (MI_WIDTH(mi) - cp->size) / 2;
691         cp->offset.y = (MI_HEIGHT(mi) - cp->size) / 2;
692         cp->color_offset = NRAND(MI_NPIXELS(mi));
693
694 #ifdef DOFONT
695         if (cp->font == None) {
696                 if ((cp->font = getFont(MI_DISPLAY(mi))) == None)
697                         return False;
698         }
699 #endif
700         cp->label = label;
701         cp->altgeom = cp->label && altgeom;
702
703         if (cp->quad == NULL) {
704                 cp->count = ABS(MI_COUNT(mi));
705                 if ((cp->quad = (apollonian_quadruple *) malloc(cp->count *
706                         sizeof (apollonian_quadruple))) == NULL) {
707                         return;
708                 }
709                 dquad(cp->count, cp->quad);
710         }
711         cp->game = NRAND(PREDEF_CIRCLE_GAMES + cp->count);
712         cp->geometry = (cp->game && cp->altgeom) ? NRAND(3) : 0;
713
714         if (cp->game < PREDEF_CIRCLE_GAMES) {   
715                 cp->c1 = examples[cp->game][0];
716                 cp->c2 = examples[cp->game][1];
717                 cp->c3 = examples[cp->game][2];
718                 cp->c4 = examples[cp->game][3];
719                 /* do not label non int */
720                 cp->label = cp->label && (cp->c4.e == (int) cp->c4.e);
721         } else { /* uses results of dquad, all int */
722                 i = cp->game - PREDEF_CIRCLE_GAMES;
723                 cp->c1.e = cp->quad[i].a;
724                 cp->c2.e = cp->quad[i].b;
725                 cp->c3.e = cp->quad[i].c;
726                 cp->c4.e = cp->quad[i].d;
727                 if (cp->geometry)
728                         cquad(&(cp->c1), &(cp->c2), &(cp->c3), &(cp->c4));
729         }
730         cp->time = 0;
731 #ifndef STANDALONE
732         MI_CLEARWINDOW(mi);
733 #endif
734         if (cp->game != 0) {
735                 double q123;
736
737                 if (cp->c1.e == 0.0 || cp->c1.e == -cp->c2.e)
738                         return;
739                 cp->c1.x = 0.0;
740                 cp->c1.y = 0.0;
741                 cp->c2.x = -(cp->c1.e + cp->c2.e) / cp->c1.e;
742                 cp->c2.y = 0;
743                 q123 = sqrt(cp->c1.e * cp->c2.e + cp->c1.e * cp->c3.e +
744                         cp->c2.e * cp->c3.e);
745 #ifdef DEBUG
746                 (void) printf("q123 = %g, ", q123);
747 #endif
748                 cp->c3.x = (cp->c1.e * cp->c1.e - q123 * q123) / (cp->c1.e *
749                         (cp->c1.e + cp->c2.e));
750                 cp->c3.y = -2.0 * q123 / (cp->c1.e + cp->c2.e);
751                 q123 = -cp->c1.e - cp->c2.e + q123;
752                 cp->c4.x = (cp->c1.e * cp->c1.e - q123 * q123) / (cp->c1.e *
753                         (cp->c1.e + cp->c2.e));
754                 cp->c4.y = -2.0 * q123 / (cp->c1.e + cp->c2.e);
755 #ifdef DEBUG
756                 (void) printf("q124 = %g\n", q123);
757                 (void) printf("%g %g %g %g %g %g %g %g\n",
758                         cp->c1.x, cp->c1.y, cp->c2.x, cp->c2.y,
759                         cp->c3.x, cp->c3.y, cp->c4.x, cp->c4.y);
760 #endif
761         }
762 #ifndef DEBUG
763         if (LRAND() & 1) {
764                 cp->c3.y = -cp->c3.y;
765                 cp->c4.y = -cp->c4.y;
766         }
767         i = NRAND(4);
768         randomize_c(i, &(cp->c1));
769         randomize_c(i, &(cp->c2));
770         randomize_c(i, &(cp->c3));
771         randomize_c(i, &(cp->c4));
772 #endif 
773
774     mi->recursion_depth = -1;
775 }
776
777 ENTRYPOINT void
778 draw_apollonian (ModeInfo * mi)
779 {
780         apollonianstruct *cp;
781
782         if (apollonians == NULL)
783                 return;
784         cp = &apollonians[MI_SCREEN(mi)];
785
786 #ifdef STANDALONE
787     if (cp->eraser) {
788       cp->eraser = erase_window (MI_DISPLAY(mi), MI_WINDOW(mi), cp->eraser);
789       return;
790     }
791 #endif
792
793         MI_IS_DRAWN(mi) = True;
794
795         if (cp->time < 5) {
796                 switch (cp->time) {
797                 case 0:
798                         p(mi, cp->c1);
799                         p(mi, cp->c2);
800                         p(mi, cp->c3);
801                         p(mi, cp->c4);
802                         break;
803                 case 1:
804                         f(mi, cp->c1, cp->c2, cp->c3, cp->c4, 0);
805                         break;
806                 case 2:
807                         f(mi, cp->c1, cp->c2, cp->c4, cp->c3, 0);
808                         break;
809                 case 3:
810                         f(mi, cp->c1, cp->c3, cp->c4, cp->c2, 0);
811                         break;
812                 case 4:
813                         f(mi, cp->c2, cp->c3, cp->c4, cp->c1, 0);
814                 }
815         }
816         if (++cp->time > MI_CYCLES(mi))
817       {
818 #ifdef STANDALONE
819         cp->eraser = erase_window (MI_DISPLAY(mi), MI_WINDOW(mi), cp->eraser);
820 #endif /* STANDALONE */
821                 init_apollonian(mi);
822       }
823 }
824
825 ENTRYPOINT void
826 reshape_apollonian(ModeInfo * mi, int width, int height)
827 {
828   XClearWindow (MI_DISPLAY (mi), MI_WINDOW(mi));
829   init_apollonian (mi);
830 }
831
832 ENTRYPOINT Bool
833 apollonian_handle_event (ModeInfo *mi, XEvent *event)
834 {
835   if (screenhack_event_helper (MI_DISPLAY(mi), MI_WINDOW(mi), event))
836     {
837       reshape_apollonian (mi, MI_WIDTH(mi), MI_HEIGHT(mi));
838       return True;
839     }
840   return False;
841 }
842
843 ENTRYPOINT void
844 release_apollonian (ModeInfo * mi)
845 {
846         if (apollonians != NULL) {
847                 int         screen;
848
849                 for (screen = 0; screen < MI_NUM_SCREENS(mi); screen++)
850                         free_apollonian(MI_DISPLAY(mi), &apollonians[screen]);
851                 (void) free((void *) apollonians);
852                 apollonians = (apollonianstruct *) NULL;
853         }
854 }
855
856 XSCREENSAVER_MODULE ("Apollonian", apollonian)
857
858 #endif /* MODE_apollonian */