Wrong approximation of PI with concurrent program

I’m trying to solve this concurrent programming problem:

The points on a unit circle centered at the origin are defined by the function f(x) = sqrt(1-x²). Recall that the area of a circle is pi*r², where r is the radius. The adaptive quadrature routine described in Lecture 1 can approximate the pi value, computing the area of the upper-right quadrant of a unit circle and then multiplying the result by 4. Develop a multithreaded program (in C using Pthreads or in Java) to compute pi for a given epsilon using a given number of processes (threads) np which is assumed to be a command-line argument (i.e., it’s given).

However, I’m having trouble getting the correct approximation of pi (I get 2.666…) when using increasingly small trapezoids to cover the area described and dividing the problem successively with processors.

Note: in the code, worker refers to a processor.

Here’s my code:

#ifndef _REENTRANT
#define _REENTRANT
#endif

#include <pthread.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdbool.h>
#include <time.h>
#include <sys/time.h>
#include <math.h>

#define MAXWORKERS 3        // Maximum number of additional allowed workers besides the main thread.

const double EPSILON = 10e-10;
int numWorkers = 0;
double startTime, endTime;
int numCreatedThreads = 1;

// Mutex lock for numCreatedThreads.
pthread_mutex_t createdThreadsLock;

struct Info {
    double a, b;     // Extreme points.
    double fa, fb;   // Evaluation of extreme points.
    double area;     // Area to be approximated.
};

double read_timer() {
    static bool initialized = false;
    static struct timeval start;
    struct timeval end;
    if(!initialized) {
        gettimeofday(&start, NULL);
        initialized = true;
    }
    gettimeofday(&end, NULL);
    return (end.tv_sec - start.tv_sec) + 1.0e-6 * (end.tv_usec - start.tv_usec);
}

void readCommandLine(int argc, char* argv[]) {
    numWorkers = (argc > 1)? atoi(argv[1]) : MAXWORKERS;
    if (numWorkers > MAXWORKERS || numWorkers < 1)
        numWorkers = MAXWORKERS;
}

double f(double x) {
    return 1 - x*x;
}

void displayResults(double piApprox) {
    printf("\n===========================RESULTS===========================\n");
    printf("Pi was approximated up to %.0f decimal places: %f\n", -log10(EPSILON), 4 * piApprox);
    printf("The execution time is %g seconds.\n", endTime - startTime);
    printf("=============================================================\n");
}

void * calculatePI(void *args) {

    struct Info *info = args;

    // Calculate new data.
    double m = (info->a + info->b) / 2;
    double fm = f(m);
    double larea = (info->fa + fm) * (m - info->a) / 2;
    double rarea = (fm + info->fb) * (info->b - m) / 2;

    // Final result that will be returned.
    double *res = malloc(sizeof(double));

    // Check for termination.
    if(fabs(larea + rarea - info->area) < EPSILON) {
        *res = larea + rarea;
        return res;
    }

    // Boolean to control number of threads working.
    bool tooManyThreads = false;
    void* resultLeft, *resultRight;
    struct Info newLeft = {info->a, m, f(info->a), f(m), larea};
    struct Info newRight = {m, info->b, f(m), f(info->b), rarea};

    // Check whether we surpass the number of allowed threads.
    pthread_mutex_lock(&createdThreadsLock);
        tooManyThreads = numCreatedThreads + 1 > MAXWORKERS;
    pthread_mutex_unlock(&createdThreadsLock);

    if(!tooManyThreads) {

        // Update numCreatedThreads atomically.
        pthread_mutex_lock(&createdThreadsLock);
            numCreatedThreads++;
        pthread_mutex_unlock(&createdThreadsLock);

        // Execute left side with new thread, and right side with current thread.
        pthread_t newThread;
        pthread_create(&newThread, NULL, calculatePI, &newLeft);
        resultRight = calculatePI(&newRight);
        pthread_join(newThread, &resultLeft);
    }

    else {
        // Calculate the area of each side recursively on the same thread.
        resultLeft = calculatePI(&newLeft);
        resultRight = calculatePI(&newRight);
    }

    double * resultL = resultLeft;
    double * resultR = resultRight;
    *res = *resultL + *resultR;
    return res;
}

int main(int argc, char* argv[]) {

    // Read command line args if any.
    readCommandLine(argc, argv);

    // Initialize createdThreads mutex;
    pthread_mutex_init(&createdThreadsLock, NULL);

    // Initialize first information.
    struct Info info = {
            .a = 0,
            .b = 1,
            .fa = f(0),
            .fb = f(1),
            .area = (f(0) + f(1)) / 2
    };

    startTime = read_timer();

    // Obtain approximation of pi/4
    double * piApproximation = (double *) calculatePI(&info);

    endTime = read_timer();
    displayResults(*piApproximation);

    return 0;
}

What am I missing?

Thanks!

Your integrand f is wrong. For unit circle, we have x*x + y*y == 1, your f instead calculates y = 1-x*x, which means what you are actually caluclating is 4 times the integral of (1-x*x) from x=0 to x=1. The integral can be analytically evaluated to be g(1)-g(0), where g(x) = x - x*x*x/3. That’s why you get 4*(1-1/3)=8/3=2.6666.

You need to fix it by integrating the correct function. Since you are integrating the first quadrant, x*x + y*y == 1 transforms to y=sqrt(1-x*x), so you should modify f to return it.

From the comment by @chux, A numerically stable way to implement it is to use (1-x)*(1+x) to calculate (1-x*x). In general, one should avoid addition and subtraction between two floating point numbers which differ much in terms of orders of magnitudes. For example, if x is just a bit larger than sqrt(DBL_EPSILON), most of the significant digits in x*x will be lost when calculating 1-x*x. Instead, in (1-x)*(1+x) most of the significant digits are preserved in the two factors before they are multiplied together. So you should modify f to the following instead,

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