Initial import

COPYING

@@ -0,0 +1,113 @@
+CC0 1.0 Universal
+
+Statement of Purpose
+
+The laws of most jurisdictions throughout the world automatically confer
+exclusive Copyright and Related Rights (defined below) upon the creator and
+subsequent owner(s) (each and all, an "owner") of an original work of
+authorship and/or a database (each, a "Work").
+
+Certain owners wish to permanently relinquish those rights to a Work for the
+purpose of contributing to a commons of creative, cultural and scientific
+works ("Commons") that the public can reliably and without fear of later
+claims of infringement build upon, modify, incorporate in other works, reuse
+and redistribute as freely as possible in any form whatsoever and for any
+purposes, including without limitation commercial purposes. These owners may
+contribute to the Commons to promote the ideal of a free culture and the
+further production of creative, cultural and scientific works, or to gain
+reputation or greater distribution for their Work in part through the use and
+efforts of others.
+
+For these and/or other purposes and motivations, and without any expectation
+of additional consideration or compensation, the person associating CC0 with a
+Work (the "Affirmer"), to the extent that he or she is an owner of Copyright
+and Related Rights in the Work, voluntarily elects to apply CC0 to the Work
+and publicly distribute the Work under its terms, with knowledge of his or her
+Copyright and Related Rights in the Work and the meaning and intended legal
+effect of CC0 on those rights.
+
+1. Copyright and Related Rights. A Work made available under CC0 may be
+protected by copyright and related or neighboring rights ("Copyright and
+Related Rights"). Copyright and Related Rights include, but are not limited
+to, the following:
+
+  i. the right to reproduce, adapt, distribute, perform, display, communicate,
+  and translate a Work;
+
+  ii. moral rights retained by the original author(s) and/or performer(s);
+
+  iii. publicity and privacy rights pertaining to a person's image or likeness
+  depicted in a Work;
+
+  iv. rights protecting against unfair competition in regards to a Work,
+  subject to the limitations in paragraph 4(a), below;
+
+  v. rights protecting the extraction, dissemination, use and reuse of data in
+  a Work;
+
+  vi. database rights (such as those arising under Directive 96/9/EC of the
+  European Parliament and of the Council of 11 March 1996 on the legal
+  protection of databases, and under any national implementation thereof,
+  including any amended or successor version of such directive); and
+
+  vii. other similar, equivalent or corresponding rights throughout the world
+  based on applicable law or treaty, and any national implementations thereof.
+
+2. Waiver. To the greatest extent permitted by, but not in contravention of,
+applicable law, Affirmer hereby overtly, fully, permanently, irrevocably and
+unconditionally waives, abandons, and surrenders all of Affirmer's Copyright
+and Related Rights and associated claims and causes of action, whether now
+known or unknown (including existing as well as future claims and causes of
+action), in the Work (i) in all territories worldwide, (ii) for the maximum
+duration provided by applicable law or treaty (including future time
+extensions), (iii) in any current or future medium and for any number of
+copies, and (iv) for any purpose whatsoever, including without limitation
+commercial, advertising or promotional purposes (the "Waiver"). Affirmer makes
+the Waiver for the benefit of each member of the public at large and to the
+detriment of Affirmer's heirs and successors, fully intending that such Waiver
+shall not be subject to revocation, rescission, cancellation, termination, or
+any other legal or equitable action to disrupt the quiet enjoyment of the Work
+by the public as contemplated by Affirmer's express Statement of Purpose.
+
+3. Public License Fallback. Should any part of the Waiver for any reason be
+judged legally invalid or ineffective under applicable law, then the Waiver
+shall be preserved to the maximum extent permitted taking into account
+Affirmer's express Statement of Purpose. In addition, to the extent the Waiver
+is so judged Affirmer hereby grants to each affected person a royalty-free,
+non transferable, non sublicensable, non exclusive, irrevocable and
+unconditional license to exercise Affirmer's Copyright and Related Rights in
+the Work (i) in all territories worldwide, (ii) for the maximum duration
+provided by applicable law or treaty (including future time extensions), (iii)
+in any current or future medium and for any number of copies, and (iv) for any
+purpose whatsoever, including without limitation commercial, advertising or
+promotional purposes (the "License"). The License shall be deemed effective as
+of the date CC0 was applied by Affirmer to the Work. Should any part of the
+License for any reason be judged legally invalid or ineffective under
+applicable law, such partial invalidity or ineffectiveness shall not
+invalidate the remainder of the License, and in such case Affirmer hereby
+affirms that he or she will not (i) exercise any of his or her remaining
+Copyright and Related Rights in the Work or (ii) assert any associated claims
+and causes of action with respect to the Work, in either case contrary to
+Affirmer's express Statement of Purpose.
+
+4. Limitations and Disclaimers.
+
+  a. No trademark or patent rights held by Affirmer are waived, abandoned,
+  surrendered, licensed or otherwise affected by this document.
+
+  b. Affirmer offers the Work as-is and makes no representations or warranties
+  of any kind concerning the Work, express, implied, statutory or otherwise,
+  including without limitation warranties of title, merchantability, fitness
+  for a particular purpose, non infringement, or the absence of latent or
+  other defects, accuracy, or the present or absence of errors, whether or not
+  discoverable, all to the greatest extent permissible under applicable law.
+
+  c. Affirmer disclaims responsibility for clearing rights of other persons
+  that may apply to the Work or any use thereof, including without limitation
+  any person's Copyright and Related Rights in the Work. Further, Affirmer
+  disclaims responsibility for obtaining any necessary consents, permissions
+  or other rights required for any use of the Work.
+
+  d. Affirmer understands and acknowledges that Creative Commons is not a
+  party to this document and has no duty or obligation with respect to this
+  CC0 or use of the Work.

Makefile

@@ -0,0 +1,37 @@
+# fractal
+# Written in 2020 by Lucas
+# CC0 1.0 Universal/Public domain - No rights reserved
+#
+# To the extent possible under law, the author(s) have dedicated all
+# copyright and related and neighboring rights to this software to the
+# public domain worldwide. This software is distributed without any
+# warranty. You should have received a copy of the CC0 Public Domain
+# Dedication along with this software. If not, see
+# <http://creativecommons.org/publicdomain/zero/1.0/>.
+.POSIX:
+.SUFFIXES:
+.SUFFIXES: .c .o
+
+.c.o:
+	${CC} ${FRACTAL_CFLAGS} -c $<
+
+FRACTAL_CFLAGS =	-O3 -ffast-math -march=native -mtune=native \
+                        ${CFLAGS}
+
+BIN =		fractal fun-gen
+FRACTAL_OBJ =	criticals.o fractal.o julia.o palette.o pnmout.o util.o
+FUN_GEN_OBJ =	fun-gen.o
+OBJ =		${FRACTAL_OBJ} ${FUN_GEN_OBJ}
+
+all: fractal
+
+julia.o: fun.h
+
+fractal: ${FRACTAL_OBJ}
+	${CC} ${LDFLAGS} -o $@ ${FRACTAL_OBJ} -lm
+
+fun-gen: ${FUN_GEN_OBJ}
+	${CC} ${LDFLAGS} -o $@ ${FUN_GEN_OBJ} -lm
+
+clean:
+	rm -f ${OBJ} ${BIN} fun.h

criticals.c

@@ -0,0 +1,70 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "criticals.h"
+#include "util.h"
+
+size_t newton_iters = 100;
+
+static size_t idx = 0;
+static long double complex criticals[100];
+
+int
+behaves_critical(long double complex (*df)(long double complex),
+    long double complex (*ddf)(long double complex), long double complex z,
+    long double complex *crit, long double bailout, long double epsilon)
+{
+	size_t i;
+	for (i = 0; i < newton_iters && cabsl(z) < bailout; i++)
+		z = z - df(z) / ddf(z);
+
+	if (cabsl(df(z)) < epsilon) {
+		*crit = z;
+		return 1;
+	}
+
+	return 0;
+}
+
+void
+add_critical_point(long double complex z)
+{
+	size_t i;
+
+	if (idx >= nelems(criticals))
+		return;
+	for (i = 0; i < idx; i++)
+		if (cabsl(z - criticals[i]) < 1.0e-9)
+			return;
+
+	criticals[idx++] = z;
+}
+
+size_t
+ncriticals(void)
+{
+	return idx;
+}
+
+long double complex
+get_critical(size_t i)
+{
+	if (i >= idx) {
+		fprintf(stderr, "get_critical: index %zu out of range", i);
+		exit(1);
+	}
+	return criticals[i];
+}

criticals.h

@@ -0,0 +1,21 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <complex.h>
+#include <stddef.h>
+
+int behaves_critical(long double complex (*)(long double complex),
+    long double complex (*)(long double complex), long double complex,
+    long double complex *, long double, long double);
+void add_critical_point(long double complex);
+size_t ncriticals(void);
+long double complex get_critical(size_t);

fractal.c

@@ -0,0 +1,78 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <errno.h>
+#include <inttypes.h>
+#include <limits.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+#include <unistd.h>
+
+#include "util.h"
+
+int palette_main(int, char *[]);
+int julia_main(int, char *[]);
+
+static void
+usage(void)
+{
+	fprintf(stderr, "Usage: %s [-s seed] <prog> [prog_args]\n",
+	    xgetprogname());
+	exit(1);
+}
+
+int
+main(int argc, char *argv[])
+{
+	char *end;
+	intmax_t v;
+	long seed;
+	int ch;
+
+	seed = time(NULL);
+	while ((ch = getopt(argc, argv, "s:")) != -1)
+		switch (ch) {
+		case 's':
+			errno = 0;
+			v = strtoimax(optarg, &end, 0);
+			if (*end != '\0' || errno != 0 || v > LONG_MAX ||
+			    v < LONG_MIN) {
+				fprintf(stderr, "-s: Invalid value \"%s\".",
+				    optarg);
+				exit(1);
+			}
+			seed = (long)v;
+			break;
+		default:
+			usage();
+		}
+	argc -= optind;
+	argv += optind;
+	optind = 1;
+
+	if (argc == 0)
+		usage();
+
+	fprintf(stderr, "seed=%ld\n", seed);
+	xsrand48(seed);
+
+	if (strcmp(argv[0], "palette") == 0)
+		return palette_main(argc, argv);
+	else if (strcmp(argv[0], "julia") == 0)
+		return julia_main(argc, argv);
+	else
+		usage();
+
+	return 1; /* UNREACHABLE */
+}

fun-gen.c

@@ -0,0 +1,296 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <errno.h>
+#include <inttypes.h>
+#include <limits.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+#define MIN(a, b) ((a) < (b) ? (a) : (b))
+#define MAX(a, b) ((a) > (b) ? (a) : (b))
+
+static char *argv0;
+
+struct poly {
+	unsigned int degree;
+	double *coefs;
+	const char *expr;
+};
+
+static void
+poly_init(struct poly *p, unsigned int degree, const char *expr)
+{
+	if ((p->coefs = malloc((degree + 1) * sizeof(double))) == NULL) {
+		fprintf(stderr, "out of memory\n");
+		exit(1);
+	}
+
+	p->degree = degree;
+	p->expr = expr;
+}
+
+static void
+poly_derivate(struct poly p, struct poly *dp, const char *expr)
+{
+	unsigned int i;
+
+	poly_init(dp, p.degree == 0 ? 0 : p.degree - 1, expr);
+	for (i = p.degree; i > 0; i--)
+		dp->coefs[i - 1] = p.coefs[i] * i;
+}
+
+static int
+poly_cmp(const void *vp, const void *vq)
+{
+	const struct poly *p = vp, *q = vq;
+	return p->degree - q->degree;
+}
+
+static void
+poly_free(struct poly *p)
+{
+	free(p->coefs);
+}
+
+static void
+print_fun_header(const char *fun_name)
+{
+	printf("static long double complex\n");
+	printf("%s(long double complex z)\n", fun_name);
+	printf("{\n");
+}
+
+static void
+print_fun_body(size_t n, struct poly *polys)
+{
+	size_t i, j;
+	unsigned int d;
+
+	printf("\tlong double complex a");
+	for (i = 0; i < n; i++)
+		printf(", %s = 0.0L", polys[i].expr);
+	printf(";\n\n");
+
+	d = 0;
+	for (i = 0; i < n; i++) {
+		for (; d <= polys[i].degree; d++) {
+			if (d == 0)
+				printf("\ta = 1.0L;\n");
+			else
+				printf("\ta *= z;\n");
+			for (j = n; j > i; j--)
+				if (polys[j - 1].coefs[d] != 0.0L)
+					printf("\t%s += %.20f * a;\n",
+					    polys[j - 1].expr,
+					    polys[j - 1].coefs[d]);
+			printf("\n");
+		}
+	}
+}
+
+static void
+print_fun_footer(void)
+{
+	printf("}\n");
+}
+
+static void
+print_funs(struct poly p, struct poly q)
+{
+	struct poly polys[6];
+	struct poly dp, ddp, dq, ddq;
+
+	polys[0] = p;
+	polys[1] = q;
+	qsort(polys, 2, sizeof(struct poly), &poly_cmp);
+
+	print_fun_header("f");
+	print_fun_body(2, polys);
+	printf("\treturn p / q;\n");
+	print_fun_footer();
+
+	poly_derivate(p, &dp, "dp");
+	poly_derivate(q, &dq, "dq");
+	polys[0] = p;
+	polys[1] = q;
+	polys[2] = dp;
+	polys[3] = dq;
+	qsort(polys, 4, sizeof(struct poly), &poly_cmp);
+
+	printf("\n");
+	print_fun_header("df");
+	print_fun_body(4, polys);
+	printf("\treturn (dp * q - p * dq) / (q * q);\n");
+	print_fun_footer();
+
+	poly_derivate(dp, &ddp, "ddp");
+	poly_derivate(dq, &ddq, "ddq");
+	polys[0] = p;
+	polys[1] = q;
+	polys[2] = dp;
+	polys[3] = dq;
+	polys[4] = ddp;
+	polys[5] = ddq;
+	qsort(polys, 6, sizeof(struct poly), &poly_cmp);
+
+	printf("\n");
+	print_fun_header("ddf");
+	print_fun_body(6, polys);
+	printf("\treturn ("
+	    "q * q * ddp "
+	    "- q * (2.0L * dp * dq + p * ddq) "
+	    "+ 2.0L * p * dq * dq"
+	    ") / (q * q * q);\n");
+	/*printf("\treturn (ddp "
+	    "- 2.0 * ((dp * q - p * dq) / (q * q)) * dq "
+	    "- (p / q) * ddq) / q;\n");*/
+	print_fun_footer();
+
+	poly_free(&dp);
+	poly_free(&dq);
+	poly_free(&ddp);
+	poly_free(&ddq);
+}
+
+static unsigned int
+parse_uint(const char *s)
+{
+	char *end;
+	uintmax_t v;
+
+	errno = 0;
+	v = strtoumax(s, &end, 10);
+	if (s == end || *end != '\0' || errno != 0 || v > UINT_MAX) {
+		fprintf(stderr, "%s: invalid value: %s\n", argv0, s);
+		exit(1);
+	}
+
+	return (unsigned int)v;
+}
+
+static double
+parse_double(const char *s)
+{
+	char *end;
+	double d;
+
+	errno = 0;
+	d = strtod(s, &end);
+	if (s == end || *end != '\0' || errno != 0 || !isfinite(d)) {
+		fprintf(stderr, "%s: invalid value: %s\n", argv0, s);
+		exit(1);
+	}
+
+	return d;
+}
+
+static void
+usage(void)
+{
+	fprintf(stderr, "Usage: %s p-degree q-degree p-coefs q-coefs\n"
+	    "    %s -n p-degree p-coefs\n", argv0, argv0);
+	exit(1);
+}
+
+static void
+do_poly_div(int argc, char *argv[])
+{
+	struct poly p, q;
+	unsigned int i, j, p_degree, q_degree;
+
+	if (argc < 2)
+		usage();
+
+	p_degree = parse_uint(argv[0]);
+	q_degree = parse_uint(argv[1]);
+	if (argc != p_degree + q_degree + 2 + 2)
+		usage();
+
+	poly_init(&p, p_degree, "p");
+	poly_init(&q, q_degree, "q");
+
+	for (j = 0, i = 1 + 1 + p.degree + 1 + q.degree;
+	    i > 1 + 1 + p.degree;
+	    i--, j++)
+		q.coefs[j] = parse_double(argv[i]);
+	for (j = 0;
+	    i > 1;
+	    i--, j++)
+		p.coefs[j] = parse_double(argv[i]);
+
+	print_funs(p, q);
+	poly_free(&p);
+	poly_free(&q);
+}
+
+static void
+do_newton(int argc, char *argv[])
+{
+	struct poly f, df, p;
+	unsigned int i, j, f_degree;
+
+	if (argc < 1)
+		usage();
+
+	f_degree = parse_uint(argv[0]);
+	if (argc != f_degree + 1 + 1)
+		usage();
+
+	poly_init(&f, f_degree, "f");
+	for (j = 0, i = 1 + f.degree;
+	    i > 0;
+	    i--, j++)
+		f.coefs[j] = parse_double(argv[i]);
+
+	poly_derivate(f, &df, "q");
+	poly_init(&p, f_degree, "p");
+	for (i = 0; i <= f.degree; i++)
+		p.coefs[i] = (i == 0 ? 0.0 : df.coefs[i - 1]) - f.coefs[i];
+
+	print_funs(p, df);
+	poly_free(&f);
+	poly_free(&df);
+	poly_free(&p);
+}
+
+int
+main(int argc, char *argv[])
+{
+	struct poly p, q;
+	char *end;
+	uintmax_t v;
+	unsigned int i, j, p_degree, q_degree;
+	int ch, nflag;
+
+	argv0 = argv[0];
+	nflag = 0;
+	while ((ch = getopt(argc, argv, "n")) != -1)
+		switch (ch) {
+		case 'n':
+			nflag = 1;
+			break;
+		default:
+			usage();
+		}
+	argc -= optind;
+	argv += optind;
+
+	if (nflag)
+		do_newton(argc, argv);
+	else
+		do_poly_div(argc, argv);
+
+	return 0;
+}

hi-res.sh

@@ -0,0 +1,34 @@
+#!/bin/sh
+# fractal
+# Written in 2020 by Lucas
+# CC0 1.0 Universal/Public domain - No rights reserved
+#
+# To the extent possible under law, the author(s) have dedicated all
+# copyright and related and neighboring rights to this software to the
+# public domain worldwide. This software is distributed without any
+# warranty. You should have received a copy of the CC0 Public Domain
+# Dedication along with this software. If not, see
+# <http://creativecommons.org/publicdomain/zero/1.0/>.
+usage()
+{
+	printf "Usage: %s width height [fractal args...]\n" "${0##*/}" >&2
+	exit 1
+}
+
+tonumber()
+{
+	printf "%u\n" "$*"
+}
+
+width=
+height=
+if [ $# -lt 2 ] || ! width=$(tonumber "$1") || ! height=$(tonumber "$2"); then
+	usage
+fi
+shift 2
+
+upwidth=$((width * 4))
+upheight=$((height * 4))
+./fractal "$@" -w $upwidth -h $upheight |
+    pamscale -width $width -height $height |
+    pamtopng >julia.png

julia.c

@@ -0,0 +1,361 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <complex.h>
+#include <errno.h>
+#include <inttypes.h>
+#include <limits.h>
+#include <math.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+#include "criticals.h"
+#include "palette.h"
+#include "pnmout.h"
+#include "util.h"
+
+#include "fun.h"
+
+struct color contour_color = { 0 };
+static long double bailout = 1e100,
+    epsilon = 1e-8,
+    sec_w = 4.0,
+    contour = 0.0005L;
+static unsigned long niters = 250;
+
+static void
+print_critical_points(void)
+{
+	long double complex z;
+	size_t i, n;
+
+	n = ncriticals();
+	fprintf(stderr, "Critical points:\n");
+	for (i = 0; i < n; i++) {
+		z = get_critical(i);
+		fprintf(stderr, "%4zu. (%.20Lf, %.20Lf)\n", i, creall(z),
+		    cimagl(z));
+	}
+}
+
+static int
+tends_to_critical(long double complex z, long double *t)
+{
+	long double d;
+	size_t i, n;
+
+	n = ncriticals();
+	for (i = 0; i < n; i++) {
+		d = cabsl(z - get_critical(i));
+		if (d < epsilon) {
+			*t = d;
+			return 1;
+		}
+	}
+
+	return 0;
+}
+
+static long
+iterate(long double complex z, long double complex c, long double *rit)
+{
+	long double complex dz = 1.0L;
+	long double t;
+	unsigned long ci, i;
+
+	for (i = 0; i < niters && cabsl(z) < bailout; i++) {
+		dz = df(z) * dz;
+		z = f(z) + c;
+
+		if (tends_to_critical(z, &t)) {
+			*rit = logl(logl(t) / logl(epsilon)) / logl(2.0L);
+			return fabsl(logl(t)) * t < sec_w * contour * cabsl(dz)
+			    ? -1 : i;
+		}
+	}
+
+	t = cabsl(z);
+	if (t >= bailout) {
+		*rit = logl(logl(t) / logl(bailout)) / logl(2.0L);
+		return fabsl(logl(t)) * t < sec_w * contour * cabsl(dz)
+		    ? -1 : i;
+	}
+
+	return -1;
+}
+
+static long
+fast_iterate(long double complex z, long double complex c, long double *rit)
+{
+	long double t;
+	unsigned long ci, i;
+
+	for (i = 0; i < niters && cabsl(z) < bailout; i++) {
+		z = f(z) + c;
+
+		if (tends_to_critical(z, &t)) {
+			*rit = logl(logl(t) / logl(epsilon)) / logl(2.0L);
+			return i;
+		}
+	}
+
+	t = cabsl(z);
+	if (t >= bailout) {
+		*rit = logl(logl(t) / logl(bailout)) / logl(2.0L);
+		return i;
+	}
+
+	return -1;
+}
+
+static void
+usage(void)
+{
+	fprintf(stderr, "Usage: %s julia [-L] [-x julia_x] [-y julia_y]\n"
+	    "    [-n niters] [-h height] [-w width] [-z zoom]\n"
+	    "    [-a center_x] [-b center_y] [-l light_intensity]\n"
+	    "    [-d density] [-c contour]\n",
+	    xgetprogname());
+	exit(1);
+}
+
+int
+julia_main(int argc, char *argv[])
+{
+	struct color *palette;
+	struct {
+		long long it;
+		long double rit;
+	} *rescache;
+	char *end;
+	struct color color;
+	long double complex crit, julia, z;
+	long double julia_x = 0.0,
+	    julia_y = 0.0,
+	    cen_x = 0.0,
+	    cen_y = 0.0,
+	    density = 28.3L,
+	    light_intensity = 84.3L,
+	    displacement = 0,
+	    zoom = 1.0,
+	    light, light_x, light_y, light_z, rit, sec_h, zx, zy;
+	intmax_t v;
+	long long it, idx, ritnum;
+	unsigned int palette_size = 720,
+	    width = 640,
+	    height = 400,
+	    i, j;
+	int ch, do_light = 0;
+
+	while ((ch = getopt(argc, argv, "a:b:c:d:h:Ll:n:w:x:y:z:")) != -1)
+		switch (ch) {
+		case 'a':
+			errno = 0;
+			cen_x = strtold(optarg, &end);
+			if (end == optarg || *end != '\0' || errno != 0) {
+				fprintf(stderr, "-a: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			break;
+		case 'b':
+			errno = 0;
+			cen_y = strtold(optarg, &end);
+			if (end == optarg || *end != '\0' || errno != 0) {
+				fprintf(stderr, "-b: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			break;
+		case 'c':
+			errno = 0;
+			contour = strtold(optarg, &end);
+			if (end == optarg || *end != '\0' || errno != 0
+			    || contour <= 0) {
+				fprintf(stderr, "-c: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			break;
+		case 'd':
+			errno = 0;
+			density = strtold(optarg, &end);
+			if (end == optarg || *end != '\0' || errno != 0) {
+				fprintf(stderr, "-d: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			break;
+		case 'h':
+			errno = 0;
+			v = strtoimax(optarg, &end, 0);
+			if (end == optarg || *end != '\0' || errno != 0 ||
+			    v <= 0 || v >= 65536) {
+				fprintf(stderr, "-h: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			height = (unsigned int)v + 1;
+			break;
+		case 'L':
+			do_light = 1;
+			break;
+		case 'l':
+			errno = 0;
+			light_intensity = strtold(optarg, &end);
+			if (end == optarg || *end != '\0' || errno != 0) {
+				fprintf(stderr, "-l: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			break;
+		case 'n':
+			errno = 0;
+			v = strtoimax(optarg, &end, 0);
+			if (end == optarg || *end != '\0' || errno != 0 ||
+			    v <= 0 || v >= ULONG_MAX) {
+				fprintf(stderr, "-n: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			niters = (unsigned long)v;
+			break;
+		case 'w':
+			errno = 0;
+			v = strtoimax(optarg, &end, 0);
+			if (end == optarg || *end != '\0' || errno != 0 ||
+			    v <= 0 || v >= 65536) {
+				fprintf(stderr, "-w: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			width = (unsigned int)v + 1;
+			break;
+		case 'x':
+			errno = 0;
+			julia_x = strtold(optarg, &end);
+			if (end == optarg || *end != '\0' || errno != 0) {
+				fprintf(stderr, "-x: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			break;
+		case 'y':
+			errno = 0;
+			julia_y = strtold(optarg, &end);
+			if (end == optarg || *end != '\0' || errno != 0) {
+				fprintf(stderr, "-y: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			break;
+		case 'z':
+			errno = 0;
+			zoom = strtold(optarg, &end);
+			if (end == optarg || *end != '\0' || errno != 0
+			    || zoom <= 0) {
+				fprintf(stderr, "-z: invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			break;
+		default:
+			usage();
+		}
+	argc -= optind;
+	argv += optind;
+	if (argc > 1)
+		usage();
+
+	palette = palette_gen(palette_size);
+
+	if ((rescache = malloc(sizeof(*rescache) * width * height)) == NULL) {
+		fprintf(stderr, "out of memory\n");
+		exit(1);
+	}
+
+	if (do_light) {
+		long double angle_xy, angle_z;
+		angle_xy = drand48() * 2.0L * 3.14159265358979323846L;
+		angle_z = (0.2L + 0.6L * drand48()) * 3.14159265358979323846L;
+		light_x = cosl(angle_xy) * sinl(angle_z);
+		light_y = sinl(angle_xy) * sinl(angle_z);
+		light_z = cosl(angle_z);
+	}
+
+	sec_h = sec_w * height / (long double)width;
+	julia = julia_x + I * julia_y;
+
+	for (i = 0; i < height; i++) {
+		zy = cen_y + (0.5 - i / (long double)height) * sec_h;
+		for (j = 0; j < width; j++) {
+			zx = cen_x + (j / (long double)width - 0.5L) * sec_w;
+			z = zx + I * zy;
+			if (behaves_critical(df, ddf, z, &crit, bailout,
+			    epsilon))
+				add_critical_point(crit);
+		}
+	}
+
+	print_critical_points();
+
+	sec_w /= zoom; sec_h /= zoom;
+	for (i = 0; i < height; i++) {
+		zy = cen_y + (0.5L - i / (long double)height) * sec_h;
+		for (j = 0; j < width; j++) {
+			zx = cen_x + (j / (long double)width - 0.5L) * sec_w;
+			z = zx + I * zy;
+			it = iterate(z, julia, &rit);
+			rescache[i * width + j].it = it;
+			rescache[i * width + j].rit = rit;
+		}
+	}
+
+	pnmout_header(stdout, width - 1, height - 1);
+	for (i = 0; i < height - 1; i++) {
+		for (j = 0; j < width - 1; j++) {
+			it = rescache[i * width + j].it;
+			rit = rescache[i * width + j].rit;
+
+			if (it < 0)
+				color = contour_color;
+			else {
+				ritnum = (it - rit) * density;
+				if (do_light) {
+					long double lxrn, lyrn, l, lx, ly, lz;
+					lxrn = rescache[i * width + j + 1].it
+					    - rescache[i * width + j + 1].rit;
+					lyrn = rescache[(i + 1) * width + j].it
+					    - rescache[(i + 1) * width + j].rit;
+					lx = -sec_h * (lxrn - it + rit);
+					ly = sec_w * (lyrn - it + rit);
+					lz = sec_w * sec_h;
+					l = sqrtl(lx * lx + ly * ly + lz * lz);
+					light = (light_x * lx + light_y * ly
+					    + light_z * lz) / l;
+					ritnum += light * light_intensity;
+				}
+				idx = (long long)ritnum
+				    + 1000000 * palette_size + displacement;
+				color = palette[idx % palette_size
+				    + (idx < 0 ? palette_size : 0)];
+			}
+
+			pnmout_pixel(stdout, color);
+		}
+	}
+
+	free(rescache);
+	free(palette);
+
+	return !(fflush(stdout) == 0);
+}

palette.c

@@ -0,0 +1,131 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <errno.h>
+#include <inttypes.h>
+#include <limits.h>
+#include <stdlib.h>
+#include <unistd.h>
+
+#include "palette.h"
+#include "pnmout.h"
+#include "util.h"
+
+#define CLAMP(x, a, b) ((x) < (a) ? (a) : (x) > (b) ? (b) : (x))
+
+static double
+bezier3(double t, double p0, double p1, double p2, double p3)
+{
+	double u = 1.0 - t;
+	return u * u * u * p0
+	    + 3.0 * u * u * t * p1
+	    + 3.0 * u * t * t * p2
+	    + t * t * t * p3;
+}
+
+static double
+bezier4(double t, double p0, double p1, double p2, double p3)
+{
+	return (1.0 - t) * bezier3(t, p0, p1, p2, p3)
+	    + t * bezier3(t, p1, p2, p3, p0);
+}
+
+struct color *
+palette_gen(size_t n)
+{
+	struct color *palette;
+	double p[12] = {
+		drand48(), drand48(), drand48(),
+		drand48(), drand48(), drand48(),
+		drand48(), drand48(), drand48(),
+		drand48(), drand48(), drand48(),
+	};
+	double r, g, b, t;
+	size_t i;
+
+	if (!(palette = malloc(n * sizeof(*palette)))) {
+		fprintf(stderr, "out of memory");
+		exit(1);
+	}
+
+	for (i = 0; i < n; i++) {
+		t = i / (double)n;
+		r = bezier4(t, p[0], p[3], p[6], p[9]);
+		g = bezier4(t, p[1], p[4], p[7], p[10]);
+		b = bezier4(t, p[2], p[5], p[8], p[11]);
+		palette[i].r = 65535.0 * CLAMP(r, 0.0, 1.0);
+		palette[i].g = 65535.0 * CLAMP(g, 0.0, 1.0);
+		palette[i].b = 65535.0 * CLAMP(b, 0.0, 1.0);
+	}
+
+	return palette;
+}
+
+static void
+usage(void)
+{
+	fprintf(stderr, "Usage: %s palette [-h height] [-w width]\n",
+	    xgetprogname());
+	exit(1);
+}
+
+int
+palette_main(int argc, char *argv[])
+{
+	struct color *palette;
+	char *end;
+	uintmax_t v;
+	unsigned int width, height;
+	int ch, i, j;
+
+	height = 100;
+	width = 720;
+	while ((ch = getopt(argc, argv, "h:w:")) != -1)
+		switch (ch) {
+		case 'h':
+			errno = 0;
+			v = strtoumax(optarg, &end, 0);
+			if (*end != '\0' || errno != 0 || v == 0
+			    || v > UINT_MAX) {
+				fprintf(stderr, "-h: Invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			height = v;
+		case 'w':
+			errno = 0;
+			v = strtoumax(optarg, &end, 0);
+			if (*end != '\0' || errno != 0 || v == 0
+			    || v > UINT_MAX) {
+				fprintf(stderr, "-w: Invalid value \"%s\"",
+				    optarg);
+				exit(1);
+			}
+			width = v;
+		default:
+			usage();
+		}
+	argc -= optind;
+	argv += optind;
+
+	if (argc > 1)
+		usage();
+
+	palette = palette_gen(width);
+	pnmout_header(stdout, width, height);
+	for (i = 0; i < height; i++)
+		for (j = 0; j < width; j++)
+			pnmout_pixel(stdout, palette[j]);
+	free(palette);
+
+	return !(fflush(stdout) == 0);
+}

palette.h

@@ -0,0 +1,21 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <stddef.h>
+#include <stdint.h>
+
+struct color {
+	uint16_t r, g, b;
+};
+
+struct color * palette_gen(size_t);
+int palette_main(int, char *[]);

pnmout.c

@@ -0,0 +1,40 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "palette.h"
+#include "pnmout.h"
+
+void
+pnmout_header(FILE *fp, unsigned int width, unsigned int height)
+{
+	if (fprintf(fp, "P6\n%u %u 255\n", width, height) < 0) {
+		perror("fprintf");
+		exit(1);
+	}
+}
+
+void
+pnmout_pixel(FILE *fp, struct color c)
+{
+	uint8_t v[3];
+	v[0] = c.r / 256;
+	v[1] = c.g / 256;
+	v[2] = c.b / 256;
+	if (fwrite(v, sizeof(v[0]), sizeof(v) / sizeof(v[0]), fp) != 3) {
+		perror("fwrite");
+		exit(1);
+	}
+
+}

pnmout.h

@@ -0,0 +1,16 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <stdio.h>
+
+void pnmout_header(FILE *, unsigned int, unsigned int);
+void pnmout_pixel(FILE *, struct color);

util.c

@@ -0,0 +1,34 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#include <stdlib.h>
+
+const char *
+xgetprogname(void)
+{
+#if defined(__OpenBSD__)
+	return getprogname();
+#else
+	extern char *__progname;
+	return __progname;
+#endif
+}
+
+void
+xsrand48(long seed)
+{
+#if defined(__OpenBSD__)
+	srand48_deterministic(seed);
+#else
+	srand48(seed);
+#endif
+}

util.h

@@ -0,0 +1,16 @@
+/*
+ * fractal
+ * Written in 2020 by Lucas
+ * CC0 1.0 Universal/Public domain - No rights reserved
+ *
+ * To the extent possible under law, the author(s) have dedicated all
+ * copyright and related and neighboring rights to this software to the
+ * public domain worldwide. This software is distributed without any
+ * warranty. You should have received a copy of the CC0 Public Domain
+ * Dedication along with this software. If not, see
+ * <http://creativecommons.org/publicdomain/zero/1.0/>.
+ */
+#define nelems(_a) (sizeof((_a)) / sizeof((_a)[0]))
+
+const char * xgetprogname(void);
+void xsrand48(long);