piWithInterrupts.cpp

Go to the documentation of this file.

// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 8 -*-
 
#include <Rcpp.h>
 
#ifdef _OPENMP
#include <omp.h>
#endif
 
#include <R_ext/Utils.h>
 
class interrupt_exception : public std::exception {
public:
    interrupt_exception(std::string message)
        : detailed_message(message)
    {};
 
    virtual ~interrupt_exception() throw() {};
 
    virtual const char* what() const throw() {
        return detailed_message.c_str();
    }
 
    std::string detailed_message;
};
 
static inline void check_interrupt_impl(void* /*dummy*/) {
    R_CheckUserInterrupt();
}
 
inline bool check_interrupt() {
    return (R_ToplevelExec(check_interrupt_impl, NULL) == FALSE);
}
 
RcppExport SEXP PiLeibniz(SEXP n, SEXP frequency)
{
    BEGIN_RCPP
 
    // cast parameters
    int n_cycles = Rcpp::as<int>(n);
    int interrupt_check_frequency = Rcpp::as<int>(frequency);
 
    // user interrupt flag
    bool interrupt = false;
 
    double pi = 0;
#ifdef _OPENMP
#pragma omp parallel for \
    shared(interrupt_check_frequency, n_cycles, interrupt)      \
    reduction(+:pi)
#endif
    for (int i=0; i<n_cycles; i+=interrupt_check_frequency) {
        // check for user interrupt
        if (interrupt) {
            continue;
        }
 
#ifdef _OPENMP
        if (omp_get_thread_num() == 0) // only in master thread!
#endif
            if (check_interrupt()) {
                interrupt = true;
            }
 
        // do actual computations
        int j_end = std::min(i+interrupt_check_frequency, n_cycles);
        for (int j=i; j<j_end; ++j) {
            double summand = 1.0 / (double)(2*j + 1);
            if (j % 2 == 0) {
                pi += summand;
            }
            else {
                pi -= summand;
            }
        }
    }
 
    // additional check, in case frequency was too large
    if (check_interrupt()) {
        interrupt = true;
    }
 
    // throw exception if interrupt occurred
    if (interrupt) {
        throw interrupt_exception("The computation of pi was interrupted.");
    }
 
    pi *= 4.0;
 
    // result list
    return Rcpp::wrap(pi);
 
    END_RCPP
}