Complex number optimization

Stop using the complex standard lib, replaced it with a variables
for the real part and imaginary part of z.

3 files changed, 51 insertions, 24 deletions

diff --git a/graphics.c b/graphics.c
index 5d3f2ac..f48ed8c 100644
--- a/graphics.c
+++ b/graphics.c
@@ -16,7 +16,7 @@
 #define IMAG_LO(state) (state->center.y - state->imag_range / 2)
 #define IMAG_HI(state) (state->center.y + state->imag_range / 2)
 
-#define IN_SET_COLOR 0x000000
+#define IN_SET_COLOR 0x050505
 #define PALETTE_START 0x000022
 #define PALETTE_END 0xd62f2f
 #define SET_DRAW_COLOR(renderer, c) ( \
@@ -24,7 +24,7 @@
 
 static void update(GState *state);
 static void event_handler(GState *state);
-static Color *create_palette(Color start, Color end);
+/* static Color *create_palette(Color start, Color end); */
 /* static void move_center(GState *state, int motion_x, int motion_y); */
 static void recenter(GState *state, int x, int y);
 static void recenter_x(GState *state, int x);
@@ -96,10 +96,13 @@ static void update(GState *state)
         {
             a = map_range((double)x, 0, state->window_w, REAL_LO(state), REAL_HI(state));
             b = map_range((double)y, 0, state->window_h, IMAG_LO(state), IMAG_HI(state));
-            int steps = mandelbrot_in_set(a + b * I);
-            if (steps == -1)
+            int frac_steps = mandelbrot_in_set(a, b);
+            if (frac_steps == -1)
                 continue;
-            color =  state->palette[steps];
+            color =  state->palette[frac_steps];
+            /* color.rgb.r = 100 * frac_steps; */
+            /* color.rgb.g = 100 * frac_steps; */
+            /* color.rgb.b = 100 * frac_steps; */
             SET_DRAW_COLOR(state->renderer, color);
             SDL_RenderDrawPoint(state->renderer, x, y);
         }
@@ -180,7 +183,7 @@ static void event_handler(GState *state)
     }
 }
 
-static Color *create_palette(Color start, Color end)
+ Color *create_palette(Color start, Color end)
 {
     int red_step = abs(end.rgb.r - start.rgb.r) / MAX_ITERATION;
     int green_step = abs(end.rgb.g - start.rgb.g) / MAX_ITERATION;
diff --git a/header.h b/header.h
index cd5a988..741919e 100644
--- a/header.h
+++ b/header.h
@@ -10,8 +10,9 @@
 # define REAL_RANGE 4.0
 # define IMAG_RANGE 4.0
 
-# define MAX_ITERATION 20
-# define ESCAPE_VALUE 2
+# define MAX_ITERATION 35
+# define ESCAPE_RADIUS 2
+# define ESCAPE_RADIUS_SQUARED (ESCAPE_RADIUS * ESCAPE_RADIUS)
 
 typedef unsigned char Byte;
 typedef int ColorHexcode;
@@ -49,13 +50,16 @@ typedef struct
 } GState;
 
 // mandelbrot.c
-int mandelbrot_in_set(double _Complex c);
+int mandelbrot_in_set(double a, double b);
 void mandelbrot_print(void);
+int *mandelbrot_array(Point center, double real_range, double imag_range,
+                      double real_len, double imag_len);
 
 // graphics.c
 GState *graphics_init(void);
 void graphics_quit(GState *state);
 void graphics_run(GState *state);
+Color *create_palette(Color start, Color end);
 
 // helper.c
 double map_range(double x, double src_lo, double src_hi, double dest_lo, double dest_hi);
diff --git a/mandelbrot.c b/mandelbrot.c
index d44b12b..92e7a64 100644
--- a/mandelbrot.c
+++ b/mandelbrot.c
@@ -14,26 +14,53 @@
 #define IN_CHAR '*'
 #define OUT_CHAR ' '
 
-int mandelbrot_in_set(double complex c)
+int mandelbrot_in_set(double ca, double cb)
 {
-    int i;
-    double complex z = 0;
-    for (i = 0; i < MAX_ITERATION; i++)
+    double zr = ca;
+    double zi = cb;
+    double zr_square;
+    double zi_square;
+    for (int n = 0; n < MAX_ITERATION; n++)
     {
-        z = cpow(z, 2) + c;
-        if (cabs(z) > ESCAPE_VALUE)
-            return i;
+        zi_square = zi * zi;
+        zr_square = zr * zr;
+        if (zr_square + zi_square > ESCAPE_RADIUS_SQUARED)
+            return n;
+            /* return n + 5 - clog(ESCAPE_RADIUS) / log(2); */
+        zi = 2 * zr * zi;
+        zr = zr_square - zi_square;
+        zi += cb;
+        zr += ca;
     }
     return -1;
 }
 
+int *mandelbrot_array(Point center, double real_range, double imag_range,
+                      double real_len, double imag_len)
+{
+    double i = center.y - imag_range / 2;
+    double r = center.x - real_range / 2;
+    int *array = malloc(sizeof(int) * imag_len * real_len);
+
+    for (int array_i = 0; array_i < imag_len;  array_i++)
+    {
+        for (int array_j = 0; array_j < real_len; array_j++)
+        {
+            array[array_i * (int)imag_len + array_j] = mandelbrot_in_set(r, i);
+            r += real_range / real_len;
+        }
+        i += imag_range / imag_len;
+    }
+    return array;
+}
+
 void mandelbrot_print(void)
 {
     for (double i = PRINT_IMAG_LO; i < PRINT_IMAG_HI; i += IMAG_AXIS_STEP)
     {
         for (double r = PRINT_REAL_LO; r < PRINT_REAL_HI; r += REAL_AXIS_STEP)
         {
-            if (mandelbrot_in_set(r + i * I) == -1)
+            if (mandelbrot_in_set(r, i) == -1)
                 putchar(IN_CHAR);
             else
                 putchar(OUT_CHAR);
@@ -42,10 +69,3 @@ void mandelbrot_print(void)
         putchar('\n');
     }
 }
-
-/*
-#define SQUARE(x) (pow((x), 2))
-#define SQUARE_CML(z) (SQUARE(creal(z)) - SQUARE(cimag(z)) + 2 * creal(z) * cimag(z) * I)
-
-double magnitude(double complex z) { return sqrt(SQUARE(creal(z)) + SQUARE(cimag(z))); }
-*/