From 9ba9c47a952ce6966b333af579bd39c636080fbc Mon Sep 17 00:00:00 2001 From: JP Appel Date: Tue, 23 Apr 2024 15:58:54 -0400 Subject: added preprocessor to change precision at compile time --- src/fractals.c | 25 +++++++++++++------------ src/fractals.h | 9 +++++---- src/grids.c | 51 ++++++++++++++++++++++++++------------------------- src/grids.h | 11 ++++++----- src/precision.h | 29 +++++++++++++++++++++++++++++ src/serial-fractals.c | 23 ++++++++++++----------- src/shared-fractals.c | 29 +++++++++++++++-------------- 7 files changed, 106 insertions(+), 71 deletions(-) create mode 100644 src/precision.h (limited to 'src') diff --git a/src/fractals.c b/src/fractals.c index 8615306..b43ad0e 100644 --- a/src/fractals.c +++ b/src/fractals.c @@ -4,6 +4,7 @@ #include #include "grids.h" +#include "precision.h" #include "fractals.h" int main(const int argc, char *argv[]) { @@ -14,11 +15,11 @@ int main(const int argc, char *argv[]) { size_t iterations = 1000; 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; //parse command line arguments @@ -35,15 +36,15 @@ 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 '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; @@ -59,15 +60,15 @@ int main(const int argc, char *argv[]) { } } - 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); } diff --git a/src/fractals.h b/src/fractals.h index d27eb00..69e5215 100644 --- a/src/fractals.h +++ b/src/fractals.h @@ -4,12 +4,13 @@ #include #include #include "grids.h" +#include "precision.h" -size_t mandelbrot(const long double complex z0, const size_t max_iterations); +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 multibrot(const long double complex z0, const size_t max_iterations, const double d); +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 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..6941d66 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 #include #include +#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..e83122a --- /dev/null +++ b/src/precision.h @@ -0,0 +1,29 @@ +/* + * 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 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 CFORMAT "%lf" + +#endif diff --git a/src/serial-fractals.c b/src/serial-fractals.c index fe420d6..5831f92 100644 --- a/src/serial-fractals.c +++ b/src/serial-fractals.c @@ -1,4 +1,5 @@ #include "fractals.h" +#include "precision.h" #include "grids.h" /* @@ -6,11 +7,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 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++; } @@ -34,11 +35,11 @@ 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 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; @@ -62,18 +63,18 @@ void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d){ * * 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 #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 +23,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; @@ -37,11 +38,11 @@ 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 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; @@ -51,7 +52,7 @@ size_t multibrot(const long double complex z0, const size_t max_iterations, cons /* * 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; @@ -68,22 +69,22 @@ void multibrot_grid(grid_t* grid, const size_t max_iterations, const double d){ * 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