# Replacing Mathematica's QuasiMonteCarlo integration in C++

I have a Mathematica program which performs some integrals in 3 or 4 dimensions using the QuasiMonteCarlo method. The problem is, it takes an annoyingly long time to run, to the point where some of these calculations can't complete in the maximum job time available on our HPC cluster. So I'm considering rewriting the program in C++, which I suspect will speed it up by a large factor.

I looked at the GSL docs and while there are sections on quasirandom sequences and regular MC integration, I don't see anything that brings them together. Also a Google search or two didn't turn up anything that looked like a widely trusted implementation. What are my options for a well-tested implementation of QMC integration in C++?

In the interest of consistency, I'd prefer to use something close to the Halton-Hammersley-Wozniakowski method that Mathematica implements, if that is an option.

• you could post the integral on Stackoverflow and we could have a look at what is going on. Note, the Mathematica is based on MKL for machine precision, which is quite efficient.
– user530
Jan 14 '12 at 7:41
• Not an answer to your question, but I was wondering if you tried Compileing the integral (to C code) before passing it to NIntegrate, i.e. is it NIntegrate that's slow or calculating the function? Using C-compiled functions might need a bit extra work on a cluster though. Jan 14 '12 at 8:39
• Compiling sounds like a very good idea, I hadn't thought of that. I'll give it a try. Each of these calculations evaluates the function about 5 million times, by my estimate, and since the entire computation is taking about 3 hours, that's 2ms per function evaluation, which seems rather slow for some purely numerical calculations. Jan 15 '12 at 0:00
• @ruebenko: I'll keep that in mind. Jan 15 '12 at 0:00
• The CUBA Library has a number of algorithms for low-dimensional problems. It even has a Mathematica interface. feynarts.de/cuba
– dls
Mar 20 '12 at 19:23