Implemented SIMD optimizations in tricorn_simd including negative height handling which needs to be fixed in normal tricorn.c

This commit is contained in:
go38coz
2026-06-26 22:38:55 +02:00
parent ee2b6b9506
commit 791231a4f5
3 changed files with 102 additions and 35 deletions
+1 -1
View File
@@ -1,5 +1,5 @@
CC = gcc
CFLAGS = -Wall -Wextra -O3 -std=c17
CFLAGS = -Wall -Wextra -O3 -std=c17 -msse4.1
DEPFLAGS = -MMD -MP
LDFLAGS = -lm
+8 -5
View File
@@ -6,6 +6,7 @@
#include <getopt.h> // for getopt_long
#include <string.h> // for strcpy
#include "tricorn.h"
#include "tricorn_simd.h"
void save_bmp(const char *filename, unsigned char *img, size_t width, ssize_t height,_Bool color);
@@ -17,7 +18,7 @@ unsigned tests(void) {
void benchmarking(void) {
// do the benchmarking
}
// REMOVE LATER: interesting test command ./tricorn -s -2.0,0.6 -r 0.0010 -n 200 -d 1920 1080 -o image.bmp -C
// REMOVE LATER: interesting test command ./tricorn -s -2.0,0.6 -r 0.0010 -n 200 -d 1920,-1080 -o image.bmp -C
const char* help_messsage = "Help Message\n\n"
"Usage:\n"
" ./tricorn [OPTIONS]\n\n"
@@ -200,7 +201,7 @@ int main(int argc, char** argv)
unsigned char tempArr[width*height*bytesPerColour]; // if not colored, rem
printf("yashar\n");
tricorn(complexA,width,height,resolution,n,colored,tempArr);
tricorn_simd(complexA,width,height,resolution,n,colored,tempArr);
printf("yashar\n");
save_bmp(output_data,tempArr,width,height,colored);
@@ -218,14 +219,16 @@ void save_bmp(const char *filename, unsigned char *img, size_t width, ssize_t he
exit(EXIT_FAILURE);
}
size_t bytesPerColor = color ? 3 : 1;
int actual_height = (height < 0) ? -height : height; // Calculate absolute height for iteration and calculation purposes
size_t row_padding = (4 - (width) * 3 % 4) % 4; // Calculate padding for each row as .bmp is 4 bytes aligned
size_t row_width = width * 3 + row_padding; // Calculate row width with padding
size_t row_padding = (4 - (width * bytesPerColor) % 4) % 4; // Calculate padding for each row as .bmp is 4 bytes aligned
size_t row_width = width * bytesPerColor + row_padding; // Calculate row width with padding
size_t image_size = actual_height * row_width;
size_t file_size = image_size + 54; // File size = Image Size + Headers
size_t bytesPerColor = color ? 3 : 1;
// File header - contains information about the file
unsigned char file_header[14] = {
+90 -26
View File
@@ -1,51 +1,115 @@
#include "tricorn.h"
#include "tricorn_simd.h"
#include <immintrin.h>
static const float convergenceLimit = 2.0f; // Mandelbrot set limit
static unsigned converge(float complex position,unsigned n);
static unsigned converge(float c_rea, float c_imag, unsigned n);
void tricorn_simd(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
height = height < 0 ? height*-1 : height;
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*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) *3;
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
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);
if(!color){ // Grayscale colors as specified by the assignment r = g = b
img[position] = value;
img[position + 1] = value;
img[position + 2] = value;
__m128 real = c_real_vector;
__m128 imag = c_imag_vector;
__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_cmplt_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;
}
else{ // Temporary colors with random multipliers
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);
}
}
}
}
// 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++) {
// Application of the given complex equation with m = 2
result = cpowf(conjf(result),2) + position;
float real2 = real * real;
float imag2 = imag * imag;
// 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)
if (real2 + imag2 > 4.0f) {
return 255 - 4 * (i % 64);
}
imag = -2.0f * real * imag + c_imag;
real = real2 - imag2 + c_real;
}
return 0;
}