the winding number for the circle map (Arnold tongue)

I have figured out why you are getting the structure you are getting. The reason has to do with your initial choice of the angle, which you set at $0$ in the Nest[] statement. The actual image is generated by choosing the mean result of iterating the map for many initial values chosen uniformly at random in $[0,1]$.

With $n = 50$ iterations and $m = 20$ trials, I obtained the following image using the following modification of your code:

f = Compile[{n, m, a, k}, Mean[Table[Nest[# + a - k Sin[2 Pi #] &,
    RandomReal[], n]/n, {j, 1, m}]]];
dat = Parallelize[Outer[f[50, 20, #2, #] &, Range[0, 1, 1/1000], 
    Range[0, 1, 1/1000]]];

I believe this is very, very close to exactly what you are looking for.

This question looks as a duplicate of these questions:

How to create internally optimized expression for computing with high WorkingPrecision?

How to work with Experimental`NumericalFunction?

An internally optimized version of the original function can be created as follows:

n = 500;
f = Experimental`CreateNumericalFunction[{a, b}, 
   Unevaluated[Nest[# + a - b Sin[2 \[Pi] #] &, 0, n]/n], {}, 
   WorkingPrecision -> 880];

How the created Experimental`NumericalFunction should be used:

f[{1/10, 1/5}]

0.00016666666666666666666666666666666666666666666666666666666666666666\
6666666666666666666666666666666666666666666666666666666666666666666666\
6666666666666666666666666666666666666666666666666666666666666666666666\
6666666666666666666666666666666666666666666666666666666666666666666666\
6666666666666666666666666666666666666666666666666666666666666666666666\
6666666666666666666666666666666666666666666666666666666666666666666666\
6666666666666666666666666666666666666666666666666666666666666666666666\
6666666666666666666666666666666666666666664685468521373229053036940442\
3962701216291153088605786704858607575162427563183707930229553541195321\
2101361759645135515775432980342040381738861309064752553797366691813555\
0017169360630698662587053257997111846801265386183928858616194265489581\
2949120884649269257779567658163657155358878262590630724222855526137604\
688700110160340115203634767418946813394731969

Here is a comparison of performance (updated):

<< GeneralUtilities`
ff[a_, b_] := Nest[# + a - b Sin[2 \[Pi] #] &, 0, n]/n;
o1 = ff[1/100, 1/500]; // AccurateTiming
o2 = ff[0.01`880, 1/500]; // AccurateTiming
o3 = f[{1/100, 1/500}]; // AccurateTiming

3.603673
0.0668274
0.068519

o2 == o3

True

As one can see from the timings, unfortunately in this concrete case Experimental`NumericalFunction does not increase performance as compared to the inexact case (o2) based on pure Nest.

It means that the only way to increase performance of computing the ArrayPlot is to use parallelization. According to the documentation, "Outer products automatically parallelize" by Parallelize so the code in the question can be modified in straightforward way:

dat = Parallelize[Outer[f[500, #2, #] &, Range[0, 1.`880, 1/500], Range[0, 1.`880, 1/500]]]; // AbsoluteTiming
ArrayPlot[dat, 
   ColorFunction -> (Blend[{Black, Blue, Green, Yellow, Red}, #] &), 
   ColorFunctionScaling -> False, DataReversed -> True]

---

But even without parallelization it is possible to plot the problematic region in a reasonable time with precision 100 (MMa 8.0.4):

f[n_, a_, b_] := Nest[# + a - b Sin[2 \[Pi] #] &, 0, n]/n;

n = 500; prec = 100;
dat = Outer[f[n, #2, #] &, N[Range[4/10, 85/100, 1/n], prec], 
    N[Range[3/10, 7/10, 1/n], prec]]; // AbsoluteTiming
ArrayPlot[dat, 
 ColorFunction -> (Blend[{Black, Blue, Green, Yellow, Red}, #] &), 
 ColorFunctionScaling -> False, DataReversed -> True]

{501.2916722, Null}

And here is the result with prec = 300:

{1109.0984368, Null}

Incidentally, here is a higher resolution picture of the image produced by heropup's answer, plotted over both positive and negative values of the two circle map parameters:

Blue represents negative, yellow represents positive.

You can view a far larger 6001x6001 pixel image (about 29MB) here.

Here is the code used to generate the upper right quarter of the image (execution time with a 4-core i5-3550 was around 90 minutes):

f = Compile[{{a, _Real}, {k, _Real}}, 
   Mean[Table[
     Nest[# + a - k Sin[2 Pi #] &, j/1600 + RandomReal[{-1, 1}/6400], 
       50]/50, {j, 1600}]], CompilationTarget -> "C", 
   RuntimeAttributes -> {Listable}];
{tim, dat} = 
  AbsoluteTiming[
   Parallelize[
    Outer[f[#2, #] &, N@Range[0, 3, 1/1000], N@Range[0, 3, 1/1000]]]];
Export["ArnoldTongue.h5", dat]

Instead of averaging over random seed values between 0 and 1, this averages uniformly over the interval $[0,1]$ with a slight random noise added in; this accelerates the convergence to a smooth final image as compared to using random seeds, and thus reduces compute time.

The other 3/4 of the image can then be recovered by symmetry and antisymmetry. The resulting data was colored and rendered in Julia.