+#ifdef USE_FAST_SQRT_HACKISH
+/* Based loosely on code from Wikipedia:
+ * https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Approximations_that_depend_on_IEEE_representation
+ */
+
+/* FAST_SQRT_EXTRA_BITS = 3: Smallest useful value
+ * = 5/6: A little bit of banding, wave_height table is on par with regular
+ * sqrt() code.
+ * = 7: No apparent difference with original @ radius = 800.
+ * = 8: One more just to be comfortable.
+ */
+
+# define FAST_SQRT_EXTRA_BITS 8
+
+union int_float
+{
+ uint32_t i;
+ float f;
+};
+
+static unsigned fast_log2(unsigned x)
+{
+ union int_float u;
+ if(!x)
+ return x;
+ u.f = x;
+ return ((u.i - 0x3f800000) >> (23 - FAST_SQRT_EXTRA_BITS)) + 1;
+}
+
+static float fast_inv_log2(unsigned x)
+{
+ union int_float u;
+ if(!x)
+ return 0.0f;
+ u.i = ((x - 1) << (23 - FAST_SQRT_EXTRA_BITS)) + 0x3f800000;
+ return u.f;
+}
+
+#endif
+
+#ifdef USE_FAST_SQRT_BIGTABLE2
+
+/* I eyeballed these figures. They could be improved. - D.O. */
+
+# define FAST_SQRT_DISCARD_BITS1 4
+/* = 5: Dot in center is almost invisible at radius = 800. */
+/* = 4: Dot in center looks OK at radius = 50. */
+
+/* 156/2740/9029 */
+/* # define FAST_SQRT_DISCARD_BITS2 8 */
+/* # define FAST_SQRT_CUTOFF 64 * 64 */
+
+/* 156/2242/5386 */
+# define FAST_SQRT_DISCARD_BITS2 9
+# define FAST_SQRT_CUTOFF 128 * 128
+
+/*
+ * This is a little faster:
+ * 44.5 FPS, 19/5000/17578
+ *
+ * # define FAST_SQRT_DISCARD_BITS1 7
+ * # define FAST_SQRT_DISCARD_BITS2 7
+ * # define FAST_SQRT_CUTOFF 0
+ *
+ * For radius = 800, FAST_SQRT_DISCARD_BITS2 =
+ * = 9/10: Approximately the original table size, some banding near origins.
+ * = 7: wave_height is 20 KB, and just fits inside a 32K L1 cache.
+ * = 6: Nearly indistinguishable from original
+ */
+
+/*
+ FAST_TABLE(x) is equivalent to, but slightly faster than:
+ x < FAST_SQRT_CUTOFF ?
+ (x >> FAST_SQRT_DISCARD_BITS1) :
+ ((x - FAST_SQRT_CUTOFF) >> FAST_SQRT_DISCARD_BITS2) +
+ (FAST_SQRT_CUTOFF >> FAST_SQRT_DISCARD_BITS1);
+*/
+
+#define FAST_TABLE(x) \
+ ((x) < FAST_SQRT_CUTOFF ? \
+ ((x) >> FAST_SQRT_DISCARD_BITS1) : \
+ (((x) + \
+ ((FAST_SQRT_CUTOFF << (FAST_SQRT_DISCARD_BITS2 - \
+ FAST_SQRT_DISCARD_BITS1)) - FAST_SQRT_CUTOFF)) >> \
+ FAST_SQRT_DISCARD_BITS2))
+
+static double fast_inv_table(unsigned x)
+{
+ return x < (FAST_SQRT_CUTOFF >> FAST_SQRT_DISCARD_BITS1) ?
+ (x << FAST_SQRT_DISCARD_BITS1) :
+ ((x - (FAST_SQRT_CUTOFF >> FAST_SQRT_DISCARD_BITS1)) <<
+ FAST_SQRT_DISCARD_BITS2) + FAST_SQRT_CUTOFF;
+}
+
+#endif
+
+static void destroy_image(Display* dpy, struct inter_context* c)
+{
+#ifdef USE_XIMAGE
+ if(c->ximage) {
+ destroy_xshm_image(dpy, c->ximage, &c->shm_info);
+ }
+#endif
+
+ if(c->threadpool.count)
+ {
+ threadpool_destroy(&c->threadpool);
+ c->threadpool.count = 0;
+ }
+}
+
+static void inter_free(Display* dpy, struct inter_context* c)
+{
+#ifndef USE_XIMAGE
+ unsigned i;
+#endif
+
+ if(c->pix_buf)
+ XFreePixmap(dpy, c->pix_buf);
+
+ if(c->copy_gc)
+ XFreeGC(dpy, c->copy_gc);
+
+ destroy_image(dpy, c);
+
+ if(c->colors <= 2)
+ free(c->pal);
+ else if(c->pal)
+ free_colors(c->screen, c->cmap, c->pal, c->colors);
+
+#ifndef USE_XIMAGE
+ for(i = 0; i != c->colors; ++i)
+ XFreeGC(dpy, c->gcs[i]);
+ free(c->gcs);
+#endif
+
+ free(c->wave_height);
+ free(c->source);
+}
+
+static void abort_on_error(int error)
+{
+ fprintf(stderr, "interference: %s\n", strerror(error));
+ exit(1);
+}
+
+static void abort_no_mem(void)
+{
+ abort_on_error(ENOMEM);
+}
+
+static void check_no_mem(Display* dpy, struct inter_context* c, void* ptr)
+{
+ if(!ptr) {
+ inter_free(dpy, c);
+ abort_no_mem();
+ }
+}
+
+static int inter_thread_create(
+ void* self_raw,
+ struct threadpool* pool,
+ unsigned id)
+{
+ struct inter_thread* self = (struct inter_thread*)self_raw;
+ const struct inter_context* c = GET_PARENT_OBJ(struct inter_context, threadpool, pool);
+
+ self->context = c;
+ self->thread_id = id;
+
+ self->result_row = malloc(c->w_div_g * sizeof(unsigned));
+ if(!self->result_row)
+ return ENOMEM;
+
+#ifdef USE_XIMAGE
+ self->row = malloc(c->w_div_g * sizeof(uint32_t));
+ if(!self->row) {
+ free(self->result_row);
+ return ENOMEM;
+ }
+#endif
+
+ return 0;
+}
+
+static void inter_thread_destroy(void* self_raw)
+{
+ struct inter_thread* self = (struct inter_thread*)self_raw;
+#ifdef USE_XIMAGE
+ free(self->row);
+#endif
+ free(self->result_row);
+}
+
+/*
+A higher performance design would have input and output queues, so that when
+worker threads finish with one frame, they can pull the next work order from
+the queue and get started on it immediately, rather than going straight to
+sleep. The current "single-buffered" design still provides reasonable
+performance at low frame rates; high frame rates are noticeably less efficient.
+*/
+
+static void inter_thread_run(void* self_raw)
+{
+ struct inter_thread* self = (struct inter_thread*)self_raw;
+ const struct inter_context* c = self->context;
+
+ int i, j, k;
+ unsigned result;
+ int dist1;
+ int g = c->grid_size;
+
+ int dx, dy, g2 = 2 * g * g;
+ int px, py, px2g;
+
+ int dist0, ddist;
+
+#ifdef USE_XIMAGE
+ unsigned img_y = g * self->thread_id;
+ void *scanline = c->ximage->data + c->ximage->bytes_per_line * g * self->thread_id;
+#endif
+
+ for(j = self->thread_id; j < c->h_div_g; j += c->threadpool.count) {
+ px = g/2;
+ py = j*g + px;
+
+ memset(self->result_row, 0, c->w_div_g * sizeof(unsigned));
+
+ for(k = 0; k < c->count; k++) {
+
+ dx = px - c->source[k].x;
+ dy = py - c->source[k].y;
+
+ dist0 = dx*dx + dy*dy;
+ ddist = -2 * g * c->source[k].x;
+
+ /* px2g = g*(px*2 + g); */
+ px2g = g2;
+
+ for(i = 0; i < c->w_div_g; i++) {
+ /*
+ * Discarded possibilities for improving performance here:
+ * 1. Using octagon-based distance estimation
+ * (Which causes giant octagons to appear.)
+ * 2. Square root approximation by reinterpret-casting IEEE floats to
+ * integers.
+ * (Which causes angles to appear when two waves interfere.)
+ */
+
+/* int_float u;
+ u.f = dx*dx + dy*dy;
+ u.i = (1 << 29) + (u.i >> 1) - (1 << 22);
+ dist = u.f; */
+
+#if defined USE_FAST_SQRT_BIGTABLE2
+ dist1 = FAST_TABLE(dist0);
+#elif defined USE_FAST_SQRT_HACKISH
+ dist1 = fast_log2(dist0);
+#else
+ dist1 = sqrt(dist0);
+#endif
+
+ if(dist1 < c->radius)
+ self->result_row[i] += c->wave_height[dist1];
+
+ dist0 += px2g + ddist;
+ px2g += g2;
+ }
+ }
+
+ for(i = 0; i < c->w_div_g; i++) {
+
+ result = self->result_row[i];
+
+ /* It's slightly faster to do a subtraction or two before calculating the
+ * modulus. - D.O. */
+ if(result >= c->colors)
+ {
+ result -= c->colors;
+ if(result >= c->colors)
+ result %= (unsigned)c->colors;
+ }
+
+#ifdef USE_XIMAGE
+ self->row[i] = c->pal[result].pixel;
+#else
+ XFillRectangle(c->dpy, TARGET(c), c->gcs[result], g*i, g*j, g, g);
+#endif /* USE_XIMAGE */
+ }
+
+#ifdef USE_XIMAGE
+ /* Fill in these `gridsize' horizontal bits in the scanline */
+ if(c->ximage->bits_per_pixel == 32)
+ {
+ uint32_t *ptr = (uint32_t *)scanline;
+ for(i = 0; i < c->w_div_g; i++) {
+ for(k = 0; k < g; k++)
+ ptr[g*i+k] = self->row[i];
+ }
+ }
+ else if(c->ximage->bits_per_pixel == 24)
+ {
+ uint8_t *ptr = (uint8_t *)scanline;
+ for(i = 0; i < c->w_div_g; i++) {
+ for(k = 0; k < g; k++) {
+ uint32_t pixel = self->row[i];
+ /* Might not work on big-endian. */
+ ptr[0] = pixel;
+ ptr[1] = (pixel & 0x0000ff00) >> 8;
+ ptr[2] = (pixel & 0x00ff0000) >> 16;
+ ptr += 3;
+ }
+ }
+ }
+ else if(c->ximage->bits_per_pixel == 16)
+ {
+ uint16_t *ptr = (uint16_t *)scanline;
+ for(i = 0; i < c->w_div_g; i++) {
+ for(k = 0; k < g; k++)
+ ptr[g*i+k] = self->row[i];
+ }
+ }
+ else if(c->ximage->bits_per_pixel == 8)
+ {
+ uint8_t *ptr = (uint8_t *)scanline;
+ for(i = 0; i < c->w_div_g; i++) {
+ for(k = 0; k < g; k++)
+ ptr[g*i+k] = self->row[i];
+ }
+ }
+ else
+ {
+ for(i = 0; i < c->w_div_g; i++) {
+ for(k = 0; k < g; k++)
+ /* XPutPixel is thread safe as long as the XImage didn't have its
+ * bits_per_pixel changed. */
+ XPutPixel(c->ximage, (g*i)+k, img_y, self->row[i]);
+ }
+ }
+
+ /* Only the first scanline of the image has been filled in; clone that
+ scanline to the rest of the `gridsize' lines in the ximage */
+ for(k = 0; k < (g-1); k++)
+ memcpy(c->ximage->data + (c->ximage->bytes_per_line * (img_y + k + 1)),
+ c->ximage->data + (c->ximage->bytes_per_line * img_y),
+ c->ximage->bytes_per_line);
+
+ scanline = (char *)scanline +
+ c->ximage->bytes_per_line * g * c->threadpool.count;
+ img_y += g * c->threadpool.count;
+
+#endif /* USE_XIMAGE */
+ }
+}
+
+/* On allocation error, c->ximage == NULL. */
+static void create_image(
+ Display* dpy,
+ struct inter_context* c,
+ const XWindowAttributes* xgwa)
+{
+#ifdef USE_XIMAGE
+
+ /* Set the width so that each thread can work on a different line. */
+ unsigned align = thread_memory_alignment(dpy) * 8 - 1;
+ unsigned wbits, w, h;
+#endif /* USE_XIMAGE */
+
+ c->w = xgwa->width;
+ c->h = xgwa->height;
+ c->w_div_g = (c->w + c->grid_size - 1) / c->grid_size;
+ c->h_div_g = (c->h + c->grid_size - 1) / c->grid_size;
+
+#ifdef USE_XIMAGE
+ w = c->w_div_g * c->grid_size;
+ h = c->h_div_g * c->grid_size;
+
+ /* The width of a scan line, in *bits*. */
+ wbits = (w * c->bits_per_pixel + align) & ~align;
+
+ /* This uses a lot more RAM than the single line approach. Users without
+ * enough RAM to fit even a single framebuffer should consider an upgrade for
+ * their 386. - D.O.
+ */
+
+ c->ximage = create_xshm_image(dpy, xgwa->visual, xgwa->depth,
+ ZPixmap, &c->shm_info,
+ wbits / c->bits_per_pixel, h);
+
+ c->shm_can_draw = True;
+
+ check_no_mem(dpy, c, c->ximage);
+#endif /* USE_XIMAGE */
+
+ {
+ static const struct threadpool_class cls =
+ {
+ sizeof(struct inter_thread),
+ inter_thread_create,
+ inter_thread_destroy
+ };
+
+ int error = threadpool_create(
+ &c->threadpool,
+ &cls,
+ dpy,
+#ifdef USE_XIMAGE
+ hardware_concurrency(dpy)
+#else
+ 1
+ /* At least two issues with threads without USE_XIMAGE:
+ * 1. Most of Xlib isn't thread safe without XInitThreads.
+ * 2. X(Un)LockDisplay would need to be called for each line, which is
+ * terrible.
+ */
+#endif
+ );
+
+ if(error) {
+ c->threadpool.count = 0; /* See the note in thread_util.h. */
+ inter_free(dpy, c);
+ abort_on_error(error);
+ }
+ }
+}
+
+static void create_pix_buf(Display* dpy, Window win, struct inter_context *c,
+ const XWindowAttributes* xgwa)
+{
+#ifdef HAVE_DOUBLE_BUFFER_EXTENSION
+ if(c->back_buf)
+ return;
+#endif /* HAVE_DOUBLE_BUFFER_EXTENSION */
+ c->pix_buf = XCreatePixmap(dpy, win, xgwa->width, xgwa->height, xgwa->depth);
+}
+
+static double float_time(void)
+{
+ struct timeval result;
+ gettimeofday(
+ &result
+#ifdef GETTIMEOFDAY_TWO_ARGS
+ , NULL
+#endif
+ );
+
+ return result.tv_usec * 1.0e-6 + result.tv_sec;
+}
+