#include "tricorn.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(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; 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; const int bytes_per_color = args->color ? 3 : 1; // Absolute height for iteration purposes const 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 }