Move inline code for optimize mlfn from tmerc to mlfn.hpp to avoid code duplication

5 files changed, 90 insertions, 94 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index fb4a104f..e5f41546 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -200,7 +200,7 @@ libproj_la_SOURCES = \
 	dmstor.cpp auth.cpp \
 	deriv.cpp ell_set.cpp ellps.cpp errno.cpp \
 	factors.cpp fwd.cpp init.cpp inv.cpp \
-	list.cpp malloc.cpp mlfn.cpp msfn.cpp proj_mdist.cpp \
+	list.cpp malloc.cpp mlfn.cpp mlfn.hpp msfn.cpp proj_mdist.cpp \
 	param.cpp phi2.cpp pr_list.cpp \
 	qsfn.cpp strerrno.cpp \
 	tsfn.cpp units.cpp ctx.cpp log.cpp zpoly1.cpp rtodms.cpp \
diff --git a/src/mlfn.cpp b/src/mlfn.cpp
index a5448e3b..80f9163b 100644
--- a/src/mlfn.cpp
+++ b/src/mlfn.cpp
@@ -2,6 +2,7 @@
 
 #include "proj.h"
 #include "proj_internal.h"
+#include "mlfn.hpp"
 
 /* meridional distance for ellipsoid and inverse
 **	8th degree - accurate to < 1e-5 meters when used in conjunction
@@ -20,46 +21,32 @@
 #define C66 .36458333333333333333
 #define C68 .00569661458333333333
 #define C88 .3076171875
-#define EPS 1e-11
-#define MAX_ITER 10
 #define EN_SIZE 5
 
 double *pj_enfn(double es) {
     double t, *en;
 
-	en = (double *) pj_malloc(EN_SIZE * sizeof (double));
-	if (nullptr==en)
-	    return nullptr;
+    en = (double *) pj_malloc(EN_SIZE * sizeof (double));
+    if (nullptr==en)
+        return nullptr;
 
     en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08)));
     en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08)));
     en[2] = (t = es * es) * (C44 - es * (C46 + es * C48));
     en[3] = (t *= es) * (C66 - es * C68);
-	en[4] = t * es * C88;
+    en[4] = t * es * C88;
 
     return en;
 }
 
-	double
-pj_mlfn(double phi, double sphi, double cphi, double *en) {
-	cphi *= sphi;
-	sphi *= sphi;
-	return(en[0] * phi - cphi * (en[1] + sphi*(en[2]
-		+ sphi*(en[3] + sphi*en[4]))));
+double
+pj_mlfn(double phi, double sphi, double cphi, const double *en) {
+    return inline_pj_mlfn(phi, sphi, cphi, en);
 }
-	double
-pj_inv_mlfn(projCtx ctx, double arg, double es, double *en) {
-	double s, t, phi, k = 1./(1.-es);
-	int i;
 
-	phi = arg;
-	for (i = MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
-		s = sin(phi);
-		t = 1. - es * s * s;
-		phi -= t = (pj_mlfn(phi, s, cos(phi), en) - arg) * (t * sqrt(t)) * k;
-		if (fabs(t) < EPS)
-			return phi;
-	}
-	pj_ctx_set_errno( ctx, PJD_ERR_NON_CONV_INV_MERI_DIST );
-	return phi;
+double
+pj_inv_mlfn(projCtx ctx, double arg, double es, const double *en) {
+    double sinphi_ignored;
+    double cosphi_ignored;
+    return inline_pj_inv_mlfn(ctx, arg, es, en, &sinphi_ignored, &cosphi_ignored);
 }
diff --git a/src/mlfn.hpp b/src/mlfn.hpp
new file mode 100644
index 00000000..26a2959f
--- /dev/null
+++ b/src/mlfn.hpp
@@ -0,0 +1,73 @@
+#ifndef MLFN_HPP
+#define MLFN_HPP
+
+/* meridional distance for ellipsoid and inverse
+**	8th degree - accurate to < 1e-5 meters when used in conjunction
+**		with typical major axis values.
+**	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
+*/
+
+inline static double inline_pj_mlfn(double phi, double sphi, double cphi, const double *en) {
+    cphi *= sphi;
+    sphi *= sphi;
+    return(en[0] * phi - cphi * (en[1] + sphi*(en[2]
+            + sphi*(en[3] + sphi*en[4]))));
+}
+
+inline static double
+inline_pj_inv_mlfn(projCtx ctx, double arg, double es, const double *en,
+                   double* sinphi, double* cosphi) {
+    const double k = 1./(1.-es);
+    constexpr double INV_MLFN_EPS = 1e-11;
+    constexpr int INV_MFN_MAX_ITER = 10;
+    double phi = arg;
+    double s = sin(phi);
+    double c = cos(phi);
+    for (int i = INV_MFN_MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
+        double t = 1. - es * s * s;
+        t = (inline_pj_mlfn(phi, s, c, en) - arg) * (t * sqrt(t)) * k;
+        phi -= t;
+        if (fabs(t) < INV_MLFN_EPS)
+        {
+            // Instead of recomputing sin(phi), cos(phi) from scratch,
+            // use sin(phi-t) and cos(phi-t) approximate formulas with
+            // 1-term approximation of sin(t) and cos(t)
+            *sinphi = s - c * t;
+            *cosphi = c + s * t;
+            return phi;
+        }
+        if (fabs(t) < 1e-3)
+        {
+            // 2-term approximation of sin(t) and cos(t)
+            // Max relative error is 4e-14 on cos(t), and 8e-15 on sin(t)
+            const double t2 = t * t;
+            const double cos_t = 1 - 0.5 * t2;
+            const double sin_t = t * (1 - (1. / 6) * t2);
+            const double s_new = s * cos_t - c * sin_t;
+            c = c * cos_t + s * sin_t;
+            s = s_new;
+        }
+        else if (fabs(t) < 1e-2)
+        {
+            // 3-term approximation of sin(t) and cos(t)
+            // Max relative error is 2e-15 on cos(t), and 2e-16 on sin(t)
+            const double t2 = t * t;
+            const double cos_t = 1 - 0.5 * t2 * (1 - (1. / 12) * t2);
+            const double sin_t = t * (1 - (1. / 6) * t2 * (1 - (1. / 20) * t2));
+            const double s_new = s * cos_t - c * sin_t;
+            c = c * cos_t + s * sin_t;
+            s = s_new;
+        }
+        else
+        {
+            s = sin(phi);
+            c = cos(phi);
+        }
+    }
+    *sinphi = s;
+    *cosphi = c;
+    pj_ctx_set_errno( ctx, PJD_ERR_NON_CONV_INV_MERI_DIST );
+    return phi;
+}
+
+#endif // MLFN_HPP
diff --git a/src/proj_internal.h b/src/proj_internal.h
index a7c1ddc8..c4af479f 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -823,8 +823,8 @@ void     *pj_dealloc_params (projCtx_t *ctx, paralist *start, int errlev);
 
 
 double *pj_enfn(double);
-double  pj_mlfn(double, double, double, double *);
-double  pj_inv_mlfn(projCtx_t *, double, double, double *);
+double  pj_mlfn(double, double, double, const double *);
+double  pj_inv_mlfn(projCtx_t *, double, double, const double *);
 double  pj_qsfn(double, double, double);
 double  pj_tsfn(double, double, double);
 double  pj_msfn(double, double, double);
diff --git a/src/projections/tmerc.cpp b/src/projections/tmerc.cpp
index 75a77ccd..754fe53f 100644
--- a/src/projections/tmerc.cpp
+++ b/src/projections/tmerc.cpp
@@ -18,7 +18,7 @@
 #include "proj.h"
 #include "proj_internal.h"
 #include <math.h>
-
+#include "mlfn.hpp"
 
 PROJ_HEAD(tmerc, "Transverse Mercator") "\n\tCyl, Sph&Ell\n\tapprox";
 PROJ_HEAD(etmerc, "Extended Transverse Mercator") "\n\tCyl, Sph";
@@ -72,13 +72,6 @@ struct pj_opaque_exact {
 /*****************************************************************************/
 
 
-inline static double inline_pj_mlfn(double phi, double sphi, double cphi, double *en) {
-    cphi *= sphi;
-    sphi *= sphi;
-    return(en[0] * phi - cphi * (en[1] + sphi*(en[2]
-            + sphi*(en[3] + sphi*en[4]))));
-}
-
 #ifdef BUILD_FMA_OPTIMIZED_VERSION
 __attribute__((target_clones("fma","default")))
 #endif
@@ -173,63 +166,6 @@ static PJ_XY approx_s_fwd (PJ_LP lp, PJ *P) {
     return xy;
 }
 
-inline static double
-inline_pj_inv_mlfn(projCtx ctx, double arg, double es, double *en,
-                   double* sinphi, double* cosphi) {
-    double phi, k = 1./(1.-es);
-    int i;
-#define EPS 1e-11
-#define MAX_ITER 10
-    phi = arg;
-    double s = sin(phi);
-    double c = cos(phi);
-    for (i = MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
-        double t = 1. - es * s * s;
-        t = (inline_pj_mlfn(phi, s, c, en) - arg) * (t * sqrt(t)) * k;
-        phi -= t;
-        if (fabs(t) < EPS)
-        {
-            // Instead of recomputing sin(phi), cos(phi) from scratch,
-            // use sin(phi-t) and cos(phi-t) approximate formulas with
-            // 1-term approximation of sin(t) and cos(t)
-            *sinphi = s - c * t;
-            *cosphi = c + s * t;
-            return phi;
-        }
-        if (fabs(t) < 1e-3)
-        {
-            // 2-term approximation of sin(t) and cos(t)
-            // Max relative error is 4e-14 on cos(t), and 8e-15 on sin(t)
-            const double t2 = t * t;
-            const double cos_t = 1 - 0.5 * t2;
-            const double sin_t = t * (1 - (1. / 6) * t2);
-            const double s_new = s * cos_t - c * sin_t;
-            c = c * cos_t + s * sin_t;
-            s = s_new;
-        }
-        else if (fabs(t) < 1e-2)
-        {
-            // 3-term approximation of sin(t) and cos(t)
-            // Max relative error is 2e-15 on cos(t), and 2e-16 on sin(t)
-            const double t2 = t * t;
-            const double cos_t = 1 - 0.5 * t2 * (1 - (1. / 12) * t2);
-            const double sin_t = t * (1 - (1. / 6) * t2 * (1 - (1. / 20) * t2));
-            const double s_new = s * cos_t - c * sin_t;
-            c = c * cos_t + s * sin_t;
-            s = s_new;
-        }
-        else
-        {
-            s = sin(phi);
-            c = cos(phi);
-        }
-    }
-    *sinphi = s;
-    *cosphi = c;
-    pj_ctx_set_errno( ctx, PJD_ERR_NON_CONV_INV_MERI_DIST );
-    return phi;
-}
-
 #ifdef BUILD_FMA_OPTIMIZED_VERSION
 __attribute__((target_clones("fma","default")))
 #endif