increased precision from double complex to long double complex

4 files changed, 38 insertions, 37 deletions

diff --git a/src/fractals.h b/src/fractals.h
index dc6783e..5506745 100644
--- a/src/fractals.h
+++ b/src/fractals.h
@@ -5,10 +5,10 @@
 #include <stdint.h>
 #include "grids.h"
 
-size_t mandelbrot(const double complex z0, const size_t max_iterations);
+size_t mandelbrot(const long double complex z0, const size_t max_iterations);
 void mandelbrot_grid(grid_t* grid,  const size_t max_iterations);
 
-size_t multibrot(const double complex z0, const size_t max_iterations, const double d);
+size_t multibrot(const long double 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 double R, const double complex z0, const double complex c, const size_t max_iterations);
+size_t julia(const double R, const long double complex z0, const long double complex c, const size_t max_iterations);
diff --git a/src/grids.c b/src/grids.c
index 89347db..eb93352 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, double complex lower_left, double complex upper_right){
+grid_t* create_grid(const size_t x, const size_t y, long double complex lower_left, long double complex upper_right){
     if(x <= 0 || y <= 0) return NULL;
 
     const size_t size = x * y;
@@ -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
  */
-double complex grid_to_complex(const grid_t* grid, const size_t index) {
+long double 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 double x_min = creal(grid->lower_left);
-    const double x_max = creal(grid->upper_right);
-    const double y_min = cimag(grid->lower_left);
-    const double y_max = cimag(grid->upper_right);
+    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 double x_step = (x_max - x_min) / (double)x_res;
-    const double y_step = (y_max - y_min) / (double)y_res;
+    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 size_t x_index = index % x_res;
     const size_t y_index = index / y_res;
 
-    const double x = x_min + x_index * x_step;
-    const double y = y_min + y_index * y_step;
+    const long double x = x_min + x_index * x_step;
+    const long double y = y_min + y_index * y_step;
 
     return x + y * I;
 }
@@ -128,11 +128,11 @@ double complex grid_to_complex(const grid_t* grid, const size_t index) {
  */
 void zoom_grid(grid_t* grid, const double magnification){
     set_grid(grid, 0);
-    const double complex upper_right = grid->upper_right;
-    const double complex lower_left = grid->lower_left;
+    const long double complex upper_right = grid->upper_right;
+    const long double complex lower_left = grid->lower_left;
 
-    const double complex center = (lower_left + upper_right) / 2.0;
-    const double complex offset = (upper_right - lower_left) / magnification;
+    const long double complex center = (lower_left + upper_right) / 2.0;
+    const long double complex offset = (upper_right - lower_left) / magnification;
 
     grid->lower_left = center - offset;
     grid->upper_right = center + offset;
@@ -163,8 +163,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(double complex), 1, file) != 1 ||
-       fwrite(&grid->upper_right, sizeof(double complex), 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){
         return GRID_WRITE_ERROR;
     }
 
@@ -270,14 +270,14 @@ grid_t* read_grid(FILE* restrict file){
         return NULL;
     }
 
-    double complex lower_left = 0;
-    double complex upper_right = 0; 
-    read_count = fread(&lower_left, sizeof(double complex), 1, file);
+    long double complex lower_left = 0;
+    long double complex upper_right = 0; 
+    read_count = fread(&lower_left, sizeof(long double complex), 1, file);
     if(read_count != 1){
         perror("Error reading file\n");
         return NULL;
     }
-    read_count = fread(&upper_right, sizeof(double complex), 1, file);
+    read_count = fread(&upper_right, sizeof(long double 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 bba6f68..1135e66 100644
--- a/src/grids.h
+++ b/src/grids.h
@@ -15,19 +15,19 @@ typedef struct {
     size_t x;
     size_t y;
     size_t size;
-    double complex lower_left;
-    double complex upper_right;
+    long double complex lower_left;
+    long double complex upper_right;
     size_t* data;
 } grid_t;
 
-grid_t* create_grid(const size_t x, const size_t y, double complex lower_left, double complex upper_right);
+grid_t* create_grid(const size_t x, const size_t y, long double complex lower_left, long double 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);
 
-double complex grid_to_complex(const grid_t* grid, const size_t index);
+long double complex grid_to_complex(const grid_t* grid, const size_t index);
 void zoom_grid(grid_t* grid, const double magnification);
 
 void print_grid_info(const grid_t* grid);
diff --git a/src/serial-fractals.c b/src/serial-fractals.c
index 9f4013d..bd4bc02 100644
--- a/src/serial-fractals.c
+++ b/src/serial-fractals.c
@@ -1,3 +1,4 @@
+#include <complex.h>
 #include "fractals.h"
 #include "grids.h"
 /*
@@ -5,11 +6,11 @@
  * 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 double complex z0, const size_t max_iterations) {
-  double complex z = z0;
+size_t mandelbrot(const long double complex z0, const size_t max_iterations) {
+  long double complex z = z0;
   size_t iteration = 0;
 
-  while (cabs(z) <= 2 && iteration < max_iterations) {
+  while (cabsl(z) <= 2 && iteration < max_iterations) {
     z = z * z + z0;
     iteration++;
   }
@@ -33,10 +34,10 @@ void mandelbrot_grid(grid_t* grid, const size_t max_iterations){
  * 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 double complex z0, const size_t max_iterations, const double d){
-    double complex z = z0;
+size_t multibrot(const long double complex z0, const size_t max_iterations, const double d){
+    long double complex z = z0;
     size_t iteration = 0;
-    while(cabs(z) <= 2 && iteration < max_iterations){
+    while(cabsl(z) <= 2 && iteration < max_iterations){
         z = cpowl(z, d) + z0;
         iteration++;
     }
@@ -59,13 +60,13 @@ void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d){
  * Computes ????? for a julia set
  * implementation of https://en.wikipedia.org/wiki/Julia_set#Pseudocode
  */
-size_t julia(const double R, const double complex z0, const double complex c, const size_t max_iterations){
+size_t julia(const double R, const long double complex z0, const long double complex c, const size_t max_iterations){
     //FIXME: I'm notsure if this is currently implemented correctly
-    if(R*R - R >= cabs(z0)) return 0;
-    double complex z = z0;
+    if(R*R - R >= cabsl(z0)) return 0;
+    long double complex z = z0;
 
     size_t iteration = 0;
-    while(cabs(z) < R && iteration < max_iterations){
+    while(cabsl(z) < R && iteration < max_iterations){
         z = z * z + c;
         iteration++;
     }