You need not do this in double precision: single precision is enough! The trouble is what you're using the precision for, and by trying to represent cosines near 1 rather than near 0 you're throwing away all the precision available to you. Navigators in days of yore—with much less precision than modern IEEE 754 binary32 single precision!—had the same problem, and what did they do?

Let's write this as a formula without all the unit conversions—you have latitudes 𝜑₁ and 𝜑₂ and longitudes 𝜆₁ and 𝜆₂ and you want the great-circle angle 𝜃, which you've decided to approach with the relation (probably derived from converting the spherical law of cosines naively to a longitude/latitude coordinate system):

cos 𝜃 = sin 𝜑₁ sin 𝜑₂ + cos 𝜑₁ cos 𝜑₂ cos(𝜆₂ − 𝜆₁).

The crux of the problem is that you're trying to find an angle whose cosine is very close to 1, and you're doing it by computing that cosine. In your example, sin 𝜑₁ sin 𝜑₂ + cos 𝜑₁ cos 𝜑₂ cos(𝜆₂ − 𝜆₁) comes out to roughly 0.999999937 = 1 − 6.3×10⁻⁸, but the single-precision floating-point numbers in [½, 1] can only discern increments of about 6×10⁻⁸, so you've essentially wasted all your precision on those leading nines. Navigators in days of yore didn't even have eight digits (for example, the 1913 tables of E.R. Hedrick and the 1800 tables of Josef de Mendoza y Ríos provide only around five digits), so how did they do it?

You probably grew up on the three trig functions sine, cosine, and tangent if you're much less than a century old, but there are others like haversine (hav 𝜃 = ½ (1 − cos 𝜃)) and exsecant (exsec 𝜃 = sec 𝜃 − 1 = 1/cos 𝜃 − 1)—invented not to torture students in school, but because they made calculations feasible with limited precision! The trick is not to compute the cosine 1 − 𝛿 for very small 𝛿, but instead compute the versine 𝛿, or the haversine ½𝛿, so you can use the full precision available to you to represent it.

The venerable haversine formula relates the latitudes 𝜑₁ and 𝜑₂ and longitudes 𝜆₁ and 𝜆₂ to the great-circle angle 𝜃 by:

hav 𝜃 = hav(𝜑₂ − 𝜑₁) + cos 𝜑₁ cos 𝜑₂ hav(𝜆₂ − 𝜆₁).

This is better for nearby angles because for small angle differences, haversine preserves most of the precision of the input to return a result near 0, while cosine throws it away to return a result near 1.

Of course, the Android math library might be missing a haversine function (pity!), but you can recover it from the trigonometric identity:

hav 𝜃 = ½ (1 − cos 𝜃) = sin²(𝜃/2).

You can always divide by two exactly in floating-point arithmetic (unless it underflows), and squaring won't amplify the error much, so the formula sin²(𝜃/2) will work reliably, and you can recover the inverse archav 𝑥 = 2 arcsin(√𝑥), using functions that probably are in the Android math library. And this identity—together with coordinate transformation between longitude/latitude and spherical triangle angles—is how you can derive the haversine formula from (e.g.) the spherical law of cosines, as the Wikipedia article shows.

Once you use the haversine formula, even if you compute everything in single precision, you get within about 30cm of the exact result in your example, or <0.014% relative error—which is smaller than the error of about 60cm arising from using a spherical approximation to the oblate spheroid of Earth!