Optimized SIMD implementation - V0

Shifted V1,V2 to V2,V3 and made V1 be the old not fully optimized SIMD implementation

.gitignore

@@ -4,3 +4,4 @@ tricorn
 *.o
 *.d
 *.data
+*.data.old

Implementierung/main.c

@@ -8,6 +8,7 @@
 #include "tricorn_V0.h"
 #include "tricorn_V1.h"
 #include "tricorn_V2.h"
+#include "tricorn_V3.h"
 
 void save_bmp(const char *filename, unsigned char *img, size_t width, ssize_t height,_Bool color);
 
@@ -208,13 +209,17 @@ int main(int argc, char** argv)
             tricorn_V0(complexA,width,height,resolution,n,colored,tempArr);
             break;
         case 1: 
-            //Uses optimized SISD version
+            // Uses unoptimized SIMD version
             tricorn_V1(complexA,width,height,resolution,n,colored,tempArr);
             break;
         case 2:
-            //Uses unoptimised SISD comparison implementation
+            //Uses optimized SISD version
             tricorn_V2(complexA,width,height,resolution,n,colored,tempArr);
             break;
+        case 3:
+            //Uses unoptimised comparison implementation
+            tricorn_V3(complexA,width,height,resolution,n,colored,tempArr);
+            break;
 
         default:
             break;

Implementierung/tricorn_V0.c

@@ -5,7 +5,7 @@ static unsigned converge(float c_real, float c_imag, unsigned n);
 
 
 /**
- * @brief Generates a Tricorn fractal image using highly optimized SSE/SIMD instructions.
+ * @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.

Implementierung/tricorn_V1.c

@@ -1,10 +1,11 @@
 #include "tricorn_V1.h"
+#include <immintrin.h>
 
-static unsigned converge(float c_real, float c_imag, unsigned n);
+static unsigned converge(float c_rea, float c_imag, unsigned n);
 
 
 /**
- * @brief Generates a Tricorn fractal image using an optimized SISD implementation.
+ * @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.
@@ -19,45 +20,94 @@ void tricorn_V1(float complex start, size_t width, ssize_t height, float res, un
     // 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);
 
-    float c_imag = c_start_imag;
+    __m128 four = _mm_set1_ps(4.0f);
-    size_t position = 0;
+    __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 = _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);
+            }
             
 
-    // Simple iteration over every pixel in the picture
+            __m128i mod64 = _mm_and_si128(iterations, _mm_set1_epi32(63));
-    for(ssize_t y = 0; y < actual_height; y++){
+            __m128i times4 = _mm_slli_epi32(mod64, 2);
+            __m128i values = _mm_sub_epi32(_mm_set1_epi32(255), times4);
 
-        // Calculate imaginary part on each row as it stays the same
+            __m128i max_iter = _mm_set1_epi32(n);
+            __m128i is_max = _mm_cmpeq_epi32(iterations, max_iter);
+            values = _mm_andnot_si128(is_max, values);
 
-        float c_real = c_start_real;
+            unsigned int value_array[4];
+            _mm_storeu_si128((__m128i*)value_array,values);
 
-        for(size_t x = 0; x < width; x++){
+            for(unsigned i = 0 ; i < 4;i++) {
-            unsigned value = converge(c_real,c_imag,n);
+                unsigned int value = value_array[i];
+                size_t position = (y*width + (x+i)) * bytesPerColor;
 
-            // Pixel value position in the output array
+                if(!color) {
-
-            if(!color){ // Grayscale colors as specified by the assignment r = g = b
                     img[position] = value;
                 }
-            else{ // Temporary colors with random multipliers
+                else {
                     img[position] = ((value*3) % 255);
-                img[position + 1] = ((value*2) % 255) ;
+                    img[position + 1] = ((value*2) % 255);
                     img[position + 2] = ((value*5) % 255);
                 }
-
-            position += bytesPerColor;
-            c_real += res;
-
             }
-        c_imag -= imag_step * 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);
+            }
+        }
+
     }
 }
 
-// Iterative application of the equation given by the assignment
+static unsigned converge(float c_real, float c_imag, unsigned n) {
-static inline unsigned converge(float c_real, float c_imag, unsigned n) {
     float real = 0.0f;
     float imag = 0.0f;
 

Implementierung/tricorn_V1.h

@@ -2,4 +2,5 @@
 #include <stddef.h>
 #include <math.h>
 #include <sys/types.h>
-void tricorn_V1(float complex  start, size_t width, ssize_t  height, float res, unsigned n, _Bool color, unsigned char* img);
+#include<stdbool.h>
+void tricorn_V1(float complex  start, size_t width, ssize_t  height, float res, unsigned n, bool color, unsigned char* img);

Implementierung/tricorn_V2.c

@@ -1,10 +1,10 @@
 #include "tricorn_V2.h"
 
-static const float convergenceLimit = 2.0f; // Mandelbrot set limit
+static unsigned converge(float c_real, float c_imag, unsigned n);
-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 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.
@@ -14,20 +14,30 @@ 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_V2(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;
-    unsigned bytesPerColor = color ? 3 : 1;
+    ssize_t imag_step = height < 0 ? -1:1;
+
+    float c_start_real = creal(start);
+    float c_start_imag = cimag(start);
+
+    float c_imag = c_start_imag;
+    size_t position = 0;
+
 
     // Simple iteration over every pixel in the picture
-    for(size_t i = 0; i < width; i++){
+    for(ssize_t y = 0; y < actual_height; y++){
-        for(ssize_t j = 0; j < height; j++){
 
-            unsigned value = converge(start + (i*res - I*j*res),n);
+        // Calculate imaginary part on each row as it stays the same
+
+        float c_real = c_start_real;
+
+        for(size_t x = 0; x < width; x++){
+            unsigned value = converge(c_real,c_imag,n);
 
             // Pixel value position in the output array
-            size_t position = (j*width +i) * bytesPerColor;
-
 
             if(!color){ // Grayscale colors as specified by the assignment r = g = b
                 img[position] =  value;
@@ -38,24 +48,29 @@ void tricorn_V2(float complex start, size_t width, ssize_t height, float res, un
                 img[position + 2] = ((value*5) % 255);
             }
 
+            position += bytesPerColor;
+            c_real += res;
+
         }
+        c_imag -= imag_step * res;
     }
 }
 
 // Iterative application of the equation given by the assignment
-static unsigned converge(float complex position,unsigned n){
+static inline 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++) {
-    for(unsigned i = 0; i < n; i++){
+        float real2 = real * real;
-         // Detection of unbounded behaviour
+        float imag2 = imag * imag;
-        // 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);
 
-        // Application of the given complex equation with m = 2
+        if (real2 + imag2 > 4.0f) {
-        result = cpowf(conjf(result),2) + position;
+            return 255 - 4 * (i % 64);
+        }
 
+        imag = -2.0f * real * imag + c_imag;
+        real = real2 - imag2 + c_real;
     }
     return 0;
 }

Implementierung/tricorn_V3.c

@@ -0,0 +1,61 @@
+#include "tricorn_V3.h"
+
+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).
+ * @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_V3(float complex start, size_t width, ssize_t height, float res, unsigned n, _Bool color, unsigned char* img){
+    
+    // Absolute height for iteration purposes
+    height = height < 0 ?  height*-1 : height;
+    unsigned bytesPerColor = color ? 3 : 1;
+
+    // Simple iteration over every pixel in the picture
+    for(size_t i = 0; i < width; i++){
+        for(ssize_t j = 0; j < height; j++){
+
+            unsigned value = converge(start + (i*res - I*j*res),n);
+
+            // Pixel value position in the output array
+            size_t position = (j*width +i) * bytesPerColor;
+
+
+            if(!color){ // Grayscale colors as specified by the assignment r = g = b
+                img[position] =  value;
+            }
+            else{ // Temporary colors with random multipliers
+                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){
+
+    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)
+            return 255-4*(i % 64);
+
+        // Application of the given complex equation with m = 2
+        result = cpowf(conjf(result),2) + position;
+       
+    }
+    return 0;
+}

Implementierung/tricorn_V3.h

@@ -0,0 +1,5 @@
+#include <complex.h>
+#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);