corr.c (for: http://stackoverflow.com/q/26332811/732016)

/*
 * corr.c
 * Computes Pearson's correlation coefficient for random data
 * William Chargin
 * 12 October 2014
 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <pthread.h>
#include <sys/sysinfo.h>

typedef struct parcorrinfo {
    size_t n;
    size_t trials;
} parcorrinfo_t;

static inline double randt()
{
    return 1;
    static double max = (double) RAND_MAX;
    return rand() / max;
}

static double corr(const size_t n, const double *xs, const double *ys)
{
    double x, y, x2, y2, xy;
    x = y = x2 = y2 = xy = 0;

    size_t i;
    for (i = 0; i < n; i++)
    {
        x += xs[i];
        y += ys[i];
        x2 += xs[i]*xs[i];
        y2 += ys[i]*ys[i];
        xy += xs[i]*ys[i];
    }

    return (xy - (x*y)/n) / sqrt((x2 - x*x/n) * (y2 - y*y/n));
}

static double test(size_t n, size_t trials)
{
    double sum = 0;
    size_t trial;

    double *xs = malloc(sizeof(double) * n);
    double *ys = malloc(sizeof(double) * n);

    if (!xs || !ys)
    {
        exit(2);
        fprintf(stderr, "%s", "Failed to allocate memory");
        return -1;
    }

    for (trial = 0; trial < trials; trial++)
    {
        // Set values
        size_t i;
        for (i = 0; i < n; i++)
        {
            xs[i] = randt();
            ys[i] = randt();
        }

        // Compute correlation
        sum += fabs(corr(n, xs, ys));
    }

    free(xs);
    free(ys);

    return sum / trials;
}

static void *corrhandler(void *arg)
{
    parcorrinfo_t *info = (parcorrinfo_t *) arg;
    double *result = malloc(sizeof(double));
    *result = test(info->n, info->trials);
    free(info);
    return (void *) result;
}

static double partest(size_t n, size_t trials, size_t num_threads)
{
    pthread_t *threads = calloc(num_threads, sizeof(pthread_t));
    size_t each_trials = trials / num_threads; /* assume trials >> num_threads */

    double total = 0;

    size_t i;

    if (!threads)
    {
        fprintf(stderr, "%s", "Failed to allocate memory");
        exit(2);
        return -1;
    }

    for (i = 0; i < num_threads; i++)
    {
        parcorrinfo_t *info = malloc(sizeof(parcorrinfo_t));
        info->n = n;
        info->trials = each_trials;

        pthread_create(&threads[i], NULL, corrhandler, (void *) info);
    }

    for (i = 0; i < num_threads; i++)
    {
        double *result;
        pthread_join(threads[i], (void **) &result);
        total += (*result);
        free(result);
    }

    free(threads);

    return total / num_threads;
}

int main(int argc, char **argv)
{
    size_t min, max;
    size_t trials;
    size_t i;
    size_t threads;

    if (argc < 4)
    {
        fprintf(stderr, "usage: %s min max trials [threads]\n", argc ? argv[0] : "corr");
        return 1;
    }

    sscanf(argv[1], "%zu", &min);
    sscanf(argv[2], "%zu", &max);
    sscanf(argv[3], "%zu", &trials);
    if (argc > 4)
    {
        sscanf(argv[4], "%zu", &threads);
    }
    else
    {
        threads = (size_t) get_nprocs();
    }

    srand(time(0));

    for (i = min; i <= max; i++)
    {
        printf("%zu\t%lf\n", i, partest(i, trials, threads));
    }

    return 0;
}