# Optimal way of performing a high-dimensional numerical integral of a sharped function

I am trying to solve the following numerical integral

fun=E^(1/48 (-113 k[1]^2+102 Cos[x1-x3] k[1] k[3]-113 k[3]^2-2 Cos[x1-x5] k[1] k[5]+102 Cos[x3-x5] k[3] k[5]-113 k[5]^2+102 Cos[x1-x7] k[1] k[7]-2 Cos[x3-x7] k[3] k[7]+102 Cos[x5-x7] k[5] k[7]-113 k[7]^2))/(36864 \[Pi]^4);
Timing[NIntegrate[k[1]k[3]k[5]k[7]fun Cos[2(x1+x3-x5-x7)],{k[1],0,Infinity},{k[3],0,Infinity},{k[5],0,Infinity},{k[7],0,Infinity},{x1,0,2Pi},{x3,0,2Pi},{x5,0,2Pi},{x7,0,2Pi}]]


which is basically a 8-dimensional gaussian function times a cosine.

The result that I'm obtaining is

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.
NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 5.956319478409388*^-7 and 8.378969793277401*^-7 for the integral and error estimates.
{43.3281,5.95632*10^-7}


My problem is that the error estimates of the integral is bigger than its value which makes me doubt about the accuracy of the result and also I would like to minimize the timing of computation.

Is there a method of computing this integral in less than 40 seconds and obtaining a more or less accurate result?

I don't mind if the result is $$6 \cdot 10^{-7}$$ or $$5 \cdot 10^{-7}$$, I just want an approximate value in the less time possible.

Timing[NIntegrate[