From 4304edc0f0477fdbfd7420796c34c34957a33a3e Mon Sep 17 00:00:00 2001 From: Lucas Date: Sat, 7 Mar 2020 22:32:57 +0000 Subject: [PATCH] Initial import --- COPYING | 113 ++++++++++++++++ Makefile | 37 ++++++ criticals.c | 70 ++++++++++ criticals.h | 21 +++ fractal.c | 78 ++++++++++++ fun-gen.c | 296 ++++++++++++++++++++++++++++++++++++++++++ hi-res.sh | 34 +++++ julia.c | 361 ++++++++++++++++++++++++++++++++++++++++++++++++++++ palette.c | 131 +++++++++++++++++++ palette.h | 21 +++ pnmout.c | 40 ++++++ pnmout.h | 16 +++ util.c | 34 +++++ util.h | 16 +++ 14 files changed, 1268 insertions(+) create mode 100644 COPYING create mode 100644 Makefile create mode 100644 criticals.c create mode 100644 criticals.h create mode 100644 fractal.c create mode 100644 fun-gen.c create mode 100644 hi-res.sh create mode 100644 julia.c create mode 100644 palette.c create mode 100644 palette.h create mode 100644 pnmout.c create mode 100644 pnmout.h create mode 100644 util.c create mode 100644 util.h diff --git a/COPYING b/COPYING new file mode 100644 index 0000000..c7a74f9 --- /dev/null +++ b/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. diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..36e75c0 --- /dev/null +++ b/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 +# . +.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 diff --git a/criticals.c b/criticals.c new file mode 100644 index 0000000..e3ad995 --- /dev/null +++ b/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 + * . + */ +#include +#include +#include + +#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]; +} diff --git a/criticals.h b/criticals.h new file mode 100644 index 0000000..2aa2417 --- /dev/null +++ b/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 + * . + */ +#include +#include + +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); diff --git a/fractal.c b/fractal.c new file mode 100644 index 0000000..722788e --- /dev/null +++ b/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 + * . + */ +#include +#include +#include +#include +#include +#include +#include +#include + +#include "util.h" + +int palette_main(int, char *[]); +int julia_main(int, char *[]); + +static void +usage(void) +{ + fprintf(stderr, "Usage: %s [-s seed] [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 */ +} diff --git a/fun-gen.c b/fun-gen.c new file mode 100644 index 0000000..cce32d0 --- /dev/null +++ b/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 + * . + */ +#include +#include +#include +#include +#include +#include +#include + +#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; +} diff --git a/hi-res.sh b/hi-res.sh new file mode 100644 index 0000000..17cb34e --- /dev/null +++ b/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 +# . +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 diff --git a/julia.c b/julia.c new file mode 100644 index 0000000..2fe7d32 --- /dev/null +++ b/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 + * . + */ +#include +#include +#include +#include +#include +#include +#include + +#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); +} diff --git a/palette.c b/palette.c new file mode 100644 index 0000000..a8a6b46 --- /dev/null +++ b/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 + * . + */ +#include +#include +#include +#include +#include + +#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); +} diff --git a/palette.h b/palette.h new file mode 100644 index 0000000..7b813cd --- /dev/null +++ b/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 + * . + */ +#include +#include + +struct color { + uint16_t r, g, b; +}; + +struct color * palette_gen(size_t); +int palette_main(int, char *[]); diff --git a/pnmout.c b/pnmout.c new file mode 100644 index 0000000..7be0db3 --- /dev/null +++ b/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 + * . + */ +#include +#include + +#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); + } + +} diff --git a/pnmout.h b/pnmout.h new file mode 100644 index 0000000..3564ce8 --- /dev/null +++ b/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 + * . + */ +#include + +void pnmout_header(FILE *, unsigned int, unsigned int); +void pnmout_pixel(FILE *, struct color); diff --git a/util.c b/util.c new file mode 100644 index 0000000..5bd7d59 --- /dev/null +++ b/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 + * . + */ +#include + +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 +} diff --git a/util.h b/util.h new file mode 100644 index 0000000..3bb406d --- /dev/null +++ b/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 + * . + */ +#define nelems(_a) (sizeof((_a)) / sizeof((_a)[0])) + +const char * xgetprogname(void); +void xsrand48(long);