After some looking around at Wikipedia and the same question/answer at StackOverflow, I figured I would take a stab at it, and try to fill in the gaps.
First off, Not sure where you got the output, but it appears to be wrong. I plotted the points in ArcMap, buffered them to the distances specified, ran intersect to on the buffers, and then captured the vertex of intersection to get the solutions. Your proposed output is the point in green. I calculated the value in the callout box, which is about 3 meters of what ArcMap gave for solution derived from the intersect.
The math on the wikipedia page isn't too bad, just need to covert your geodetic coordinates to the cartesian ECEF, which can be found here. the a/x +h terms can be replaced by the authalic sphere radius, if you aren't using an ellipsoid.
Probably easiest just give you some well(?) documented code, so here it is in python
import math
import numpy
#assuming elevation = 0
earthR = 6371
LatA = 37.418436
LonA = -121.963477
DistA = 0.265710701754
LatB = 37.417243
LonB = -121.961889
DistB = 0.234592423446
LatC = 37.418692
LonC = -121.960194
DistC = 0.0548954278262
#using authalic sphere
#if using an ellipsoid this step is slightly different
#Convert geodetic Lat/Long to ECEF xyz
# 1. Convert Lat/Long to radians
# 2. Convert Lat/Long(radians) to ECEF
xA = earthR *(math.cos(math.radians(LatA)) * math.cos(math.radians(LonA)))
yA = earthR *(math.cos(math.radians(LatA)) * math.sin(math.radians(LonA)))
zA = earthR *(math.sin(math.radians(LatA)))
xB = earthR *(math.cos(math.radians(LatB)) * math.cos(math.radians(LonB)))
yB = earthR *(math.cos(math.radians(LatB)) * math.sin(math.radians(LonB)))
zB = earthR *(math.sin(math.radians(LatB)))
xC = earthR *(math.cos(math.radians(LatC)) * math.cos(math.radians(LonC)))
yC = earthR *(math.cos(math.radians(LatC)) * math.sin(math.radians(LonC)))
zC = earthR *(math.sin(math.radians(LatC)))
P1 = numpy.array([xA, yA, zA])
P2 = numpy.array([xB, yB, zB])
P3 = numpy.array([xC, yC, zC])
#from wikipedia
#transform to get circle 1 at origin
#transform to get circle 2 on x axis
ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1))
i = numpy.dot(ex, P3 - P1)
ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex))
ez = numpy.cross(ex,ey)
d = numpy.linalg.norm(P2 - P1)
j = numpy.dot(ey, P3 - P1)
#from wikipedia
#plug and chug using above values
x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d)
y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i/j)*x)
# only one case shown here
z = numpy.sqrt(pow(DistA,2) - pow(x,2) - pow(y,2))
#triPt is an array with ECEF x,y,z of trilateration point
triPt = P1 + x*ex + y*ey + z*ez
#convert back to lat/long from ECEF
#convert to degrees
lat = math.degrees(math.asin(triPt[2] / earthR))
lon = math.degrees(math.atan2(triPt[1],triPt[0]))
print lat, lon
Yes, you will get these kinds of errors with a global Mercator projection: it is accurate at the equator and the distortion increases exponentially with latitude away from the equator. The distance distortion is exactly 2 (100%) at 60 degrees latitude. At your test latitudes (64.14 degrees) I compute a distortion of 2.294, exactly agreeing with the ratio 904/394 = 2.294. (Earlier I computed 2.301 but that was based on a sphere, not the WGS84 ellipsoid. The difference (of 0.3%) gives us a sense of the accuracy you might gain from using an ellipsoid-based projection versus the sphere-based Haversine formula.)
There is no such thing as a global projection that yields highly accurate distances everywhere. That's one reason the UTM zone system is used!
One solution is to use spherical geometry for all your calculations, but you have rejected that (which is reasonable if you're going to be doing complex operations, but the decision might be worth revisiting).
Another solution is to adapt the projection to the points being compared. For example, you could safely use a transverse Mercator (as in the UTM system) with a meridian lying near the center of the region of interest. Moving the meridian is a simple thing to do: just subtract the meridian's longitude from all the longitudes and use a single TM projection centered at the Prime Meridian (with a scale factor of 1, rather than the 0.9996 of the UTM system). For your work this will tend to be more accurate than using UTM itself. It will give correct angles (TM is conformal) and will be remarkably accurate for points separated by only a few tens of kilometers: expect better than six-digit accuracy. In fact, I would be inclined to attribute any small differences between these adapted-TM distances and the Haversine distances to the difference between the ellipsoid (used for the TM projection) and the sphere (used by Haversine), rather than to distortion in the projection.
Best Answer
When you find distance using latitude and longitude, the tool would be calculating Geodetic distance between the two points, so it does not need a projection. However it will need the datum or ellipsoid, say WGS84 etc. Since it is not asking for a datum, it signifies a default system datum which generally is WGS84, the world datum.
Once you apply a projection system, then the co-ordinate system changes to planar co-ordinates and the distance can be calculated by applying Cartesian co-ordinate geometry. The distance if measured in planar co-ordinates will change if you change the projection system. But the distance between two points when using Geodetic distance will change only when you change the datum/ellipsoid, projection will have little difference or may be nil.