#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>

/* Allow compiler -DTHREADCOUNT=4 but have a default */
#ifndef THREADCOUNT
#define THREADCOUNT 16
#endif

/** Defines a particular task to handle */
typedef struct {
    size_t from;
    size_t to;
    double (*getnum)(size_t);
} task_description;

/**
 * Function invoked by each new thread.
 * The argument must be a task_description *;
 * the return value is a malloced double *.
 */
void *sum_array(void *args) {
    task_description *task = (task_description *)args;
    printf("  Summing from %zd to %zd...\n", 
        task->from, task->to);
    double work = 0;
    for(size_t i=task->from; i<task->to; i+=1) {
        work += task->getnum(i);
    }
    printf("  ... sub-sum from %zd to %zd = %g\n", 
        task->from, task->to, work);
    double *sum = malloc(sizeof(double));
    *sum = work;
    return (void *)sum;
}

/** A simple reciprocal function */
double fraction(size_t i) {
    return 1.0/(i+1);
}

/**
 * Sum all fractions 1/n from 1 to pow(2,-28)
 * in THREADCOUNT parallel threads
 */
int main(int argc, const char *argv[]) {
    // set up task sizes to take a few seconds on 2019-era laptops
    size_t max = 1<<28;
    size_t step = max / THREADCOUNT;

    // store per-thread information (don't re-use, memory is shared)
    pthread_t id[THREADCOUNT];
    task_description tasks[THREADCOUNT];

    // spawn the threads
    for(int i=0; i<THREADCOUNT; i+=1) {
        tasks[i].from = i*step;
        tasks[i].to = (i+1)*step;
        tasks[i].getnum = fraction;
        pthread_create(id+i, NULL, sum_array, tasks+i);
    }

    // wait for and combine the results
    double result = 0;
    for(int i=0; i<THREADCOUNT; i+=1) {
        void *ans;
        pthread_join(id[i], &ans);
        result += *(double *)ans;
        free(ans); // was malloced in just-joined thread
    }

    printf("The sum is %g\n", result);
    return 0;
}