GRA-Project/Implementierung/tricorn_V4.c

#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 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_V4(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;
    _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;

    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
}