Parallel C code for finding the min and max bounds of a function of the form float f(float,float)

getBounds.c

#include <stdio.h>
#include <pthread.h>
#include <float.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

#include "xoshiro128plus.c"
#include "splitmix64.c"

#define NUM_SAMPLES 1000
size_t num_samples = NUM_SAMPLES;

float test_function(float x, float y)
{
     return x*y;
}

void init_random_number(void* state, size_t bytes)
{
    for( size_t i = 0, n = bytes/sizeof(uint64_t); i < n; ++i )
    {
       *((uint64_t*)((uint32_t*)state+i*2)) = splitmix64_next();
    }
}

float random_number( uint32_t* state )
{
    int32_t u = next(state)>>4u;
    return (float)u / (float)(1u<<28u);
}

// find min and max bounds on all possible input for function with two inputs
void* get_function_bounds( void* state )
{
    float* ret = (float*)malloc(6 * sizeof(float));
    float min_x1,min_x2,min_y;
    float max_x1,max_x2,max_y;
    min_x1 = min_x2 = min_y = FLT_MAX;
    max_x1 = max_x2 = max_y = -FLT_MAX;
    for( size_t i = 0; i < num_samples ; i++ )
    {
        // check random inputs
        float x1 = random_number(state);
        float x2 = random_number(state);
        float y = test_function(x1,x2);
        if( x1 < min_x1 ) min_x1 = x1;
        if( x1 > max_x1 ) max_x1 = x1;
        if( x2 < min_x2 ) min_x2 = x2;
        if( x2 > max_x2 ) max_x2 = x2;
        if( y < min_y ) min_y = y;
        if( y > max_y ) max_y = y;
    }

    ret[0] = min_x1;
    ret[1] = max_x1;
    ret[2] = min_x2;
    ret[3] = max_x2;
    ret[4] = min_y;
    ret[5] = max_y;

    return ret;
}

// usage: getBounds <number of samples> [<thread count>]
int main(int argc, char** argv)
{
    num_samples = NUM_SAMPLES;
    size_t num_threads = 1;
    if( argc > 1 )
        num_samples = atoi(argv[1]);
    if( argc > 2 )
        num_threads = atoi(argv[2]);

    splitmix64_x = time(NULL);
    
    // run get_function_bounds on each thread
    pthread_t* threads = (pthread_t*)malloc(num_threads * sizeof(pthread_t));
    float* returnBuffer = NULL;
    float* bounds = (float*)malloc(num_threads * 6 * sizeof(float));
    uint32_t* states = (uint32_t*)malloc(num_threads * sizeof(uint32_t) * 4u);
    for( size_t i = 0; i < num_threads; i++ )
    {
        if( i == 0 ) init_random_number( states+i*4, sizeof(uint32_t)*4u );
        else
        {
            memcpy(states+i,states+i-1,sizeof(uint32_t)*4u);
            jump(states+i);
        }
        pthread_create(&threads[i], NULL, get_function_bounds, (void*)states+i*4);
    }
    for( size_t i = 0; i < num_threads; i++ )
    {
        pthread_join(threads[i], (void**)&returnBuffer);
        memcpy(bounds + (i*6), returnBuffer, 6 * sizeof(float));
        free(returnBuffer);
    }
    
    free(threads);
    free(states);
    
    // find min and max bounds on all possible input for function with two inputs
    float min_x1 = FLT_MAX;
    float max_x1 = -FLT_MAX;
    float min_x2 = FLT_MAX;
    float max_x2 = -FLT_MAX;
    float min_y = FLT_MAX;
    float max_y = -FLT_MAX;
    for( size_t i = 0; i < num_threads; i++ )
    {
        if( bounds[i*6] < min_x1 ) min_x1 = bounds[i*6];
        if( bounds[i*6+1] > max_x1 ) max_x1 = bounds[i*6+1];
        if( bounds[i*6+2] < min_x2 ) min_x2 = bounds[i*6+2];
        if( bounds[i*6+3] > max_x2 ) max_x2 = bounds[i*6+3];
        if( bounds[i*6+4] < min_y ) min_y = bounds[i*6+4];
        if( bounds[i*6+5] > max_y ) max_y = bounds[i*6+5];
    }
    printf("y = f(x1,x2)\n-------------------\n");
    printf("min_x1: %f\n", min_x1);
    printf("max_x1: %f\n", max_x1);
    printf("min_x2: %f\n", min_x2);
    printf("max_x2: %f\n", max_x2);
    printf("min_y:  %f\n", min_y);
    printf("max_y:  %f\n", max_y);
    
    free(bounds);
    return 0;
}