The Madhava-Leibniz approximation to Pi

Pthreads versions, using a busy-wait loop and using a mutex to protect the critical section.

Madhava-Leibniz-pthread.c

Basic pthreads version, with no protection against a race condition when accessing the shared global sum variable.

/*
| Pi approximation, based on the Madhava-Leibnitz approximation
| 2020-10-12
*/
#include <stdio.h>
#include <stdlib.h> // strtod(), atoi()
#include "timing.h" // TIMING_initTime(), TIMING_reporttime()
#include <math.h>   // log10()

#include <unistd.h> // gethostname()
#define MAXNAME 256 // should be enough...

#include <pthread.h>

#define LABEL "Pthreads"

struct Params {
    long double start, stop;
    int rank;
} ;


// long double only suppports 35 digits
#define PI35 3.1415926535897932384626433832795029L
//           -.---|----|----|----|----|----|----|

unsigned nthreads;
long double sum = 0.0L;

void *terms(void *vp)
{
    struct Params *params = vp;
    long double factor = ((unsigned long)(params->start - 1) % 4) ? -1.0L : 1.0L ;
    for (long double x = params->start; x < params->stop; x += 2) {

        sum += factor * 1.0L / x;

        factor = -factor;
    }
    return NULL;
}


int main(int argc, char **argv)
{
    char hostname[MAXNAME];
    gethostname(hostname, MAXNAME);
    hostname[MAXNAME-1] = '\0';     // guarantee a terminating null character

    TIMING_initTime(0, hostname, 0);

    // printf("Starting...\n");

    long double nterms;
    if (argc == 3) {
        nthreads = atoi(argv[1]);
        nterms = strtold(argv[2], NULL);
    } else {
        fprintf(stderr, "usage: %s <nthreads> <nterms>\n", argv[0]);
        return 1;
    }

    pthread_t threads[ nthreads ];
    struct Params arguments[ nthreads ];

    //long double first_term = 1.0L, last_term = 2.0 * nterms;
    int terms_per_thread = nterms / nthreads;   // what if nthreads not a factor of nterms ?
    for (int n = 0; n < nthreads; n++) {
        arguments[n].start = n * (2*terms_per_thread) + 1;
        arguments[n] .stop = (n+1) * (2*terms_per_thread) + 1;
        arguments[n].rank = n;

        printf("thread %d:  start %Lf,  stop %Lf\n",
            arguments[n].rank, arguments[n].start, arguments[n].stop);
    }

    // This main() IS thread 0...
    for (int n = 1; n < nthreads; n++) {
        //pthread_create( &threads[n], NULL, terms, &arguments[n] );
        pthread_create( (threads+n), NULL, terms, (arguments+n) );
    }

    //------------------------------------------------------------------
    // Sum up the terms:
    terms( (arguments+0) );

    for (int n = 1; n < nthreads; n++) {
        pthread_join( *(threads+n), NULL);
    }

    // Calculate the approximation:
    long double pi_approx = 4.0L * sum;

    //------------------------------------------------------------------
    // Report the calculation results:
    long double error = pi_approx - PI35;
    if (error < 0.0L)
        error = -error;
    int width = abs(log10(error)) + 2;  // show this many digits...
    printf("#---- %s -- %Lg terms, %i threads ----\n",
        LABEL, nterms, nthreads);
    printf("pi35bits: %.*Lf\n  approx: %.*Lf\nerror=%.6Le\n#--------\n",
        width, PI35,  width, pi_approx,  error);

    TIMING_reportTime("Serial run", 0);

    return 0;
}

Madhava-Leibniz-busywait.c

Add "busy-wait" protection to the critical section that accesses the shared global sum variable.

/*
| Pi approximation, based on the Madhava-Leibnitz approximation
| 2020-10-12
*/
#include <stdio.h>
#include <stdlib.h> // strtod(), atoi()
#include "timing.h" // TIMING_initTime(), TIMING_reporttime()
#include <math.h>   // log10()

#include <unistd.h> // gethostname()
#define MAXNAME 256 // should be enough...

#include <pthread.h>


#define LABEL "Busywait"

struct Params {
    long double start, stop;
    int rank;
} ;


// long double only suppports 35 digits
#define PI35 3.1415926535897932384626433832795029L
//           -.---|----|----|----|----|----|----|

unsigned nthreads;
long double sum = 0.0L;
int wait_flag = 0;

void *terms(void *vp)
{
    struct Params *params = vp;
    long double factor = ((unsigned long)(params->start - 1) % 4) ? -1.0L : 1.0L ;
    for (long double x = params->start; x < params->stop; x += 2) {

        while (wait_flag != params->rank) {
            ;
        }
        sum += factor * 1.0L / x;
        wait_flag = (wait_flag + 1) % nthreads;

        factor = -factor;
    }
    return NULL;
}


int main(int argc, char **argv)
{
    char hostname[MAXNAME];
    gethostname(hostname, MAXNAME);
    hostname[MAXNAME-1] = '\0';     // guarantee a terminating null character

    TIMING_initTime(0, hostname, 0);

    // printf("Starting...\n");

    long double nterms;
    if (argc == 3) {
        nthreads = atoi(argv[1]);
        nterms = strtold(argv[2], NULL);
    } else {
        fprintf(stderr, "usage: %s <nthreads> <nterms>\n", argv[0]);
        return 1;
    }

    pthread_t threads[ nthreads ];
    struct Params arguments[ nthreads ];
    
    //long double first_term = 1.0L, last_term = 2.0 * nterms;
    int terms_per_thread = nterms / nthreads;   // what if nthreads not a factor of nterms ?
    for (int n = 0; n < nthreads; n++) {
        arguments[n].start = n * (2*terms_per_thread) + 1;
        arguments[n] .stop = (n+1) * (2*terms_per_thread) + 1;
        arguments[n].rank = n;

//        printf("thread %d:  start %Lf,  stop %Lf\n",
//            arguments[n].rank, arguments[n].start, arguments[n].stop);
    }

    // This main() IS thread 0...
    for (int n = 1; n < nthreads; n++) {
        //pthread_create( &threads[n], NULL, terms, &arguments[n] );
        pthread_create( (threads+n), NULL, terms, (arguments+n) );
    }

    //------------------------------------------------------------------
    // Sum up the terms:
    terms( (arguments+0) );

    for (int n = 1; n < nthreads; n++) {
        pthread_join( *(threads+n), NULL);
    }

    // Calculate the approximation:
    long double pi_approx = 4.0L * sum;

    //------------------------------------------------------------------
    // Report the calculation results:
    long double error = pi_approx - PI35;
    if (error < 0.0L)
        error = -error;
    int width = abs(log10(error)) + 2;  // show this many digits...
    printf("#---- %s -- %Lg terms, %i threads ----\n",
        LABEL, nterms, nthreads);
    printf("pi35bits: %.*Lf\n  approx: %.*Lf\nerror=%.6Le\n#--------\n",
        width, PI35,  width, pi_approx,  error);

    TIMING_reportTime("Serial run", 0);

    return 0;
}

Madhava-Leibniz-mutex.c

Replace the "busy-wait" protection with a "mutex" global object.

/*
| Pi approximation, based on the Madhava-Leibnitz approximation
| 2020-10-12
*/
#include <stdio.h>
#include <stdlib.h> // strtod(), atoi()
#include "timing.h" // TIMING_initTime(), TIMING_reporttime()
#include <math.h>   // log10()

#include <unistd.h> // gethostname()
#define MAXNAME 256 // should be enough...

#include <pthread.h>


#define LABEL "Mutex"

struct Params {
    long double start, stop;
    int rank;
} ;


// long double only suppports 35 digits
#define PI35 3.1415926535897932384626433832795029L
//           -.---|----|----|----|----|----|----|

unsigned nthreads;
long double sum = 0.0L;
int wait_flag = 0;
pthread_mutex_t mutex;

void *terms(void *vp)
{
    struct Params *params = vp;
    long double factor = ((unsigned long)(params->start - 1) % 4) ? -1.0L : 1.0L ;
    for (long double x = params->start; x < params->stop; x += 2) {
//        while (wait_flag != params->rank) {
//            ;
//        }
        pthread_mutex_lock( &mutex );
        sum += factor * 1.0L / x;
        pthread_mutex_unlock( &mutex );
//        wait_flag = (wait_flag + 1) % nthreads;

        factor = -factor;
    }
    return NULL;
}


int main(int argc, char **argv)
{
    char hostname[MAXNAME];
    gethostname(hostname, MAXNAME);
    hostname[MAXNAME-1] = '\0';     // guarantee a terminating null character

    TIMING_initTime(0, hostname, 0);

    // printf("Starting...\n");

    long double nterms;
    if (argc == 3) {
        nthreads = atoi(argv[1]);
        nterms = strtold(argv[2], NULL);
    } else {
        fprintf(stderr, "usage: %s <nthreads> <nterms>\n", argv[0]);
        return 1;
    }

    pthread_t threads[ nthreads ];
    struct Params arguments[ nthreads ];

    //long double first_term = 1.0L, last_term = 2.0 * nterms;
    int terms_per_thread = nterms / nthreads;   // what if nthreads not a factor of nterms ?
    for (int n = 0; n < nthreads; n++) {
        arguments[n].start = n * terms_per_thread + 1;
        arguments[n] .stop = (n+1) * terms_per_thread + 1;
        arguments[n].rank = n;

//        printf("thread %d:  start %Lf,  stop %Lf\n",
//            arguments[n].rank, arguments[n].start, arguments[n].stop);
    }

    pthread_mutex_init(&mutex, NULL);

    // This main() IS thread 0...
    for (int n = 1; n < nthreads; n++) {
        //pthread_create( &threads[n], NULL, terms, &arguments[n] );
        pthread_create( (threads+n), NULL, terms, (arguments+n) );
    }

    //------------------------------------------------------------------
    // Sum up the terms:
    terms( (arguments+0) );

    for (int n = 1; n < nthreads; n++) {
        pthread_join( *(threads+n), NULL);
    }

    pthread_mutex_destroy(&mutex);

    // Calculate the approximation:
    long double pi_approx = 4.0L * sum;

    //------------------------------------------------------------------
    // Report the calculation results:
    long double error = pi_approx - PI35;
    if (error < 0.0L)
        error = -error;
    int width = abs(log10(error)) + 2;  // show this many digits...
    printf("#---- %s -- %Lg terms, %i threads ----\n",
        LABEL, nterms, nthreads);
    printf("pi35bits: %.*Lf\n  approx: %.*Lf\nerror=%.6Le\n#--------\n",
        width, PI35,  width, pi_approx,  error);

    TIMING_reportTime("Serial run", 0);

    return 0;
}