Restructured versions and file names - Changed main implementation

This commit is contained in:
Yashar Alan Zenker
2026-07-04 20:13:20 +02:00
parent 9c28454b7e
commit b14cf7a79c
11 changed files with 384 additions and 380 deletions
+9 -6
View File
@@ -5,12 +5,15 @@
#include <unistd.h> // for getopt
#include <getopt.h> // for getopt_long
#include <string.h> // for strcpy
#include "tricorn_V0.h"
// All implementations
#include "tricorn.h" // Final implmementation
#include "tricorn_V1.h"
#include "tricorn_V2.h"
#include "tricorn_V3.h"
#include "tricorn_V4.h"
void save_bmp(const char *filename, unsigned char *img, size_t width, ssize_t height,_Bool color);
unsigned tests(void) {
@@ -207,11 +210,11 @@ int main(int argc, char** argv)
switch (version_to_use)
{
case 0:
// Uses highly optimized SIMD version
tricorn_V0(complexA,width,height,resolution,n,colored,tempArr);
// Uses highly optimized SIMD and threaded version / main implementation
tricorn(complexA,width,height,resolution,n,colored,tempArr);
break;
case 1:
// Uses unoptimized SIMD version
// Uses unoptimised comparison implementation
tricorn_V1(complexA,width,height,resolution,n,colored,tempArr);
break;
case 2:
@@ -219,11 +222,11 @@ int main(int argc, char** argv)
tricorn_V2(complexA,width,height,resolution,n,colored,tempArr);
break;
case 3:
// Uses unoptimised comparison implementation
// Uses unoptimized SIMD version
tricorn_V3(complexA,width,height,resolution,n,colored,tempArr);
break;
case 4:
// Uses SIMD threaded version
// Uses highly optimized SIMD
tricorn_V4(complexA,width,height,resolution,n,colored,tempArr);
break;
default:
+202
View File
@@ -0,0 +1,202 @@
#include "tricorn.h"
static void* compute(void * args);
static unsigned converge(float c_real, float c_imag, unsigned n);
static size_t get_thread_count();
typedef struct {
unsigned n;
_Bool color;
_Bool flipped;
float start_r;
float start_i;
float res;
size_t start_row;
size_t stop_row;
size_t width;
unsigned char* img;
} compute_args;
/**
* @brief Generates a Tricorn fractal image using an optimized SISD implementation.
* @param start The starting complex coordinate (bottom-left or top-left depending on height sign).
* @param width The width of the target image in pixels.
* @param height The height of the target image. A negative height reverses the Y-axis calculation.
* @param res The step size (resolution) per pixel in the complex plane.
* @param n The maximum number of mathematical iterations per pixel.
* @param color If true, generates a 24-bit RGB image. If false, writes a grayscale image.
* @param img Pointer to the pre-allocated memory buffer where the bitmap data will be stored.
*/
void tricorn(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img){
const ssize_t min_rows_per_thread = 10;
const _Bool flipped = height < 0;
const size_t abs_height = flipped ? -height : height;
// Calculating the threadcount so that it is limited to a reasonable distribution where thread starting doesnt outway calculation cost
size_t thread_count = get_thread_count();
if(height >= min_rows_per_thread)
while((height/thread_count) < min_rows_per_thread) thread_count--;
else thread_count = 1;
const size_t y_step = (height+thread_count-1)/thread_count;
const float start_r = crealf(start);
const float start_i = cimagf(start);
const float directional_res = (flipped ? -1.0f: 1.0f)*res;
pthread_t threads[thread_count];
for(size_t i = 0; i < thread_count; i++){
compute_args* args = malloc(sizeof(*args));
args->start_r =start_r;
args->start_i =start_i-i*directional_res*y_step;
args->start_row = i*y_step;
args->stop_row = abs_height < (args->start_row+y_step) ? abs_height: (args->start_row+y_step);
args->width = width;
args-> flipped = flipped;
args->color = color;
args->res = res;
args->img = img;
args->n = n;
if(pthread_create(&threads[i],NULL,compute,args)){
compute(args);
free(args);
args = NULL;
}
}
for(size_t i = 0; i < thread_count; i++){
pthread_join(threads[i],NULL);
}
return;
}
static void* compute(void * vargs){
compute_args *args = vargs;
vargs = NULL;
const int bytes_per_color = args->color ? 3 : 1;
// Absolute height for iteration purposes
const size_t simd_width = (args->width / 4) * 4;
float c_start_real = args->start_r;
__m128 four = _mm_set1_ps(4.0f);
__m128 minus_two = _mm_set1_ps(-2.0f);
float c_imag = args->start_i;
__m128 c_imag_vector = _mm_set1_ps(c_imag);
float four_res = 4 * args->res;
float vertical_step = (args->flipped? -1.0f: 1.0f) * args->res;
__m128 vertical_step_vector = _mm_set1_ps(vertical_step);
for(size_t y = args->start_row; y < args->stop_row;y++) {
size_t x = 0;
__m128 c_real_vector = _mm_setr_ps(c_start_real,c_start_real + 1.0f*args->res,c_start_real + 2.0f*args->res, c_start_real + 3.0f*args->res);
for(; x < simd_width;x += 4) { //SIMD Loop
__m128 real = _mm_setzero_ps();
__m128 imag = _mm_setzero_ps();
__m128i iterations = _mm_setzero_si128();
__m128i active_mask = _mm_set1_epi32(-1);
for(unsigned i = 0; i < args->n;i++) {
__m128 real2 = _mm_mul_ps(real,real);
__m128 imag2 = _mm_mul_ps(imag,imag);
__m128 unbounded = _mm_cmple_ps(_mm_add_ps(real2,imag2),four);
__m128i unbounded_i = _mm_castps_si128(unbounded);
active_mask = _mm_and_si128(active_mask,unbounded_i);
if(_mm_test_all_zeros(active_mask,active_mask)) {
break;
}
iterations = _mm_sub_epi32(iterations,active_mask);
imag = _mm_add_ps(c_imag_vector,_mm_mul_ps(minus_two, _mm_mul_ps(real,imag)));
real = _mm_add_ps(_mm_sub_ps(real2,imag2),c_real_vector);
}
__m128i mod64 = _mm_and_si128(iterations, _mm_set1_epi32(63));
__m128i times4 = _mm_slli_epi32(mod64, 2);
__m128i values = _mm_sub_epi32(_mm_set1_epi32(255), times4);
__m128i max_iter = _mm_set1_epi32(args->n);
__m128i is_max = _mm_cmpeq_epi32(iterations, max_iter);
values = _mm_andnot_si128(is_max, values);
unsigned int value_array[4];
_mm_storeu_si128((__m128i*)value_array,values);
for(unsigned i = 0 ; i < 4;i++) {
unsigned int value = value_array[i];
size_t position = (y*args->width + (x+i)) * bytes_per_color;
if(!args->color) {
args->img[position] = value;
}
else {
args->img[position] = ((value*3) % 255);
args->img[position + 1] = ((value*2) % 255);
args->img[position + 2] = ((value*5) % 255);
}
}
c_real_vector = _mm_add_ps(c_real_vector,_mm_set1_ps(four_res));
}
for(;x < args->width; x++) {
float c_real = c_start_real + (x*args->res);
unsigned int value = converge(c_real,c_imag,args->n);
size_t position = (y*args->width + x) * bytes_per_color;
if(!args->color) {
args->img[position] = value;
}
else {
args->img[position] = ((value*3) % 255);
args->img[position + 1] = ((value*2) % 255);
args->img[position + 2] = ((value*5) % 255);
}
}
c_imag -= vertical_step;
c_imag_vector = _mm_sub_ps(c_imag_vector,vertical_step_vector);
}
free(args);
args = NULL;
return NULL;
}
// Iterative application of the equation given by the assignment
static unsigned converge(float c_real, float c_imag, unsigned n) {
float real = 0.0f;
float imag = 0.0f;
for (unsigned i = 0; i < n; i++) {
float real2 = real * real;
float imag2 = imag * imag;
if (real2 + imag2 > 4.0f) {
return 255 - 4 * (i % 64);
}
imag = -2.0f * real * imag + c_imag;
real = real2 - imag2 + c_real;
}
return 0;
}
static inline size_t get_thread_count(){
#ifdef WIN32
SYSTEM_INFO sysinfo;
GetSystemInfo(&sysinfo);
return sysinfo.dwNumberOfProcessors;
#else
return sysconf(_SC_NPROCESSORS_ONLN);
#endif
}
+13
View File
@@ -0,0 +1,13 @@
#include <complex.h>
#include <stddef.h>
#include <math.h>
#include <sys/types.h>
#include <immintrin.h>
#include <pthread.h>
#ifdef WIN32
#include <sysinfoapi.h>
#else
#include <unistd.h>
#endif
void tricorn(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img);
-130
View File
@@ -1,130 +0,0 @@
#include "tricorn_V0.h"
#include <immintrin.h>
static unsigned converge(float c_real, float c_imag, unsigned n);
/**
* @brief Generates a Tricorn fractal image using highly optimized SSE/SIMD instructions and optimized code.
* @param start The starting complex coordinate (bottom-left or top-left depending on height sign).
* @param width The width of the target image in pixels.
* @param height The height of the target image. A negative height reverses the Y-axis calculation.
* @param res The step size (resolution) per pixel in the complex plane.
* @param n The maximum number of mathematical iterations per pixel.
* @param color If true, generates a 24-bit RGB image. If false, writes a grayscale image.
* @param img Pointer to the pre-allocated memory buffer where the bitmap data will be stored.
*/
void tricorn_V0(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img){
int bytesPerColor = color ? 3 : 1;
// Absolute height for iteration purposes
ssize_t actual_height = height < 0 ? (ssize_t)(-height) : height;
ssize_t imag_step = height < 0 ? -1:1;
size_t simd_width = (width / 4) * 4;
float c_start_real = creal(start);
float c_start_imag = cimag(start);
__m128 four = _mm_set1_ps(4.0f);
__m128 minus_two = _mm_set1_ps(-2.0f);
float c_imag = c_start_imag;
__m128 c_imag_vector = _mm_set1_ps(c_imag);
float four_res = 4 * res;
float vertical_step = imag_step * res;
__m128 vertical_step_vector = _mm_set1_ps(vertical_step);
for(ssize_t y = 0; y < actual_height;y++) {
size_t x = 0;
__m128 c_real_vector = _mm_setr_ps(c_start_real,c_start_real + 1.0f*res,c_start_real + 2.0f*res, c_start_real + 3.0f*res);
for(; x < simd_width;x += 4) { //SIMD Loop
__m128 real = _mm_setzero_ps();
__m128 imag = _mm_setzero_ps();
__m128i iterations = _mm_setzero_si128();
__m128i active_mask = _mm_set1_epi32(-1);
for(unsigned i = 0; i < n;i++) {
__m128 real2 = _mm_mul_ps(real,real);
__m128 imag2 = _mm_mul_ps(imag,imag);
__m128 unbounded = _mm_cmple_ps(_mm_add_ps(real2,imag2),four);
__m128i unbounded_i = _mm_castps_si128(unbounded);
active_mask = _mm_and_si128(active_mask,unbounded_i);
if(_mm_test_all_zeros(active_mask,active_mask)) {
break;
}
iterations = _mm_sub_epi32(iterations,active_mask);
imag = _mm_add_ps(c_imag_vector,_mm_mul_ps(minus_two, _mm_mul_ps(real,imag)));
real = _mm_add_ps(_mm_sub_ps(real2,imag2),c_real_vector);
}
__m128i mod64 = _mm_and_si128(iterations, _mm_set1_epi32(63));
__m128i times4 = _mm_slli_epi32(mod64, 2);
__m128i values = _mm_sub_epi32(_mm_set1_epi32(255), times4);
__m128i max_iter = _mm_set1_epi32(n);
__m128i is_max = _mm_cmpeq_epi32(iterations, max_iter);
values = _mm_andnot_si128(is_max, values);
unsigned int value_array[4];
_mm_storeu_si128((__m128i*)value_array,values);
for(unsigned i = 0 ; i < 4;i++) {
unsigned int value = value_array[i];
size_t position = (y*width + (x+i)) * bytesPerColor;
if(!color) {
img[position] = value;
}
else {
img[position] = ((value*3) % 255);
img[position + 1] = ((value*2) % 255);
img[position + 2] = ((value*5) % 255);
}
}
c_real_vector = _mm_add_ps(c_real_vector,_mm_set1_ps(four_res));
}
for(;x < width; x++) {
float c_real = c_start_real + (x*res);
unsigned int value = converge(c_real,c_imag,n);
size_t position = (y*width + x) * bytesPerColor;
if(!color) {
img[position] = value;
}
else {
img[position] = ((value*3) % 255);
img[position + 1] = ((value*2) % 255);
img[position + 2] = ((value*5) % 255);
}
}
c_imag -= vertical_step;
c_imag_vector = _mm_sub_ps(c_imag_vector,vertical_step_vector);
}
}
static unsigned converge(float c_real, float c_imag, unsigned n) {
float real = 0.0f;
float imag = 0.0f;
for (unsigned i = 0; i < n; i++) {
float real2 = real * real;
float imag2 = imag * imag;
if (real2 + imag2 > 4.0f) {
return 255 - 4 * (i % 64);
}
imag = -2.0f * real * imag + c_imag;
real = real2 - imag2 + c_real;
}
return 0;
}
-6
View File
@@ -1,6 +0,0 @@
#include <complex.h>
#include <stddef.h>
#include <math.h>
#include <sys/types.h>
#include<stdbool.h>
void tricorn_V0(float complex start, size_t width, ssize_t height, float res, unsigned n, bool color, unsigned char* img);
+25 -90
View File
@@ -1,11 +1,10 @@
#include "tricorn_V1.h"
#include <immintrin.h>
static unsigned converge(float c_rea, float c_imag, unsigned n);
static const float convergenceLimit = 2.0f; // Mandelbrot set limit
static unsigned converge(float complex position,unsigned n);
/**
* @brief Generates a Tricorn fractal image using highly optimized SSE/SIMD instructions, however the algoriithm is not fully optimized.
* @brief Generates a Tricorn fractal image using the standard, unoptimized scalar implementation for benchmarking comparison (Version 2).
* @param start The starting complex coordinate (bottom-left or top-left depending on height sign).
* @param width The width of the target image in pixels.
* @param height The height of the target image. A negative height reverses the Y-axis calculation.
@@ -15,112 +14,48 @@ static unsigned converge(float c_rea, float c_imag, unsigned n);
* @param img Pointer to the pre-allocated memory buffer where the bitmap data will be stored.
*/
void tricorn_V1(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img){
int bytesPerColor = color ? 3 : 1;
// Absolute height for iteration purposes
ssize_t actual_height = height < 0 ? (ssize_t)(-height) : height;
ssize_t imag_step = height < 0 ? -1:1;
size_t simd_width = (width / 4) * 4;
float imag_direction = height < 0 ? -1.0f : 1.0f;
unsigned bytesPerColor = color ? 3 : 1;
height = height < 0? -height : height;
float c_start_real = creal(start);
float c_start_imag = cimag(start);
// Simple iteration over every pixel in the picture
for(size_t i = 0; i < width; i++){
for(ssize_t j = 0; j < height; j++){
__m128 four = _mm_set1_ps(4.0f);
__m128 minus_two = _mm_set1_ps(-2.0f);
unsigned value = converge(start + (i*res - I* imag_direction*j*res),n);
for(ssize_t y = 0; y < actual_height;y++) {
float c_imag = c_start_imag - (y*imag_step) * res;
__m128 c_imag_vector = _mm_set1_ps(c_imag);
size_t x = 0;
for(; x < simd_width;x += 4) { //SIMD Loop
float c_real = c_start_real + (x*res);
__m128 c_real_vector = _mm_setr_ps(c_real,c_real+1.0f*res,c_real+2.0f*res,c_real+3.0f*res);
__m128 real = _mm_setzero_ps();
__m128 imag = _mm_setzero_ps();
__m128i iterations = _mm_setzero_si128();
__m128i active_mask = _mm_set1_epi32(-1);
for(unsigned i = 0; i < n;i++) {
__m128 real2 = _mm_mul_ps(real,real);
__m128 imag2 = _mm_mul_ps(imag,imag);
__m128 unbounded = _mm_cmple_ps(_mm_add_ps(real2,imag2),four);
__m128i unbounded_i = _mm_castps_si128(unbounded);
active_mask = _mm_and_si128(active_mask,unbounded_i);
if(_mm_test_all_zeros(active_mask,active_mask)) {
break;
}
iterations = _mm_sub_epi32(iterations,active_mask);
imag = _mm_add_ps(c_imag_vector,_mm_mul_ps(minus_two, _mm_mul_ps(real,imag)));
real = _mm_add_ps(_mm_sub_ps(real2,imag2),c_real_vector);
}
// Pixel value position in the output array
size_t position = (j*width +i) * bytesPerColor;
__m128i mod64 = _mm_and_si128(iterations, _mm_set1_epi32(63));
__m128i times4 = _mm_slli_epi32(mod64, 2);
__m128i values = _mm_sub_epi32(_mm_set1_epi32(255), times4);
__m128i max_iter = _mm_set1_epi32(n);
__m128i is_max = _mm_cmpeq_epi32(iterations, max_iter);
values = _mm_andnot_si128(is_max, values);
unsigned int value_array[4];
_mm_storeu_si128((__m128i*)value_array,values);
for(unsigned i = 0 ; i < 4;i++) {
unsigned int value = value_array[i];
size_t position = (y*width + (x+i)) * bytesPerColor;
if(!color) {
if(!color){ // Grayscale colors as specified by the assignment r = g = b
img[position] = value;
}
else {
else{ // Temporary colors with random multipliers
img[position] = ((value*3) % 255);
img[position + 1] = ((value*2) % 255) ;
img[position + 2] = ((value*5) % 255);
}
}
}
for(;x < width; x++) {
float c_real = c_start_real + (x*res);
unsigned int value = converge(c_real,c_imag,n);
size_t position = (y*width + x) * bytesPerColor;
if(!color) {
img[position] = value;
}
else {
img[position] = ((value*3) % 255);
img[position + 1] = ((value*2) % 255);
img[position + 2] = ((value*5) % 255);
}
}
}
}
}
static unsigned converge(float c_real, float c_imag, unsigned n) {
float real = 0.0f;
float imag = 0.0f;
// Iterative application of the equation given by the assignment
static unsigned converge(float complex position,unsigned n){
complex float result = 0;
for(unsigned i = 0; i < n; i++){
float real2 = real * real;
float imag2 = imag * imag;
if (real2 + imag2 > 4.0f) {
// Detection of unbounded behaviour
// If the magnitude of the result vector is greater than the convergence limit, we know that the function will grow unbounded
// Current convergence limit taken from mandelbrot set
if(sqrtf(powf (crealf(result),2)+powf(cimagf(result),2)) > convergenceLimit)
return 255-4*(i % 64);
}
imag = -2.0f * real * imag + c_imag;
real = real2 - imag2 + c_real;
// Application of the given complex equation with m = 2
result = cpowf(conjf(result),2) + position;
}
return 0;
}
+1 -2
View File
@@ -2,5 +2,4 @@
#include <stddef.h>
#include <math.h>
#include <sys/types.h>
#include<stdbool.h>
void tricorn_V1(float complex start, size_t width, ssize_t height, float res, unsigned n, bool color, unsigned char* img);
void tricorn_V1(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img);
+91 -26
View File
@@ -1,10 +1,11 @@
#include "tricorn_V3.h"
#include <immintrin.h>
static unsigned converge(float c_rea, float c_imag, unsigned n);
static const float convergenceLimit = 2.0f; // Mandelbrot set limit
static unsigned converge(float complex position,unsigned n);
/**
* @brief Generates a Tricorn fractal image using the standard, unoptimized scalar implementation for benchmarking comparison (Version 2).
* @brief Generates a Tricorn fractal image using highly optimized SSE/SIMD instructions, however the algoriithm is not fully optimized.
* @param start The starting complex coordinate (bottom-left or top-left depending on height sign).
* @param width The width of the target image in pixels.
* @param height The height of the target image. A negative height reverses the Y-axis calculation.
@@ -14,48 +15,112 @@ static unsigned converge(float complex position,unsigned n);
* @param img Pointer to the pre-allocated memory buffer where the bitmap data will be stored.
*/
void tricorn_V3(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img){
int bytesPerColor = color ? 3 : 1;
float imag_direction = height < 0 ? -1.0f : 1.0f;
unsigned bytesPerColor = color ? 3 : 1;
height = height < 0? -height : height;
// Absolute height for iteration purposes
ssize_t actual_height = height < 0 ? (ssize_t)(-height) : height;
ssize_t imag_step = height < 0 ? -1:1;
size_t simd_width = (width / 4) * 4;
// Simple iteration over every pixel in the picture
for(size_t i = 0; i < width; i++){
for(ssize_t j = 0; j < height; j++){
float c_start_real = creal(start);
float c_start_imag = cimag(start);
unsigned value = converge(start + (i*res - I* imag_direction*j*res),n);
__m128 four = _mm_set1_ps(4.0f);
__m128 minus_two = _mm_set1_ps(-2.0f);
// Pixel value position in the output array
size_t position = (j*width +i) * bytesPerColor;
for(ssize_t y = 0; y < actual_height;y++) {
float c_imag = c_start_imag - (y*imag_step) * res;
__m128 c_imag_vector = _mm_set1_ps(c_imag);
size_t x = 0;
for(; x < simd_width;x += 4) { //SIMD Loop
float c_real = c_start_real + (x*res);
__m128 c_real_vector = _mm_setr_ps(c_real,c_real+1.0f*res,c_real+2.0f*res,c_real+3.0f*res);
__m128 real = _mm_setzero_ps();
__m128 imag = _mm_setzero_ps();
__m128i iterations = _mm_setzero_si128();
__m128i active_mask = _mm_set1_epi32(-1);
for(unsigned i = 0; i < n;i++) {
__m128 real2 = _mm_mul_ps(real,real);
__m128 imag2 = _mm_mul_ps(imag,imag);
__m128 unbounded = _mm_cmple_ps(_mm_add_ps(real2,imag2),four);
__m128i unbounded_i = _mm_castps_si128(unbounded);
active_mask = _mm_and_si128(active_mask,unbounded_i);
if(_mm_test_all_zeros(active_mask,active_mask)) {
break;
}
iterations = _mm_sub_epi32(iterations,active_mask);
imag = _mm_add_ps(c_imag_vector,_mm_mul_ps(minus_two, _mm_mul_ps(real,imag)));
real = _mm_add_ps(_mm_sub_ps(real2,imag2),c_real_vector);
}
if(!color){ // Grayscale colors as specified by the assignment r = g = b
__m128i mod64 = _mm_and_si128(iterations, _mm_set1_epi32(63));
__m128i times4 = _mm_slli_epi32(mod64, 2);
__m128i values = _mm_sub_epi32(_mm_set1_epi32(255), times4);
__m128i max_iter = _mm_set1_epi32(n);
__m128i is_max = _mm_cmpeq_epi32(iterations, max_iter);
values = _mm_andnot_si128(is_max, values);
unsigned int value_array[4];
_mm_storeu_si128((__m128i*)value_array,values);
for(unsigned i = 0 ; i < 4;i++) {
unsigned int value = value_array[i];
size_t position = (y*width + (x+i)) * bytesPerColor;
if(!color) {
img[position] = value;
}
else{ // Temporary colors with random multipliers
else {
img[position] = ((value*3) % 255);
img[position + 1] = ((value*2) % 255);
img[position + 2] = ((value*5) % 255);
}
}
}
}
// Iterative application of the equation given by the assignment
static unsigned converge(float complex position,unsigned n){
for(;x < width; x++) {
float c_real = c_start_real + (x*res);
unsigned int value = converge(c_real,c_imag,n);
size_t position = (y*width + x) * bytesPerColor;
if(!color) {
img[position] = value;
}
else {
img[position] = ((value*3) % 255);
img[position + 1] = ((value*2) % 255);
img[position + 2] = ((value*5) % 255);
}
}
}
}
static unsigned converge(float c_real, float c_imag, unsigned n) {
float real = 0.0f;
float imag = 0.0f;
complex float result = 0;
for (unsigned i = 0; i < n; i++) {
// Detection of unbounded behaviour
// If the magnitude of the result vector is greater than the convergence limit, we know that the function will grow unbounded
// Current convergence limit taken from mandelbrot set
if(sqrtf(powf (crealf(result),2)+powf(cimagf(result),2)) > convergenceLimit)
float real2 = real * real;
float imag2 = imag * imag;
if (real2 + imag2 > 4.0f) {
return 255 - 4 * (i % 64);
}
// Application of the given complex equation with m = 2
result = cpowf(conjf(result),2) + position;
imag = -2.0f * real * imag + c_imag;
real = real2 - imag2 + c_real;
}
return 0;
}
+2 -1
View File
@@ -2,4 +2,5 @@
#include <stddef.h>
#include <math.h>
#include <sys/types.h>
void tricorn_V3(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img);
#include<stdbool.h>
void tricorn_V3(float complex start, size_t width, ssize_t height, float res, unsigned n, bool color, unsigned char* img);
+29 -101
View File
@@ -1,22 +1,8 @@
#include "tricorn_V4.h"
static void* compute(void * args);
static unsigned converge(float c_real, float c_imag, unsigned n);
static size_t get_thread_count();
typedef struct {
unsigned n;
_Bool color;
_Bool flipped;
float start_r;
float start_i;
float res;
size_t start_row;
size_t stop_row;
size_t width;
unsigned char* img;
} compute_args;
/**
* @brief Generates a Tricorn fractal image using an optimized SIMD + Threading implementation.
* @param start The starting complex coordinate (bottom-left or top-left depending on height sign).
@@ -28,81 +14,36 @@ typedef struct {
* @param img Pointer to the pre-allocated memory buffer where the bitmap data will be stored.
*/
void tricorn_V4(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img){
int bytesPerColor = color ? 3 : 1;
const ssize_t min_rows_per_thread = 10;
_Bool flipped = height < 0;
size_t abs_height = flipped ? -height : height;
size_t thread_count = get_thread_count();
if(height >= min_rows_per_thread)
while((height/thread_count) < min_rows_per_thread) thread_count--;
else thread_count = 1;
const size_t y_step = (height+thread_count-1)/thread_count;
const float start_r = crealf(start);
const float start_i = cimagf(start);
float directional_res = (flipped ? -1.0f: 1.0f)*res;
pthread_t threads[thread_count];
for(size_t i = 0; i < thread_count; i++){
compute_args* args = malloc(sizeof(*args));
args->start_r =start_r;
args->start_i =start_i-i*directional_res*y_step;
args->start_row = i*y_step;
args->stop_row = abs_height < (args->start_row+y_step) ? abs_height: (args->start_row+y_step);
args->width = width;
args-> flipped = flipped;
args->color = color;
args->res = res;
args->img = img;
args->n = n;
if(pthread_create(&threads[i],NULL,compute,args)){
compute(args);
free(args);
args = NULL;
}
}
for(size_t i = 0; i < thread_count; i++){
pthread_join(threads[i],NULL);
}
return;
}
static void* compute(void * vargs){
compute_args *args = vargs;
vargs = NULL;
int bytes_per_color = args->color ? 3 : 1;
// Absolute height for iteration purposes
size_t simd_width = (args->width / 4) * 4;
ssize_t actual_height = height < 0 ? (ssize_t)(-height) : height;
ssize_t imag_step = height < 0 ? -1:1;
size_t simd_width = (width / 4) * 4;
float c_start_real = creal(start);
float c_start_imag = cimag(start);
float c_start_real = args->start_r;
__m128 four = _mm_set1_ps(4.0f);
__m128 minus_two = _mm_set1_ps(-2.0f);
float c_imag = args->start_i;
float c_imag = c_start_imag;
__m128 c_imag_vector = _mm_set1_ps(c_imag);
float four_res = 4 * args->res;
float vertical_step = (args->flipped? -1.0f: 1.0f) * args->res;
float four_res = 4 * res;
float vertical_step = imag_step * res;
__m128 vertical_step_vector = _mm_set1_ps(vertical_step);
for(size_t y = args->start_row; y < args->stop_row;y++) {
for(ssize_t y = 0; y < actual_height;y++) {
size_t x = 0;
__m128 c_real_vector = _mm_setr_ps(c_start_real,c_start_real + 1.0f*args->res,c_start_real + 2.0f*args->res, c_start_real + 3.0f*args->res);
__m128 c_real_vector = _mm_setr_ps(c_start_real,c_start_real + 1.0f*res,c_start_real + 2.0f*res, c_start_real + 3.0f*res);
for(; x < simd_width;x += 4) { //SIMD Loop
__m128 real = _mm_setzero_ps();
__m128 imag = _mm_setzero_ps();
__m128i iterations = _mm_setzero_si128();
__m128i active_mask = _mm_set1_epi32(-1);
for(unsigned i = 0; i < args->n;i++) {
for(unsigned i = 0; i < n;i++) {
__m128 real2 = _mm_mul_ps(real,real);
__m128 imag2 = _mm_mul_ps(imag,imag);
@@ -126,7 +67,7 @@ static void* compute(void * vargs){
__m128i times4 = _mm_slli_epi32(mod64, 2);
__m128i values = _mm_sub_epi32(_mm_set1_epi32(255), times4);
__m128i max_iter = _mm_set1_epi32(args->n);
__m128i max_iter = _mm_set1_epi32(n);
__m128i is_max = _mm_cmpeq_epi32(iterations, max_iter);
values = _mm_andnot_si128(is_max, values);
@@ -135,44 +76,40 @@ static void* compute(void * vargs){
for(unsigned i = 0 ; i < 4;i++) {
unsigned int value = value_array[i];
size_t position = (y*args->width + (x+i)) * bytes_per_color;
size_t position = (y*width + (x+i)) * bytesPerColor;
if(!args->color) {
args->img[position] = value;
if(!color) {
img[position] = value;
}
else {
args->img[position] = ((value*3) % 255);
args->img[position + 1] = ((value*2) % 255);
args->img[position + 2] = ((value*5) % 255);
img[position] = ((value*3) % 255);
img[position + 1] = ((value*2) % 255);
img[position + 2] = ((value*5) % 255);
}
}
c_real_vector = _mm_add_ps(c_real_vector,_mm_set1_ps(four_res));
}
for(;x < args->width; x++) {
float c_real = c_start_real + (x*args->res);
unsigned int value = converge(c_real,c_imag,args->n);
for(;x < width; x++) {
float c_real = c_start_real + (x*res);
unsigned int value = converge(c_real,c_imag,n);
size_t position = (y*args->width + x) * bytes_per_color;
size_t position = (y*width + x) * bytesPerColor;
if(!args->color) {
args->img[position] = value;
if(!color) {
img[position] = value;
}
else {
args->img[position] = ((value*3) % 255);
args->img[position + 1] = ((value*2) % 255);
args->img[position + 2] = ((value*5) % 255);
img[position] = ((value*3) % 255);
img[position + 1] = ((value*2) % 255);
img[position + 2] = ((value*5) % 255);
}
}
c_imag -= vertical_step;
c_imag_vector = _mm_sub_ps(c_imag_vector,vertical_step_vector);
}
free(args);
args = NULL;
return NULL;
}
// Iterative application of the equation given by the assignment
static unsigned converge(float c_real, float c_imag, unsigned n) {
float real = 0.0f;
float imag = 0.0f;
@@ -190,12 +127,3 @@ static unsigned converge(float c_real, float c_imag, unsigned n) {
}
return 0;
}
static inline size_t get_thread_count(){
#ifdef WIN32
SYSTEM_INFO sysinfo;
GetSystemInfo(&sysinfo);
return sysinfo.dwNumberOfProcessors;
#else
return sysconf(_SC_NPROCESSORS_ONLN);
#endif
}
+2 -8
View File
@@ -2,12 +2,6 @@
#include <stddef.h>
#include <math.h>
#include <sys/types.h>
#include<stdbool.h>
#include <immintrin.h>
#include <pthread.h>
#ifdef WIN32
#include <sysinfoapi.h>
#else
#include <unistd.h>
#endif
void tricorn_V4(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img);
void tricorn_V4(float complex start, size_t width, ssize_t height, float res, unsigned n, bool color, unsigned char* img);