Merge branch 'tricorn_simd'

2026-06-28 21:02:13 +02:00
parent dd32e43309 a3d829533c
commit fb11fa9302
5 changed files with 144 additions and 10 deletions
@@ -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

@@ -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"
@@ -198,13 +199,24 @@ int main(int argc, char** argv)
        bytesPerColour = 1;
    }
    
-    unsigned char tempArr[width*height*bytesPerColour]; // if not colored, rem
+    ssize_t absolute_height = height >= 0 ? height : height * - 1;
+    unsigned char *tempArr = malloc(width*absolute_height*bytesPerColour);// allocate image array on the heap to support large images as the stack is usually only 8 MB big
+    if(tempArr == NULL) {
+        exit(EXIT_FAILURE);
+    }
+    
    printf("yashar\n");
-    tricorn(complexA,width,height,resolution,n,colored,tempArr);
+    if(version_to_use == 0) {
+        tricorn_simd(complexA,width,height,resolution,n,colored,tempArr);
+    } 
+    else {
+        tricorn(complexA,width,height,resolution,n,colored,tempArr);
+    }
    printf("yashar\n");
    save_bmp(output_data,tempArr,width,height,colored);
    
-    free(output_data); // 
+    free(output_data);
+    free(tempArr);
    return 0;
 }

@@ -218,14 +230,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] = { 
@@ -1,7 +1,7 @@
 #include "tricorn.h"

 static const float convergenceLimit = 2.0f; // Mandelbrot set limit
-unsigned converge(float complex position,unsigned n);
+static unsigned converge(float complex position,unsigned n);

 void tricorn(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img){
    
@@ -34,7 +34,7 @@ void tricorn(float complex start, size_t width, ssize_t height, float res, unsig
 }

 // Iterative application of the equation given by the assignment
-unsigned converge(float complex position,unsigned n){
+static unsigned converge(float complex position,unsigned n){

    complex float result  = 0;
    for(unsigned i = 0; i < n; i++){
@@ -0,0 +1,115 @@
+#include "tricorn_simd.h"
+#include <immintrin.h>
+
+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
+    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);
+
+    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 = 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_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);
+                }
+            }
+        }
+
+        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;
+
+    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; 
+}
@@ -0,0 +1,5 @@
+#include <complex.h>
+#include <stddef.h>
+#include <math.h>
+#include <sys/types.h>
+void tricorn_simd(float complex  start, size_t width, ssize_t  height, float res, unsigned n, _Bool color, unsigned char* img);