Multithreaded π calculation/estimation

pi.c

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

/* Estimate pi using the Monte Carlo algorithm */
double pi_thread(int n) {
    int n1 = 0;
    double x, y;
    unsigned int seed;

    #pragma omp parallel reduction(+:n1) private(x,y,seed)
    {
        seed = omp_get_thread_num();

        #pragma omp for
        for (int i = 0; i < n; ++i)
        {
            x = (double)rand_r(&seed)/RAND_MAX;
            y = (double)rand_r(&seed)/RAND_MAX;
            if ((x*x + y*y) <= 1.0) n1++;
        }
    }

    return 4.0*n1/n;
}

double f(double x) {
    return sqrt(1.0 - x*x);
}

/* Estimate pi using the integral */
double pi_thread_area(int n) {
    double step = 1.0/n;
    double area = 0.0;

    #pragma omp parallel for reduction(+:area)
    for (int i = 0; i < n; i++)
        area += f(i*step) * step;

    return 4.0*area;
}

int main(int argc, char const *argv[])
{
    if (argc < 2) {
        fprintf(stderr, "Usage: %s <num_exp_total>\n", argv[0]);
        return 1;
    }

    /* Number of iterations */
    int n = atoi(argv[1]);

    printf("Result: %.9f\n", pi_thread_area(n));
    printf("Result: %.9f\n", pi_thread(n));

    return 0;
}