It is more correct to solve this problem not by increasing the accuracy of calculations, but by transforming the formula to a numerically stable form. That is, this is a mathematical problem. The answer to it is as follows.

S = R*arcctg ((b-a)/sqrt(4*a*b+c^2)),

where (x is the latitude, and y is the longitude):

a = cos(x1)*cos(x2)*sin^2((y1-y2)/2),

b = cos(x1-x2)-a,

c = sin(x1-x2).

The most convenient way to calculate arccotangent in C++ is to use the atan2() function, passing the denominator and numerator of the expression under the arcctg symbol as arguments.