Improve numerical precision of inverse spherical mercator

Complements f2b3604

src/PJ_merc.c

@@ -20,6 +20,14 @@ static double logtanpfpim1(double x) {       /* log(tan(x/2 + M_FORTPI)) */
     return log(tan(M_FORTPI + .5 * x));
 }
 
+static double hpip2atanexp(double x) {       /* M_HALFPI - 2. * atan(exp(x)) */
+    /* atan(exp(x))  can be approximated by (exp(x)-1)/2 + pi/4 */
+    if (fabs(x) <= DBL_EPSILON)
+        return -expm1(x);
+
+    return M_HALFPI - 2. * atan(exp(x));
+}
+
 static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
     XY xy = {0.0,0.0};
     if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
@@ -49,15 +57,15 @@ static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
     if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL) {
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return lp;
-}
+    }
     lp.lam = xy.x / P->k0;
     return lp;
 }
 
 
 static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
     LP lp = {0.0,0.0};
-    lp.phi = M_HALFPI - 2. * atan(exp(-xy.y / P->k0));
+    lp.phi = hpip2atanexp(-xy.y / P->k0);
     lp.lam = xy.x / P->k0;
     return lp;
 }

src/pj_math.c

@@ -69,6 +69,15 @@ double pj_log1p(double x) {
     return z == 0 ? x : x * log(y) / z;
 }
 
+/* Compute exp(x)-1 accurately */
+double pj_expm1(double x) {
+    /* For x << 1, Taylor expansion of exp(x) gives exp(x) ~ 1 + x */
+    if (fabs(x) <= DBL_EPSILON) {
+        return x;
+    }
+    return exp(x) - 1.0;
+}
+
 /* Compute asinh(x) accurately */
 double pj_asinh(double x) {
     double y = fabs(x); /* Enforce odd parity */

src/proj_math.h

@@ -47,13 +47,15 @@ extern "C" {
 
 double  pj_hypot(double x, double y);
 double  pj_log1p(double x);
+double  pj_expm1(double x);
 double  pj_asinh(double x);
 double  pj_round(double x);
 long    pj_lround(double x);
 int     pj_isnan(double x);
 
 #define hypot   pj_hypot
 #define log1p   pj_log1p
+#define expm1   pj_expm1
 #define asinh   pj_asinh
 #define round   pj_round
 #define lround  pj_lround