diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/fractals.c | 131 | ||||
| -rw-r--r-- | src/fractals.h | 30 | ||||
| -rw-r--r-- | src/grids.c | 51 | ||||
| -rw-r--r-- | src/grids.h | 11 | ||||
| -rw-r--r-- | src/precision.h | 31 | ||||
| -rw-r--r-- | src/serial-fractals.c | 110 | ||||
| -rw-r--r-- | src/shared-fractals.c | 116 |
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); } } |