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