(Somewhat elderly topic, I realize, but an interesting and useful one.)
If the earth were a sphere, you could get an analytically exact answer by splitting the region into triangles and computing the area of each triangle. You'd do that by computing the angle at each vertex, summing it up, and seeing by how much it exceeds 180 degrees = pi radians. If the "spherical excess" is theta degrees, the area is theta * r^2, where r=earth radius in your desired units.
This has two flaws. If the area is small, it'll be hard to avoid roundoff error (theta+pi, the total of the three vertex angles, will be very close to pi; subtract to get theta, and you're subtracting two nearly-equal quantities.) Also, the earth is not a sphere.
Best I can see for small areas would be to compute the sides of the triangles using, say, the method of Vicenty. (This allows one to compute quite precise distances and angles on an ellipsoid.) Then compute the area of a planar triangle with those sides. The error will be vanishingly low for small triangles.
For larger triangles... well, the best I can come up with is a recursive method that splits them into lots of small triangles :
area of a triangle :
If all three sides are below some threshhold in length, compute the area of a planar triangle with those sides, and
Otherwise, find the midpoint of the longest side. (The Vicenty algorithm makes that possible.) Say that point M lies on AB, opposite vertex C. Create two triangles, AMC and CMB, each with nearly half the area of the starting triangle ABC. Return the sum of the areas of AMC and CMB.