Q:

Problem with NIntegrate over a highly-oscillatory integrand

I'm trying to numerically evaluate the integral $$\int_{a}^{b}\mathop{\mathrm{d}x}\int_{x}^{b}\frac{\sin(x-y)}{xy}\mathop{\mathrm{d}y}$$
using Mathematica. To do that, I the function
Si2[a_, b_] := NIntegrate[Sin[x - y]/(x y), {x, a, b}, {y, x, b},
  AccuracyGoal -> 25, PrecisionGoal -> 25, WorkingPrecision -> 40,
  MaxRecursion -> 1000000, Method -> "InterpolationPointsSubdivision"];

However, running Si2[.1,1] gives the error

NIntegrate::inumr: The integrand Sin[x-y]/(x y) has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1},{0,1}}.

However, I'm not integrating over $x=y=0$ (which is an obvious singularity). 
Two notes:

This is an example of usage. In practice, I need to evaluate this function for parameters much closer to $a=0$ (e.g. $\log_{10}(a)\sim-6$).
Note that I use the InterpolationPointsSubdivision method because I saw in various answers that it is a good method to evaluate numerically a highly-oscillatory integrand. I tried to use few other methods, but got the same error.

Any advice? Thanks!

A:

Here is @JM's answer. First integrate the interior integral:
i1[a_, b_] = Integrate[Sin[x-y]/y, {y, x, b}, Assumptions -> 0<x<b]

(CosIntegral[b] - CosIntegral[x]) Sin[x] + 
   Cos[x] (-SinIntegral[b] + SinIntegral[x])

Then, use this integral:
int[a_, b_, opts:OptionsPattern[NIntegrate]] := NIntegrate[i1[a, b]/x, {x, a, b}, opts]

Reproducing previous results:
int[.1, 1, WorkingPrecision->40] //AbsoluteTiming
int[10^-6, 1, WorkingPrecision->40] //AbsoluteTiming

{0.413422, -0.7031515701355189970289164059287652664106}
{0.882268, -11.18694435428611911727667338594134335609}

A:

Something like that? 
Si2[a_, b_] := 
 NIntegrate[Sin[x - y]/(x y), {x, a, b}, {y, x, b}, 
   AccuracyGoal -> 25, PrecisionGoal -> 25, WorkingPrecision -> 40, 
   MaxRecursion -> 100000000, {Method -> "QuasiMonteCarlo"}] // Quiet
Si2[0.1, 1]

The result is 
-0.7030662536921237781417087351943610040114
and it's not very slow 
In[65]:= % // AbsoluteTiming

Out[65]= {0.007766, -0.7030662536921237781417087351943610040114}

so you can try to increase maxrecursion, etc etc

