4 files changed, 31 insertions, 7 deletions
diff --git a/src/projections/lcc.cpp b/src/projections/lcc.cpp
index 5eee0d14..55d28b80 100644
--- a/src/projections/lcc.cpp
+++ b/src/projections/lcc.cpp
@@ -105,6 +105,9 @@ PJ *PROJECTION(lcc) {
 
         m1 = pj_msfn(sinphi, cosphi, P->es);
         ml1 = pj_tsfn(Q->phi1, sinphi, P->e);
+        if( ml1 == 0 ) {
+            return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
+        }
         if (secant) { /* secant cone */
             sinphi = sin(Q->phi2);
             Q->n = log(m1 / pj_msfn(sinphi, cos(Q->phi2), P->es));
@@ -112,7 +115,11 @@ PJ *PROJECTION(lcc) {
                 // Not quite, but es is very close to 1...
                 return pj_default_destructor(P, PJD_ERR_ECCENTRICITY_IS_ONE_OR_GREATER);
             }
-            Q->n /= log(ml1 / pj_tsfn(Q->phi2, sinphi, P->e));
+            const double ml2 = pj_tsfn(Q->phi2, sinphi, P->e);
+            if( ml2 == 0 ) {
+                return pj_default_destructor(P, PJD_ERR_LAT_1_OR_2_ZERO_OR_90);
+            }
+            Q->n /= log(ml1 / ml2);
         }
         Q->c = (Q->rho0 = m1 * pow(ml1, -Q->n) / Q->n);
         Q->rho0 *= (fabs(fabs(P->phi0) - M_HALFPI) < EPS10) ? 0. :
diff --git a/src/projections/stere.cpp b/src/projections/stere.cpp
index 9836f341..fd9f9827 100644
--- a/src/projections/stere.cpp
+++ b/src/projections/stere.cpp
@@ -55,11 +55,18 @@ static PJ_XY e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
     }
 
     switch (Q->mode) {
-    case OBLIQ:
-        A = Q->akm1 / (Q->cosX1 * (1. + Q->sinX1 * sinX +
-           Q->cosX1 * cosX * coslam));
+    case OBLIQ: {
+        const double denom = Q->cosX1 * (1. + Q->sinX1 * sinX +
+           Q->cosX1 * cosX * coslam);
+        if( denom == 0 ) {
+            proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+            return proj_coord_error().xy;
+        }
+        A = Q->akm1 / denom;
         xy.y = A * (Q->cosX1 * sinX - Q->sinX1 * cosX * coslam);
-        goto xmul; /* but why not just  xy.x = A * cosX; break;  ? */
+        xy.x = A * cosX;
+        break;
+    }
 
     case EQUIT:
         /* avoid zero division */
@@ -69,7 +76,6 @@ static PJ_XY e_forward (PJ_LP lp, PJ *P) {          /* Ellipsoidal, forward */
             A = Q->akm1 / (1. + cosX * coslam);
             xy.y = A * sinX;
         }
-xmul:
         xy.x = A * cosX;
         break;
 
diff --git a/src/projections/tmerc.cpp b/src/projections/tmerc.cpp
index c91c5174..bb56f8ae 100644
--- a/src/projections/tmerc.cpp
+++ b/src/projections/tmerc.cpp
@@ -188,6 +188,10 @@ static PJ_LP approx_s_inv (PJ_XY xy, PJ *P) {
     double h, g;
 
     h = exp(xy.x / static_cast<struct pj_opaque_approx*>(P->opaque)->esp);
+    if( h == 0 ) {
+        proj_errno_set(P, PJD_ERR_INVALID_X_OR_Y);
+        return proj_coord_error().lp;
+    }
     g = .5 * (h - 1. / h);
     h = cos (P->phi0 + xy.y / static_cast<struct pj_opaque_approx*>(P->opaque)->esp);
     lp.phi = asin(sqrt((1. - h * h) / (1. + g * g)));
diff --git a/src/projections/vandg.cpp b/src/projections/vandg.cpp
index 89620356..c669f8fa 100644
--- a/src/projections/vandg.cpp
+++ b/src/projections/vandg.cpp
@@ -80,7 +80,14 @@ static PJ_LP s_inverse (PJ_XY xy, PJ *P) {           /* Spheroidal, inverse */
     al = c1 / c3 - THIRD * c2 * c2;
     m = 2. * sqrt(-THIRD * al);
     d = C2_27 * c2 * c2 * c2 + (c0 * c0 - THIRD * c2 * c1) / c3;
-    if (((t = fabs(d = 3. * d / (al * m))) - TOL) <= 1.) {
+    const double al_mul_m = al * m;
+    if( al_mul_m == 0 ) {
+        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+        return proj_coord_error().lp;
+    }
+    d = 3. * d /al_mul_m;
+    t = fabs(d);
+    if ((t - TOL) <= 1.) {
         d = t > 1. ? (d > 0. ? 0. : M_PI) : acos(d);
         lp.phi = M_PI * (m * cos(d * THIRD + PI4_3) - THIRD * c2);
         if (xy.y < 0.) lp.phi = -lp.phi;