Speed-up inverse merc (spherical and ellipsoidal formulations) for co…

…ordinates with same northing

Same rationale as in OSGeo#2043
When converting by batch of coordinates of 2000 consecutive same northings,
the speed-up over master is 7 times for ellipsoidal formulation, and
about 40% for the spherical formulation (and thus webmerc)
On top of OSGeo#2052, the speed-up for ellipsoidal would also be about 40%

src/projections/merc.cpp

@@ -10,6 +10,13 @@
 PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts=";
 PROJ_HEAD(webmerc, "Web Mercator / Pseudo Mercator") "\n\tCyl, Ell\n\t";
 
+namespace { // anonymous namespace
+struct merc_data {
+    double  y;
+    double  phi;
+};
+}
+
 #define EPS10 1.e-10
 static double logtanpfpim1(double x) {       /* log(tan(x/2 + M_FORTPI)) */
     if (fabs(x) <= DBL_EPSILON) {
@@ -36,7 +43,7 @@ static PJ_XY merc_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward
     if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return xy;
-}
+    }
     xy.x = P->k0 * lp.lam;
     xy.y = P->k0 * logtanpfpim1(lp.phi);
     return xy;
@@ -45,23 +52,47 @@ static PJ_XY merc_s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward
 
 static PJ_LP merc_e_inverse (PJ_XY xy, PJ *P) {          /* Ellipsoidal, inverse */
     PJ_LP lp = {0.0,0.0};
-    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;
-}
+    merc_data *Q = static_cast<merc_data*>(P->opaque);
+    if( xy.y == Q->y ) {
+        lp.phi = Q->phi;
+    } else {
+        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;
+        }
+        Q->y = xy.y;
+        Q->phi = lp.phi;
+    }
     lp.lam = xy.x / P->k0;
     return lp;
 }
 
 
 static PJ_LP merc_s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
     PJ_LP lp = {0.0,0.0};
-    lp.phi = atan(sinh(xy.y / P->k0));
+    merc_data *Q = static_cast<merc_data*>(P->opaque);
+    if( xy.y == Q->y ) {
+        lp.phi = Q->phi;
+    } else {
+        lp.phi = atan(sinh(xy.y / P->k0));
+        Q->y = xy.y;
+        Q->phi = lp.phi;
+    }
     lp.lam = xy.x / P->k0;
     return lp;
 }
 
 
+static PJ* setup_private_member(PJ* P)
+{
+    merc_data *Q = static_cast<merc_data*>(pj_calloc (1, sizeof (merc_data)));
+    if (nullptr==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = Q;
+    return P;
+}
+
+
 PJ *PROJECTION(merc) {
     double phits=0.0;
     int is_phits;
@@ -72,6 +103,10 @@ PJ *PROJECTION(merc) {
             return pj_default_destructor(P, PJD_ERR_LAT_TS_LARGER_THAN_90);
     }
 
+    P = setup_private_member(P);
+    if( !P )
+        return nullptr;
+
     if (P->es != 0.0) { /* ellipsoid */
         if (is_phits)
             P->k0 = pj_msfn(sin(phits), cos(phits), P->es);
@@ -91,6 +126,10 @@ PJ *PROJECTION(merc) {
 
 PJ *PROJECTION(webmerc) {
 
+    P = setup_private_member(P);
+    if( !P )
+        return nullptr;
+
     /* Overriding k_0 with fixed parameter */
     P->k0 = 1.0;
 

test/gie/builtins.gie

@@ -3233,10 +3233,12 @@ expect  -222638.981586547 -110579.965218249
 direction inverse
 accept  200 100
 expect  0.001796631 0.000904369
-accept  200 -100
-expect  0.001796631 -0.000904369
+# Make sure this test immediately follows the previous one so they use the same
+# northing, to test the related optimization
 accept  -200 100
 expect  -0.001796631 0.000904369
+accept  200 -100
+expect  0.001796631 -0.000904369
 accept  -200 -100
 expect  -0.001796631 -0.000904369
 
@@ -3256,10 +3258,12 @@ expect  -223402.144255274 -111706.743574944
 direction inverse
 accept  200 100
 expect  0.001790493 0.000895247
-accept  200 -100
-expect  0.001790493 -0.000895247
+# Make sure this test immediately follows the previous one so they use the same
+# northing, to test the related optimization
 accept  -200 100
 expect  -0.001790493 0.000895247
+accept  200 -100
+expect  0.001790493 -0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247