Merge branch 'main' of github.com:jpappel/complex-fractals

7 files changed, 390 insertions, 90 deletions

diff --git a/src/fractals.c b/src/fractals.c
index 8615306..269bfaf 100644
--- a/src/fractals.c
+++ b/src/fractals.c
@@ -2,28 +2,70 @@
 #include <unistd.h>
 #include <stdlib.h>
 #include <sys/ioctl.h>
+#include <getopt.h>
 
 #include "grids.h"
+#include "precision.h"
 #include "fractals.h"
 
+
+void print_usage(FILE* file, char* program_name){
+    fprintf(file, "Usage: %s [-v] [-i iterations] [-x x_res] [-y y_res] [-z magnification] [-l lower_left] [-u upper_right]\n", program_name);
+}
+
+void print_help(){
+    //TODO: add help
+}
+
+void print_info(){
+    #ifdef EXTENDED_PRECISION
+    printf("Compiled with long double float precision\n");
+    #endif
+    #ifndef EXTENDED_PRECISION
+    printf("Compiled with double float precision\n");
+    #endif
+}
+
 int main(const int argc, char *argv[]) {
     struct winsize w;
     ioctl(STDOUT_FILENO, TIOCGWINSZ, &w);
 
     //default values
-    size_t iterations = 1000;
+    size_t iterations = 100;
     size_t x_res = w.ws_col;
     size_t y_res = w.ws_row;
-    double re_lower_left = -2;
-    double im_lower_left = -2;
-    double re_upper_right = 2;
-    double im_upper_right = 2;
-    double magnification = 1;
+    CBASE re_lower_left = -2;
+    CBASE im_lower_left = -2;
+    CBASE re_upper_right = 2;
+    CBASE im_upper_right = 2;
+    CBASE magnification = 1;
     bool verbose = false;
+    bool print = false;
+    //TODO: allocate adequate size buffer
+    bool output_to_file = false;
+    char* output_filename = "fractal.grid";
+
+    // TODO: have output format option
+    // TODO: have output file
+
+
+    static struct option long_options[] = {
+        {"iterations", required_argument, NULL, 'i'},
+        {"x-res", required_argument, NULL, 'x'},
+        {"y-res", required_argument, NULL, 'y'},
+        {"lower-left", required_argument, NULL, 'l'},
+        {"upper-right", required_argument, NULL, 'u'},
+        {"magnification", required_argument, NULL, 'z'},
+        {"output", required_argument, NULL, 'o'},
+        {"print", no_argument, NULL, 'p'},
+        {"verbose", no_argument, NULL, 'v'},
+        {"help", no_argument, NULL, 'h'},
+        {0, 0, 0, 0} // Termination element
+    };
 
     //parse command line arguments
     int opt;
-    while((opt = getopt(argc, argv, "vhi:x:y:l:u:z:")) != -1){
+    while((opt = getopt_long(argc, argv, "i:x:y:l:u:z:o:pvh", long_options, NULL)) != -1){
         switch(opt){
             case 'i':
                 iterations = strtoull(optarg, NULL, 10);
@@ -35,15 +77,22 @@ int main(const int argc, char *argv[]) {
                 y_res = strtoull(optarg, NULL, 10);
                 break;
             case 'l':
-                sscanf(optarg, "%lf+%lfi", &re_lower_left, &im_lower_left);
+                sscanf(optarg, CFORMAT"+"CFORMAT"i", &re_lower_left, &im_lower_left);
                 break;
             case 'u':
-                sscanf(optarg, "%lf+%lfi", &re_upper_right, &im_upper_right);
+                sscanf(optarg, CFORMAT"+"CFORMAT"i", &re_upper_right, &im_upper_right);
+                break;
+            case 'o':
+                //TODO: check if within length
+                //TODO: 
+                break;
+            case 'p':
+                print = true;
                 break;
             case 'z':
-                sscanf(optarg, "%lf", &magnification);
+                sscanf(optarg, CFORMAT, &magnification);
                 if(magnification <= 0){
-                    fprintf(stderr, "Invalid magnification %f, exitting\n", magnification);
+                    fprintf(stderr, "Invalid magnification "CFORMAT", exitting\n", magnification);
                     return 1;
                 }
                 break;
@@ -51,41 +100,67 @@ int main(const int argc, char *argv[]) {
                 verbose = true;
                 break;
             case 'h':
-                fprintf(stderr, "Usage: %s [-v] [-i iterations] [-x x_res] [-y y_res] [-z magnification] [-l lower_left] [-u upper_right]\n", argv[0]);
+                print_usage(stdout, argv[0]);
                 return 0;
             default:
-                fprintf(stderr, "Usage: %s [-v] [-i iterations] [-x x_res] [-y y_res] [-z magnification] [-l lower_left] [-u upper_right]\n", argv[0]);
+                print_usage(stderr, argv[0]);
                 return 1;
         }
     }
 
-    const double complex lower_left = re_lower_left + im_lower_left * I;
-    const double complex upper_right = re_upper_right + im_upper_right * I;
+    const CBASE complex lower_left = re_lower_left + im_lower_left * I;
+    const CBASE complex upper_right = re_upper_right + im_upper_right * I;
 
     grid_t* grid = create_grid(x_res, y_res, lower_left, upper_right);
     if(!grid) return 1;
 
 
     if(magnification != 1){
-        if(verbose) printf("Magnification: %f\n", magnification);
+        if(verbose) printf("Magnification: "CFORMAT"\n", magnification);
         zoom_grid(grid, magnification);
     }
 
+    // params for different fractals
+    const double degree = 3.0;
+    const CBASE complex constant = 0.285L + 0.01L*I;
+    // const CBASE complex constant = -0.835L -0.321L* I;
+    const double radius = 100;
 
-    //const long double complex c = 0.285L + 0.01L*I;
-    // const long double complex c = -0.835L -0.321L* I;
-    // const double R = 100;
-    // julia_grid(grid, iterations, c, R);
-    mandelbrot_grid(grid, iterations);
-    //multibrot_grid(grid, iterations, 3);
+    enum fractal f = BURNING_SHIP;
+    switch(f){
+        case MANDELBROT:
+            mandelbrot_grid(grid, iterations);
+            break;
+        case TRICORN:
+            tricorn_grid(grid, iterations);
+            break;
+        case MULTIBROT:
+            multibrot_grid(grid, iterations, degree);
+            break;
+        case MULTICORN:
+            multicorn_grid(grid, iterations, degree);
+            break;
+        case BURNING_SHIP:
+            burning_ship_grid(grid, iterations);
+            break;
+        case JULIA:
+            julia_grid(grid, iterations, constant, radius);
+            break;
+        default:
+            //TODO: update fractal type
+            fprintf(stderr, "Unrecognized fractal\n");
+            return 1;
+    }
 
     if(verbose)print_grid_info(grid);
-    print_grid(grid, iterations);
+    if(print)print_grid(grid, iterations);
 
-    // //write grid to file
-    // FILE* write_file = fopen("test.grid", "wb");
-    // write_grid(write_file , grid);
-    // fclose(write_file);
+    //write grid to file
+    if(output_to_file){
+        FILE* write_file = fopen("test.grid", "wb");
+        write_grid(write_file , grid);
+        fclose(write_file);
+    }
     //
     // //attempt to read grid from file
     // FILE* read_file = fopen("test2.grid", "rb");
@@ -93,4 +168,6 @@ int main(const int argc, char *argv[]) {
     // fclose(read_file);
     //
     // printf("Grids are %s equal\n", grid_equal(grid, grid2) ? "exactly" :"not exactly");
+
+    return 0;
 }
diff --git a/src/fractals.h b/src/fractals.h
index d27eb00..2b296eb 100644
--- a/src/fractals.h
+++ b/src/fractals.h
@@ -4,12 +4,32 @@
 #include <stddef.h>
 #include <stdint.h>
 #include "grids.h"
+#include "precision.h"
 
-size_t mandelbrot(const long double complex z0, const size_t max_iterations);
-void mandelbrot_grid(grid_t* grid,  const size_t max_iterations);
+enum fractal {
+    MANDELBROT, // IMPLEMENTED IN SERIAL SHARED
+    TRICORN, // IMPLEMENTED IN SERIAL SHARED
+    MULTIBROT, // IMPLEMENTED IN SERIAL SHARED
+    MULTICORN, // IMPLEMENTED IN SERIAL SHARED
+    BURNING_SHIP, // IMPLEMENTED IN SERIAL SHARED
+    //NEWTON,  // MIGHT NEVER BE IMPLEMENTED, REQUIRES SPECIAL COLORING
+    JULIA // IMPLEMENTED IN SERIAL SHARED
+};
 
-size_t multibrot(const long double complex z0, const size_t max_iterations, const double d);
+size_t mandelbrot(const CBASE complex z0, const size_t max_iterations);
+void mandelbrot_grid(grid_t* grid, const size_t max_iterations);
+
+size_t tricorn(const CBASE complex z0, const size_t max_iterations);
+void tricorn_grid(grid_t* grid, const size_t max_iterations);
+
+size_t burning_ship(const CBASE complex z0, const size_t max_iterations);
+void burning_ship_grid(grid_t* grid, const size_t max_iterations);
+
+size_t multibrot(const CBASE complex z0, const size_t max_iterations, const double d);
 void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d);
 
-size_t julia(const long double complex z0, const long double complex c, const size_t max_iterations, const double R);
-void julia_grid(grid_t* grid, const size_t max_iterations, const long double complex c, const double R);
+size_t multicorn(const CBASE complex z0, const size_t max_iterations, const double d);
+void multicorn_grid(grid_t* grid, const size_t max_iterations, const double d);
+
+size_t julia(const CBASE complex z0, const CBASE complex c, const size_t max_iterations, const double R);
+void julia_grid(grid_t* grid, const size_t max_iterations, const CBASE complex c, const double R);
diff --git a/src/grids.c b/src/grids.c
index 54acecc..09c5c49 100644
--- a/src/grids.c
+++ b/src/grids.c
@@ -7,7 +7,7 @@
 /*
  * Creates a grid for storing the results of the escape algorithm
  */
-grid_t* create_grid(const size_t x, const size_t y, long double complex lower_left, long double complex upper_right){
+grid_t* create_grid(const size_t x, const size_t y, CBASE complex lower_left, CBASE complex upper_right){
     if(x <= 0 || y <= 0) return NULL;
 
     const size_t size = x * y;
@@ -83,7 +83,7 @@ bool grid_equal(const grid_t* grid1, const grid_t* grid2){
 /*
  * Checks if two grids have a given maximum difference
  */
-bool grid_allclose(const grid_t *grid1, const grid_t *grid2, const size_t max_error){
+bool grid_allclose(const grid_t* restrict grid1, const grid_t* restrict grid2, const size_t max_error){
     if(grid1->x != grid2->x || grid1->y != grid2->y ||
             grid1->lower_left != grid2->lower_left || grid1->upper_right != grid2->upper_right){
         return false;
@@ -101,22 +101,22 @@ bool grid_allclose(const grid_t *grid1, const grid_t *grid2, const size_t max_er
 /*
  * Converts a grid point into the corresponding complex number
  */
-long double complex grid_to_complex(const grid_t* grid, const size_t index) {
+CBASE complex grid_to_complex(const grid_t* grid, const size_t index) {
     const size_t x_res = grid->x;
     const size_t y_res = grid->y;
-    const long double x_min = creal(grid->lower_left);
-    const long double x_max = creal(grid->upper_right);
-    const long double y_min = cimag(grid->lower_left);
-    const long double y_max = cimag(grid->upper_right);
+    const CBASE x_min = CREAL(grid->lower_left);
+    const CBASE x_max = CREAL(grid->upper_right);
+    const CBASE y_min = CIMAG(grid->lower_left);
+    const CBASE y_max = CIMAG(grid->upper_right);
 
-    const long double x_step = (x_max - x_min) / (double)x_res;
-    const long double y_step = (y_max - y_min) / (double)y_res;
+    const CBASE x_step = (x_max - x_min) / (double)x_res;
+    const CBASE y_step = (y_max - y_min) / (double)y_res;
 
     const size_t x_index = index % x_res;
     const size_t y_index = index / y_res;
 
-    const long double x = x_min + x_index * x_step;
-    const long double y = y_min + y_index * y_step;
+    const CBASE x = x_min + x_index * x_step;
+    const CBASE y = y_min + y_index * y_step;
 
     return x + y * I;
 }
@@ -126,13 +126,13 @@ long double complex grid_to_complex(const grid_t* grid, const size_t index) {
  *
  * Resets all grid values to 0
  */
-void zoom_grid(grid_t* grid, const double magnification){
+void zoom_grid(grid_t* restrict grid, const CBASE magnification){
     set_grid(grid, 0);
-    const long double complex upper_right = grid->upper_right;
-    const long double complex lower_left = grid->lower_left;
+    const CBASE complex upper_right = grid->upper_right;
+    const CBASE complex lower_left = grid->lower_left;
 
-    const long double complex center = (lower_left + upper_right) / 2.0;
-    const long double complex offset = (upper_right - lower_left) / magnification;
+    const CBASE complex center = (lower_left + upper_right) / 2.0;
+    const CBASE complex offset = (upper_right - lower_left) / magnification;
 
     grid->lower_left = center - offset;
     grid->upper_right = center + offset;
@@ -166,8 +166,8 @@ int write_grid(FILE* restrict file, const grid_t *grid){
 
     if(fwrite(&grid->x, sizeof(size_t), 1, file) != 1 ||
        fwrite(&grid->y, sizeof(size_t), 1, file) != 1 ||
-       fwrite(&grid->lower_left, sizeof(long double complex), 1, file) != 1 ||
-       fwrite(&grid->upper_right, sizeof(long double complex), 1, file) != 1){
+       fwrite(&grid->lower_left, sizeof(CBASE complex), 1, file) != 1 ||
+       fwrite(&grid->upper_right, sizeof(CBASE complex), 1, file) != 1){
         return GRID_WRITE_ERROR;
     }
 
@@ -190,8 +190,8 @@ void print_grid_info(const grid_t* grid){
     printf("x\t%zu\n", grid->x);
     printf("y\t%zu\n", grid->y);
     printf("size\t%zu\n", grid->size);
-    printf("lower_left\t%f + %fI\n", creal(grid->lower_left), cimag(grid->lower_left));
-    printf("upper_right\t%f + %fI\n", creal(grid->upper_right), cimag(grid->upper_right));
+    printf("lower_left\t"CFORMAT"+ "CFORMAT"I\n", CREAL(grid->lower_left), CIMAG(grid->lower_left));
+    printf("upper_right\t"CFORMAT"+ "CFORMAT"I\n", CREAL(grid->upper_right), CIMAG(grid->upper_right));
 
     printf("Data is %s NULL\n", grid->data ? "not" : "");
 }
@@ -204,6 +204,8 @@ void print_grid(const grid_t* grid, const size_t iterations){
     const size_t x_res = grid->x;
     const size_t* data = grid->data;
 
+    //TODO: set values in output buffer rather than multiple printf calls
+    //      the buffer needs to be larger to hold newlines
     // char* output_buffer = malloc(size);
     // if(!output_buffer){
     //     fprintf(stderr, "Failed to allocate output buffer for %zu points\n");
@@ -215,7 +217,6 @@ void print_grid(const grid_t* grid, const size_t iterations){
     size_t last_bin = iterations - bin_width;
     char point;
 
-    //TODO: set values in output buffer rather than multiple printf calls
     for(size_t i = 0; i < size; i++){
         const size_t value = data[i];
         if(value == iterations){
@@ -270,14 +271,14 @@ grid_t* read_grid(FILE* restrict file){
         return NULL;
     }
 
-    long double complex lower_left = 0;
-    long double complex upper_right = 0; 
-    read_count = fread(&lower_left, sizeof(long double complex), 1, file);
+    CBASE complex lower_left = 0;
+    CBASE complex upper_right = 0; 
+    read_count = fread(&lower_left, sizeof(CBASE complex), 1, file);
     if(read_count != 1){
         perror("Error reading file\n");
         return NULL;
     }
-    read_count = fread(&upper_right, sizeof(long double complex), 1, file);
+    read_count = fread(&upper_right, sizeof(CBASE complex), 1, file);
     if(read_count != 1){
         perror("Error reading file\n");
         return NULL;
diff --git a/src/grids.h b/src/grids.h
index 1135e66..8448fa6 100644
--- a/src/grids.h
+++ b/src/grids.h
@@ -4,6 +4,7 @@
 #include <stdio.h>
 #include <stdbool.h>
 #include <complex.h>
+#include "precision.h"
 
 //grid write errors
 #define GRID_NO_DATA 1
@@ -15,20 +16,20 @@ typedef struct {
     size_t x;
     size_t y;
     size_t size;
-    long double complex lower_left;
-    long double complex upper_right;
+    CBASE complex lower_left;
+    CBASE complex upper_right;
     size_t* data;
 } grid_t;
 
-grid_t* create_grid(const size_t x, const size_t y, long double complex lower_left, long double complex upper_right);
+grid_t* create_grid(const size_t x, const size_t y, CBASE complex lower_left, CBASE complex upper_right);
 void set_grid(grid_t* grid, const size_t val);
 grid_t* copy_grid(const grid_t* grid);
 void free_grid(grid_t* grid);
 bool grid_equal(const grid_t* grid1, const grid_t* grid2);
 bool grid_allclose(const grid_t* grid1, const grid_t* grid2, const size_t max_error);
 
-long double complex grid_to_complex(const grid_t* grid, const size_t index);
-void zoom_grid(grid_t* grid, const double magnification);
+CBASE complex grid_to_complex(const grid_t* grid, const size_t index);
+void zoom_grid(grid_t* grid, const CBASE magnification);
 
 void print_grid_info(const grid_t* grid);
 void print_grid(const grid_t* grid, const size_t iterations);
diff --git a/src/precision.h b/src/precision.h
new file mode 100644
index 0000000..4bca354
--- /dev/null
+++ b/src/precision.h
@@ -0,0 +1,31 @@
+/*
+ * Header that allows conditional compilation into double floating point precisions ot extended double floating point precision
+ *
+ * CUDA does not support long double so cuda-fractals may fail to compile or have unexpected behavior if linked with code that was compiled with EXTENDED_PRECISION
+ */
+#pragma once
+
+#ifdef EXTENDED_PRECISION
+
+#define CBASE long double
+#define CREAL creall
+#define CIMAG cimagl
+#define CPOW cpowl
+#define CONJ conjl
+#define CABS cabsl
+#define RABS fabsl
+#define CFORMAT "%Lf"
+
+#endif
+#ifndef EXTENDED_PRECISION
+
+#define CBASE double
+#define CREAL creal
+#define CIMAG cimag
+#define CPOW cpow
+#define CONJ conj
+#define CABS cabs
+#define RABS fabs
+#define CFORMAT "%lf"
+
+#endif
diff --git a/src/serial-fractals.c b/src/serial-fractals.c
index fe420d6..4e081c3 100644
--- a/src/serial-fractals.c
+++ b/src/serial-fractals.c
@@ -1,16 +1,18 @@
 #include "fractals.h"
+#include <math.h>
+#include "precision.h"
 #include "grids.h"
 
 /*
- * Computes the number of iterations it takes for a point z0 to diverge
+ * Computes the number of iterations it takes for a point z0 to become unbounded
  * if the return value is equal to max_iterations, the point lies within the mandelbrot set
  * This is an implementation the escape algorithm
  */
-size_t mandelbrot(const long double complex z0, const size_t max_iterations) {
-  long double complex z = z0;
+size_t mandelbrot(const CBASE complex z0, const size_t max_iterations) {
+  CBASE complex z = z0;
   size_t iteration = 0;
 
-  while (cabsl(z) <= 2 && iteration < max_iterations) {
+  while (CABS(z) <= 2 && iteration < max_iterations) {
     z = z * z + z0;
     iteration++;
   }
@@ -30,15 +32,71 @@ void mandelbrot_grid(grid_t* grid, const size_t max_iterations){
 }
 
 /*
- * Computes the number of iterations it takes for a point z0 to diverge
+ * Computes the number of iterations it takes for a point z0 to become unbounded
+ * if the return value is equal to max_iterations, the point lies within the tricorn set
+ * This is nearly identical to mandelbrot, except for the complex conjugate
+ */
+size_t tricorn(const CBASE complex z0, const size_t max_iterations){
+    CBASE complex z = z0;
+    size_t iteration = 0;
+    while(CABS(z) <= 2 && iteration < max_iterations){
+        z = CONJ(z * z) + z0;
+        iteration++;
+    }
+    return iteration;
+}
+
+/*
+ * Fills a grid with tricorn values
+ */
+void tricorn_grid(grid_t* grid, const size_t max_iterations){
+    const size_t size = grid->size;
+    size_t* data = grid->data;
+
+    for(size_t i = 0; i < size; i++){
+        data[i] = tricorn(grid_to_complex(grid, i), max_iterations);
+    }
+}
+
+/*
+ * Computes the number of iterations it takes for a point z0 to become unbounded
+ * if the return value is equal to max_iterations, the point lies within the burningship set (oh no! I hope they have fire safety gear)
+ */
+size_t burning_ship(const CBASE complex z0, const size_t max_iterations) {
+  CBASE complex z = z0;
+  CBASE complex z_mod;
+  size_t iteration = 0;
+
+  while (CABS(z) <= 2 && iteration < max_iterations) {
+    z_mod = RABS(CREAL(z)) + RABS(CIMAG(z))*I;
+    z = z_mod * z_mod + z0;
+    iteration++;
+  }
+  return iteration;
+}
+
+/*
+ * Fills a grid with burning_ship values
+ */
+void burning_ship_grid(grid_t* grid, const size_t max_iterations){
+    const size_t size = grid->size;
+    size_t* data = grid->data;
+
+    for(size_t i = 0; i < size; i++){
+        data[i] = burning_ship(grid_to_complex(grid, i), max_iterations);
+    }
+}
+
+/*
+ * Computes the number of iterations it takes for a point z0 to become unbounded
  * if the return value is equal to max_iterations, the point lies within the multibrot set
  * This is implementation closely matches mandelbrot, but uses cpow which might degrade performance.
  */
-size_t multibrot(const long double complex z0, const size_t max_iterations, const double d){
-    long double complex z = z0;
+size_t multibrot(const CBASE complex z0, const size_t max_iterations, const double d){
+    CBASE complex z = z0;
     size_t iteration = 0;
-    while(cabsl(z) <= 2 && iteration < max_iterations){
-        z = cpowl(z, d) + z0;
+    while(CABS(z) <= 2 && iteration < max_iterations){
+        z = CPOW(z, d) + z0;
         iteration++;
     }
     return iteration;
@@ -57,23 +115,49 @@ void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d){
 }
 
 /*
+ * Computes the number ofiterations it takes for a point z0 to become unbounded
+ * if the return value is equal to max_iterations, the point lies within the multicorn set
+ * This function is to tricorn as multibrot is to mandelbrot
+ */
+size_t multicorn(const CBASE complex z0, const size_t max_iterations, const double d){
+    CBASE complex z = z0;
+    size_t iteration = 0;
+    while(CABS(z) <= 2 && iteration < max_iterations){
+        z = CONJ(CPOW(z, d)) + z0;
+        iteration++;
+    }
+    return iteration;
+}
+
+/*
+ * Fills a grid with multicorn values
+ */
+void multicorn_grid(grid_t* grid, const size_t max_iterations, const double d){
+    const size_t size = grid->size;
+    size_t* data = grid->data;
+    for(size_t i = 0; i < size; i ++){
+        data[i] = multicorn(grid_to_complex(grid, i), max_iterations, d);
+    }
+}
+
+/*
  * Computes ????? for a julia set
  * implementation of https://en.wikipedia.org/wiki/Julia_set#Pseudocode
  *
  * This behaves weirdly, needs a very small number of iterations to be visibile
  */
-size_t julia(const long double complex z0, const long double complex c, const size_t max_iterations, const double R){
-    long double complex z = z0;
+size_t julia(const CBASE complex z0, const CBASE complex c, const size_t max_iterations, const double R){
+    double complex z = z0;
 
     size_t iteration = 0;
-    while(cabsl(z) < R && iteration < max_iterations){
+    while(CABS(z) < R && iteration < max_iterations){
         z = z * z + c;
         iteration++;
     }
     return iteration;
 }
 
-void julia_grid(grid_t* grid, const size_t max_iterations, const long double complex c, const double R){
+void julia_grid(grid_t* grid, const size_t max_iterations, const CBASE complex c, const double R){
     const size_t size = grid->size;
     size_t* data = grid->data;
     for(size_t i = 0; i <size; i++){
diff --git a/src/shared-fractals.c b/src/shared-fractals.c
index 49a034e..92c8c9a 100644
--- a/src/shared-fractals.c
+++ b/src/shared-fractals.c
@@ -1,17 +1,19 @@
 #include "fractals.h"
 
 #include <omp.h>
+#include <math.h>
 #include "fractals.h"
+#include "precision.h"
 
 /*
  * Computes the number of iterations it takes for a point z0 to diverge
  * if the return value is equal to max_iterations, the point lies within the mandelbrot set
  * This should be identical to the version found in serial-fractals.c
  */
-size_t mandelbrot(const long double complex z0, const size_t max_iterations){
-    long double complex z = z0;
+size_t mandelbrot(const CBASE complex z0, const size_t max_iterations){
+    CBASE complex z = z0;
     size_t iteration = 0;
-    while(cabsl(z) <= 2 && iteration < max_iterations){
+    while(CABS(z) <= 2 && iteration < max_iterations){
         z = z*z + z0;
         iteration++;
     }
@@ -22,7 +24,7 @@ size_t mandelbrot(const long double complex z0, const size_t max_iterations){
 /*
  * Fills a grid with mandelbrot values
  */
-void mandelbrot_grid(grid_t* grid, const size_t max_iterations){
+void mandelbrot_grid(grid_t* restrict grid, const size_t max_iterations){
     const size_t size = grid->size;
     size_t* data = grid->data;
 
@@ -33,25 +35,82 @@ void mandelbrot_grid(grid_t* grid, const size_t max_iterations){
 }
 
 /*
+ * Computes the number of iterations it takes for a point z0 to become unbounded
+ * if the return value is equal to max_iterations, the point lies within the tricorn set
+ * This is nearly identical to mandelbrot, except for the complex conjugate
+ */
+size_t tricorn(const CBASE complex z0, const size_t max_iterations){
+    CBASE complex z = z0;
+    size_t iteration = 0;
+    while(CABS(z) <= 2 && iteration < max_iterations){
+        z = CONJ(z * z) + z0;
+        iteration++;
+    }
+    return iteration;
+}
+
+/*
+ * Fills a grid with tricorn values
+ */
+void tricorn_grid(grid_t* grid, const size_t max_iterations){
+    const size_t size = grid->size;
+    size_t* data = grid->data;
+
+    #pragma omp parallel for default(none) shared(data, size, grid, max_iterations) schedule(dynamic)
+    for(size_t i = 0; i < size; i++){
+        data[i] = tricorn(grid_to_complex(grid, i), max_iterations);
+    }
+}
+
+/*
+ * Computes the number of iterations it takes for a point z0 to become unbounded
+ * if the return value is equal to max_iterations, the point lies within the burningship set (oh no! I hope they have fire safety gear)
+ */
+size_t burning_ship(const CBASE complex z0, const size_t max_iterations) {
+  CBASE complex z = z0;
+  CBASE complex z_mod;
+  size_t iteration = 0;
+
+  while (CABS(z) <= 2 && iteration < max_iterations) {
+    z_mod = RABS(CREAL(z)) + RABS(CIMAG(z))*I;
+    z = z_mod * z_mod + z0;
+    iteration++;
+  }
+  return iteration;
+}
+
+/*
+ * Fills a grid with burning_ship values
+ */
+void burning_ship_grid(grid_t* grid, const size_t max_iterations){
+    const size_t size = grid->size;
+    size_t* data = grid->data;
+
+    #pragma omp parallel for default(none) shared(data, size, grid, max_iterations) schedule(dynamic)
+    for(size_t i = 0; i < size; i++){
+        data[i] = burning_ship(grid_to_complex(grid, i), max_iterations);
+    }
+}
+
+/*
  * Computes the number of iterations it takes for a point z0 to diverge
  * if the return value is equal to max_iterations, the point lies within the multibrot set
  * This should be identical to the version found in serial-fractals.c
  */
-size_t multibrot(const long double complex z0, const size_t max_iterations, const double d){
-    long double complex z = z0;
+size_t multibrot(const CBASE complex z0, const size_t max_iterations, const double d){
+    CBASE complex z = z0;
     size_t iteration = 0;
-    while(cabsl(z) <= 2 && iteration < max_iterations){
-        z = cpowl(z, d) + z0;
+    while(CABS(z) <= 2 && iteration < max_iterations){
+        z = CPOW(z, d) + z0;
         iteration++;
     }
     return iteration;
 }
 
-
 /*
  * Fills a grid with multibrot values
  */
-void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d){
+void multibrot_grid(grid_t* restrict grid, const size_t max_iterations, const double d){
     const size_t size = grid->size;
     size_t* data = grid->data;
 
@@ -62,28 +121,55 @@ void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d){
 }
 
 /*
+ * Computes the number ofiterations it takes for a point z0 to become unbounded
+ * if the return value is equal to max_iterations, the point lies within the multicorn set
+ * This function is to tricorn as multibrot is to mandelbrot
+ */
+size_t multicorn(const CBASE complex z0, const size_t max_iterations, const double d){
+    CBASE complex z = z0;
+    size_t iteration = 0;
+    while(CABS(z) <= 2 && iteration < max_iterations){
+        z = CONJ(CPOW(z, d)) + z0;
+        iteration++;
+    }
+    return iteration;
+}
+
+/*
+ * Fills a grid with multicorn values
+ */
+void multicorn_grid(grid_t* grid, const size_t max_iterations, const double d){
+    const size_t size = grid->size;
+    size_t* data = grid->data;
+    #pragma omp parallel for default(none) shared(data, size, grid, max_iterations, d) schedule(dynamic)
+    for(size_t i = 0; i < size; i ++){
+        data[i] = multicorn(grid_to_complex(grid, i), max_iterations, d);
+    }
+}
+
+/*
  * Computes ????? for a julia set
  * implementation of https://en.wikipedia.org/wiki/Julia_set#Pseudocode
  *
  * This behaves weirdly, needs a very small number of iterations to be visibile
  */
 
-size_t julia(const long double complex z0, const long double complex c, const size_t max_iterations, const double R){
-    long double complex z = z0;
+size_t julia(const CBASE complex z0, const CBASE complex c, const size_t max_iterations, const double R){
+    double complex z = z0;
 
     size_t iteration = 0;
-    while(cabsl(z) < R && iteration < max_iterations){
+    while(CABS(z) < R && iteration < max_iterations){
         z = z * z + c;
         iteration++;
     }
     return iteration;
 }
 
-void julia_grid(grid_t* grid, const size_t max_iterations, const long double complex c, const double R){
+void julia_grid(grid_t* restrict grid, const size_t max_iterations, const CBASE complex c, const double R){
     const size_t size = grid->size;
     size_t* data = grid->data;
     #pragma omp parallel for default(none) shared(data, size, grid, max_iterations, c, R) schedule(dynamic)
-    for(size_t i = 0; i <size; i++){
+    for(size_t i = 0; i < size; i++){
         data[i] = julia(grid_to_complex(grid, i), c, max_iterations, R);
     }
 }