calcular integrales dobles en R rápidamente

Estoy buscando una solución para una integral doble que sea más rápida que

integrate(function(y) { sapply(y, function(y) { integrate(function(x) myfun(x,y), llim, ulim)$value }) }, llim, ulim) 

con eg

 myfun <- function(x,y) cos(x+y) llim <- -0.5 ulim <- 0.5 

Encontré un documento antiguo que se refería a un progtwig de FORTRAN llamado quad2d , pero no pude encontrar nada más que algunas páginas de ayuda para matlab para el rest. Así que estoy buscando una biblioteca C o FORTRAN que pueda hacer integrales dobles rápidamente (es decir, sin el bucle sapply), y que se puede llamar desde R. Todas las ideas son muy apreciadas, siempre y cuando todas sean compatibles con GPL.

Si la solución implica llamar a otras funciones de las bibliotecas que ya vienen con R, también me gustaría saber de ellas.

El paquete de cubature hace la integración 2D (y ND) usando un algoritmo adaptativo. Debería superar los enfoques más simples para la mayoría de los integrandos.

El paquete de pracma que Joshua señaló contiene una versión de quad2d .

 quad2d(myfun, llim, ulim, llim, ulim) 

Esto le da la misma respuesta, dentro de la tolerancia numérica, como su bucle, utilizando la función de ejemplo.

Por defecto, con su función de ejemplo, quad2d es más lento que el bucle. Si bajas n , puedes hacerlo más rápido, pero supongo que depende de cuán suave sea tu función y de cuánta precisión estés dispuesto a sacrificar por velocidad.