ref: c52df87d07198a3140227a03a9ddddddd418777a
author: sirjofri <[email protected]>
date: Sat Dec 9 13:02:20 EST 2023
adds files
--- /dev/null
+++ b/complex.c
@@ -1,0 +1,27 @@
+#include <u.h>
+#include <libc.h>
+#include "complex.h"
+
+Complex
+cadd(Complex *a, Complex *b)
+{+ Complex r;
+ r.r = a->r + b->r;
+ r.i = a->i + b->i;
+ return r;
+}
+
+Complex
+cmul(Complex *a, Complex *b)
+{+ Complex r;
+ r.r = a->r*b->r - a->i*b->i;
+ r.i = a->r*b->i + a->i*b->r;
+ return r;
+}
+
+Complex
+cpow2(Complex *a)
+{+ return cmul(a, a);
+}
\ No newline at end of file
--- /dev/null
+++ b/complex.h
@@ -1,0 +1,9 @@
+typedef struct Complex Complex;
+struct Complex {+ double r;
+ double i;
+};
+
+Complex cadd(Complex *a, Complex *b);
+Complex cmul(Complex *a, Complex *b);
+Complex cpow2(Complex *a);
--- /dev/null
+++ b/fractal.c
@@ -1,0 +1,96 @@
+#include <u.h>
+#include <libc.h>
+#include <draw.h>
+#include <memdraw.h>
+#include "complex.h"
+#include "gen.h"
+
+void
+usage(void)
+{+ fprint(2, "usage: %s [-j a,b] [-s x,y] [-h]\n", argv0);
+ exits("usage");+}
+
+void
+parsejulia(Complex *c, char *str)
+{+ if (!str)
+ return;
+
+ char *args[2];
+ int n = getfields(str, args, 2, 1, ",");
+ if (n != 2)
+ usage();
+
+ c->r = atof(args[0]);
+ c->i = atof(args[1]);
+}
+
+void
+parsesize(int *sizex, int *sizey, char *str)
+{+ if (!str)
+ return;
+
+ char *args[2];
+ int n = getfields(str, args, 2, 1, "x");
+ if (n != 2)
+ usage();
+
+ *sizex = atoi(args[0]);
+ *sizey = atoi(args[1]);
+}
+
+void
+main(int argc, char **argv)
+{+ int sizex = 256;
+ int sizey = 256;
+ char *s;
+ SetParams params;
+ DrawParams dparams;
+ Memimage *img;
+
+ params.iterations = 10;
+ dparams.scale = 3.;
+ dparams.offsetx = 0.;
+ dparams.offsety = 0.;
+
+ ARGBEGIN {+ case 'j':
+ s = EARGF(usage());
+ parsejulia(¶ms.c, s);
+ params.julia = 1;
+ break;
+ case 's':
+ s = EARGF(usage());
+ parsesize(&sizex, &sizey, s);
+ params.c.r = 0.;
+ params.c.i = 0.;
+ break;
+ case 'z':
+ dparams.scale = atof(EARGF(usage()));
+ break;
+ case 'x':
+ dparams.offsetx = atof(EARGF(usage()));
+ break;
+ case 'y':
+ dparams.offsety = atof(EARGF(usage()));
+ break;
+ case 'i':
+ params.iterations = atoi(EARGF(usage()));
+ break;
+ case 'h':
+ usage();
+ } ARGEND;
+
+ memimageinit();
+ img = generate(sizex, sizey, params, dparams);
+ if (!img)
+ sysfatal("error: %r");+
+ writememimage(1, img);
+ freememimage(img);
+ exits(nil);
+}
--- /dev/null
+++ b/gen.c
@@ -1,0 +1,147 @@
+#include <u.h>
+#include <libc.h>
+#include <draw.h>
+#include <memdraw.h>
+#include "complex.h"
+#include "gen.h"
+
+int is∈(Complex *c);
+double len(Complex *c)
+{+ Complex r;
+ r.r = c->r * c->r;
+ r.r = r.r < 0.00001 && r.r > 0.00001 ? 0.00001 : r.r;
+ r.i = c->i * c->i;
+ r.i = r.i < 0.00001 && r.i > 0.00001 ? 0.00001 : r.i;
+ double ret = r.r + r.i;
+ ret = ret < 3200000. ? ret : 3200000.;
+ return sqrt(ret);
+};
+
+float lerp(float a, float b, float s)
+{+ return s*b + (1.-s)*a;
+}
+
+Memimage*
+generate(int x, int y, SetParams params, DrawParams dparams)
+{+ if (x < 8 || y < 8) {+ werrstr("invalid image size");+ return nil;
+ }
+
+ Rectangle r;
+ r.min.x = r.min.y = 0;
+ r.max.x = x;
+ r.max.y = y;
+ Memimage *img = allocmemimage(r, XRGB32);
+ if (!img) {+ return nil;
+ }
+
+ for (int nx = 0; nx < x; nx++) {+ for (int ny = 0; ny < y; ny++) {+ int valid, red, green, blue;
+ Complex c;
+
+ float uvx = (float)nx/(float)x;
+ float uvy = (float)ny/(float)y;
+ uvx -= .5;
+ uvy -= .5;
+ uvx *= dparams.scale;
+ uvy *= dparams.scale;
+ uvx += dparams.offsetx;
+ uvy += dparams.offsety;
+
+ c.r = uvx;
+ c.i = uvy;
+
+ valid = calc(&c, c, params);
+ if (!valid)
+ goto err;
+
+ if (is∈(&c)) {+ red = 0;
+ green = 0;
+ blue = (int)(len(&c) * 255.);
+ } else {+ float fmix = (float)valid/(float)params.iterations;
+ red = lerp(0., 255., pow(fmix, 0.5));
+ green = lerp(255., 0., pow(fmix, 0.2)) * 0.5;
+ blue = 0;
+ }
+
+ if (red > 255)
+ red = 255;
+ if (green > 255)
+ green = 255;
+ if (blue > 255)
+ blue = 255;
+
+ int arri = (ny*x + nx) * 4;
+ img->data->bdata[arri+0] = blue;
+ img->data->bdata[arri+1] = green;
+ img->data->bdata[arri+2] = red;
+ img->data->bdata[arri+3] = 255;
+ }
+ }
+
+ return img;
+
+err:
+ freememimage(img);
+ return nil;
+}
+
+int
+calc(Complex *out, Complex c, SetParams params)
+{+ Complex t, f;
+
+ if (params.iterations < 1) {+ werrstr("iterations < 1");+ return 0;
+ }
+
+ if (params.julia) {+ /* calculate julia */
+ t = c;
+ f = params.c;
+ } else {+ /* calculate mandelbrot */
+ t = params.c;
+ f = c;
+ }
+
+ int iter = -1;
+
+ for (int i = 0; i < params.iterations; i++) {+ // abort if values grow too large
+ if (t.r > 64. || t.i > 64.)
+ break;
+ if (iter < 0 && t.r >= 2. && t.i >= 2.)
+ iter = i;
+ t = cpow2(&t);
+ t = cadd(&t, &f);
+ }
+
+ *out = t;
+ return iter > 0 ? iter : 1;
+}
+
+int
+∈(Complex c, SetParams params)
+{+ Complex r;
+ int valid = calc(&r, c, params);
+ if (!valid)
+ return -1;
+ return is∈(&r);
+}
+
+int
+is∈(Complex *c)
+{+ return (c->r <= 2. && c->i <= 2.);
+}
--- /dev/null
+++ b/gen.h
@@ -1,0 +1,32 @@
+typedef struct SetParams SetParams;
+struct SetParams {+ int julia;
+ Complex c;
+ int iterations;
+};
+
+typedef struct DrawParams DrawParams;
+struct DrawParams {+ float scale;
+ float offsetx;
+ float offsety;
+};
+
+/* generate
+ * nil → err
+ * * → image
+ */
+Memimage* generate(int x, int y, SetParams params, DrawParams dparams);
+
+/* iterate
+ * 0 → err
+ * >0 → out is valid, value is num iterations
+ */
+int calc(Complex *out, Complex c, SetParams params);
+
+/* is c inside ℂ
+ * -1 → err
+ * 0 → c ∉ ℂ
+ * 1 → c ∈ ℂ
+ */
+int ∈(Complex c, SetParams params);
--- /dev/null
+++ b/mkfile
@@ -1,0 +1,26 @@
+</$objtype/mkfile
+
+TARG=fractal
+OFILES=fractal.$O gen.$O complex.$O
+HFILES=complex.h gen.h
+
+all:V: $O.out
+
+test:V: julia.bit mandel.bit iter.bit size.bit
+
+julia.bit: $O.out
+ $O.out -j 0.28,0.01 > $target
+
+mandel.bit: $O.out
+ $O.out -x -.5 > $target
+
+iter.bit: $O.out
+ $O.out -x -.5 -i 30 > $target
+
+size.bit: $O.out
+ $O.out -x -.5 -i 30 -s 1024x1024 > $target
+
+</sys/src/cmd/mkone
+
+# for profiling. WARNING: this changes the result for some reason!
+#LDFLAGS=-p