Added: version 5 for testing, example test of inlining conversion function included but slower

Implementierung/main.c

@@ -12,6 +12,7 @@
 #include "tricorn_V2.h"
 #include "tricorn_V3.h"
 #include "tricorn_V4.h"
+#include "tricorn_V5.h"
 
 
 void save_bmp(const char *filename, unsigned char *img, size_t width, ssize_t height,_Bool color);
@@ -229,6 +230,10 @@ int main(int argc, char** argv)
             // Uses highly optimized SIMD
             tricorn_V4(complexA,width,height,resolution,n,colored,tempArr);
             break;
+        case 5:
+            // Experimental cleanups and inlines / for testing
+            tricorn_V5(complexA,width,height,resolution,n,colored,tempArr);
+            break;
         default:
             break;
     }

Implementierung/tricorn_V5.c

@@ -0,0 +1,216 @@
+#include "tricorn_V5.h"
+// Argument structs passed to the threads for the caluclations
+typedef struct{
+    unsigned n;
+    _Bool color;
+    _Bool flipped;
+    float res;
+    size_t width;
+    unsigned char* img;
+
+} const_args;
+typedef struct {
+    float start_r;
+    float start_i;
+    size_t start_row;
+    size_t stop_row;
+    const_args* consts;
+} compute_args;
+
+static  void* compute(void * args);
+static unsigned converge(float c_real, float c_imag, unsigned n);
+static inline void simd_converge(unsigned int (*value_array)[] ,__m128 c_real_vector,__m128 c_imag_vector,__m128 four,__m128 minus_two,const_args* c);
+static size_t get_thread_count();
+
+/**
+ * @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_V5(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;
+
+    const_args* consts = malloc(sizeof(*consts));
+    consts->width = width;
+    consts->flipped = flipped;
+    consts->color = color;
+    consts->res = res;
+    consts->img = img;
+    consts->n = n;
+
+    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->consts = consts;
+        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);
+    }
+
+    free(consts);
+    consts = NULL;
+    return;
+}
+
+static void* compute(void * vargs){
+
+    compute_args *args = vargs;
+    const_args *c = args->consts;
+    vargs = NULL;
+
+
+    const int bytes_per_color = c->color ? 3 : 1;
+    // Absolute height for iteration purposes
+    const size_t simd_width = (c->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 * c->res;
+    float vertical_step = (c->flipped? -1.0f: 1.0f) * c->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*c->res,c_start_real + 2.0f*c->res, c_start_real + 3.0f*c->res);
+        for(; x < simd_width;x += 4) { //SIMD Loop
+
+            unsigned int value_array[4];
+            simd_converge(&value_array,c_real_vector,c_imag_vector,four,minus_two,c);
+
+            for(unsigned i = 0 ; i < 4;i++) {
+                unsigned int value = value_array[i];
+                size_t position = (y*c->width + (x+i)) * bytes_per_color;
+
+                if(!c->color) {
+                    c->img[position] = value;
+                }
+                else {
+                    c->img[position] = ((value*3) % 255);
+                    c->img[position + 1] = ((value*2) % 255);
+                    c->img[position + 2] = ((value*5) % 255);
+                }
+            }
+            c_real_vector = _mm_add_ps(c_real_vector,_mm_set1_ps(four_res));
+        }
+
+        for(;x < c->width; x++) {
+            float c_real = c_start_real + (x*c->res);
+            unsigned int value = converge(c_real,c_imag,c->n);
+
+            size_t position = (y*c->width + x) * bytes_per_color;
+
+            if(!c->color) {
+                c->img[position] = value;
+            }
+            else {
+                c->img[position] = ((value*3) % 255);
+                c->img[position + 1] = ((value*2) % 255);
+                c->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;
+}
+static inline void simd_converge(unsigned int (*value_array)[] ,__m128 c_real_vector,__m128 c_imag_vector,__m128 four,__m128 minus_two,const_args* c){
+    __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 < c->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(c->n);
+    __m128i is_max = _mm_cmpeq_epi32(iterations, max_iter);
+    values = _mm_andnot_si128(is_max, values);
+    _mm_storeu_si128((__m128i*)value_array,values);
+
+}
+// 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; 
+}
+// Platform independent function returning the amount of threads available on the systems cpu
+static inline size_t get_thread_count(){
+    #ifdef WIN32
+    SYSTEM_INFO sysinfo;
+    GetSystemInfo(&sysinfo);
+    return sysinfo.dwNumberOfProcessors;
+    #else
+    return sysconf(_SC_NPROCESSORS_ONLN);
+    #endif
+}

Implementierung/tricorn_V5.h

@@ -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_V5(float complex  start, size_t width, ssize_t  height, float res, unsigned n, _Bool color, unsigned char* img);