Merge pull request #519 from kbevers/avoid-zero-division

Handle zero division errors reported by OSS-Fuzz

9 files changed, 204 insertions, 165 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index d23734a2..524657b2 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -30,8 +30,8 @@
 #define PJ_LIB__
 #include <projects.h>
 
-# define EPS10	1.e-10
-# define TOL7	1.e-7
+# define EPS10  1.e-10
+# define TOL7   1.e-7
 
 PROJ_HEAD(aea, "Albers Equal Area") "\n\tConic Sph&Ell\n\tlat_1= lat_2=";
 PROJ_HEAD(leac, "Lambert Equal Area Conic") "\n\tConic, Sph&Ell\n\tlat_1= south";
@@ -42,39 +42,39 @@ PROJ_HEAD(leac, "Lambert Equal Area Conic") "\n\tConic, Sph&Ell\n\tlat_1= south"
 # define EPSILON 1.0e-7
 # define TOL 1.0e-10
 static double phi1_(double qs, double Te, double Tone_es) {
-	int i;
-	double Phi, sinpi, cospi, con, com, dphi;
-
-	Phi = asin (.5 * qs);
-	if (Te < EPSILON)
-		return( Phi );
-	i = N_ITER;
-	do {
-		sinpi = sin (Phi);
-		cospi = cos (Phi);
-		con = Te * sinpi;
-		com = 1. - con * con;
-		dphi = .5 * com * com / cospi * (qs / Tone_es -
-		   sinpi / com + .5 / Te * log ((1. - con) /
-		   (1. + con)));
-		Phi += dphi;
-	} while (fabs(dphi) > TOL && --i);
-	return( i ? Phi : HUGE_VAL );
+    int i;
+    double Phi, sinpi, cospi, con, com, dphi;
+
+    Phi = asin (.5 * qs);
+    if (Te < EPSILON)
+        return( Phi );
+    i = N_ITER;
+    do {
+        sinpi = sin (Phi);
+        cospi = cos (Phi);
+        con = Te * sinpi;
+        com = 1. - con * con;
+        dphi = .5 * com * com / cospi * (qs / Tone_es -
+           sinpi / com + .5 / Te * log ((1. - con) /
+           (1. + con)));
+        Phi += dphi;
+    } while (fabs(dphi) > TOL && --i);
+    return( i ? Phi : HUGE_VAL );
 }
 
 
 struct pj_opaque {
-	double	ec;
-	double	n;
-	double	c;
-	double	dd;
-	double	n2;
-	double	rho0;
-	double	rho;
-	double	phi1;
-	double	phi2;
-	double	*en;
-	int		ellips;
+    double  ec;
+    double  n;
+    double  c;
+    double  dd;
+    double  n2;
+    double  rho0;
+    double  rho;
+    double  phi1;
+    double  phi2;
+    double  *en;
+    int     ellips;
 };
 
 
@@ -82,42 +82,42 @@ struct pj_opaque {
 static XY e_forward (LP lp, PJ *P) {   /* Ellipsoid/spheroid, forward */
     XY xy = {0.0,0.0};
     struct pj_opaque *Q = P->opaque;
-	if ((Q->rho = Q->c - (Q->ellips ? Q->n * pj_qsfn(sin(lp.phi),
-		P->e, P->one_es) : Q->n2 * sin(lp.phi))) < 0.) F_ERROR
-	Q->rho = Q->dd * sqrt(Q->rho);
-	xy.x = Q->rho * sin( lp.lam *= Q->n );
-	xy.y = Q->rho0 - Q->rho * cos(lp.lam);
-	return xy;
+    if ((Q->rho = Q->c - (Q->ellips ? Q->n * pj_qsfn(sin(lp.phi),
+        P->e, P->one_es) : Q->n2 * sin(lp.phi))) < 0.) F_ERROR
+    Q->rho = Q->dd * sqrt(Q->rho);
+    xy.x = Q->rho * sin( lp.lam *= Q->n );
+    xy.y = Q->rho0 - Q->rho * cos(lp.lam);
+    return xy;
 }
 
 
 static LP e_inverse (XY xy, PJ *P) {   /* Ellipsoid/spheroid, inverse */
     LP lp = {0.0,0.0};
     struct pj_opaque *Q = P->opaque;
-	if( (Q->rho = hypot(xy.x, xy.y = Q->rho0 - xy.y)) != 0.0 ) {
-		if (Q->n < 0.) {
-			Q->rho = -Q->rho;
-			xy.x = -xy.x;
-			xy.y = -xy.y;
-		}
-		lp.phi =  Q->rho / Q->dd;
-		if (Q->ellips) {
-			lp.phi = (Q->c - lp.phi * lp.phi) / Q->n;
-			if (fabs(Q->ec - fabs(lp.phi)) > TOL7) {
-				if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL)
-					I_ERROR
-			} else
-				lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
-		} else if (fabs(lp.phi = (Q->c - lp.phi * lp.phi) / Q->n2) <= 1.)
-			lp.phi = asin(lp.phi);
-		else
-			lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
-		lp.lam = atan2(xy.x, xy.y) / Q->n;
-	} else {
-		lp.lam = 0.;
-		lp.phi = Q->n > 0. ? M_HALFPI : - M_HALFPI;
-	}
-	return lp;
+    if( (Q->rho = hypot(xy.x, xy.y = Q->rho0 - xy.y)) != 0.0 ) {
+        if (Q->n < 0.) {
+            Q->rho = -Q->rho;
+            xy.x = -xy.x;
+            xy.y = -xy.y;
+        }
+        lp.phi =  Q->rho / Q->dd;
+        if (Q->ellips) {
+            lp.phi = (Q->c - lp.phi * lp.phi) / Q->n;
+            if (fabs(Q->ec - fabs(lp.phi)) > TOL7) {
+                if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL)
+                    I_ERROR
+            } else
+                lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
+        } else if (fabs(lp.phi = (Q->c - lp.phi * lp.phi) / Q->n2) <= 1.)
+            lp.phi = asin(lp.phi);
+        else
+            lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
+        lp.lam = atan2(xy.x, xy.y) / Q->n;
+    } else {
+        lp.lam = 0.;
+        lp.phi = Q->n > 0. ? M_HALFPI : - M_HALFPI;
+    }
+    return lp;
 }
 
 
@@ -141,47 +141,49 @@ static void freeup (PJ *P) {
 
 
 static PJ *setup(PJ *P) {
-	double cosphi, sinphi;
-	int secant;
+    double cosphi, sinphi;
+    int secant;
     struct pj_opaque *Q = P->opaque;
 
     P->inv = e_inverse;
     P->fwd = e_forward;
 
-	if (fabs(Q->phi1 + Q->phi2) < EPS10) E_ERROR(-21);
-	Q->n = sinphi = sin(Q->phi1);
-	cosphi = cos(Q->phi1);
-	secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
-	if( (Q->ellips = (P->es > 0.))) {
-		double ml1, m1;
-
-		if (!(Q->en = pj_enfn(P->es))) E_ERROR_0;
-		m1 = pj_msfn(sinphi, cosphi, P->es);
-		ml1 = pj_qsfn(sinphi, P->e, P->one_es);
-		if (secant) { /* secant cone */
-			double ml2, m2;
-
-			sinphi = sin(Q->phi2);
-			cosphi = cos(Q->phi2);
-			m2 = pj_msfn(sinphi, cosphi, P->es);
-			ml2 = pj_qsfn(sinphi, P->e, P->one_es);
-			Q->n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
-		}
-		Q->ec = 1. - .5 * P->one_es * log((1. - P->e) /
-			(1. + P->e)) / P->e;
-		Q->c = m1 * m1 + Q->n * ml1;
-		Q->dd = 1. / Q->n;
-		Q->rho0 = Q->dd * sqrt(Q->c - Q->n * pj_qsfn(sin(P->phi0),
-			P->e, P->one_es));
-	} else {
-		if (secant) Q->n = .5 * (Q->n + sin(Q->phi2));
-		Q->n2 = Q->n + Q->n;
-		Q->c = cosphi * cosphi + Q->n2 * sinphi;
-		Q->dd = 1. / Q->n;
-		Q->rho0 = Q->dd * sqrt(Q->c - Q->n2 * sin(P->phi0));
-	}
-
-	return P;
+    if (fabs(Q->phi1 + Q->phi2) < EPS10) E_ERROR(-21);
+    Q->n = sinphi = sin(Q->phi1);
+    cosphi = cos(Q->phi1);
+    secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
+    if( (Q->ellips = (P->es > 0.))) {
+        double ml1, m1;
+
+        if (!(Q->en = pj_enfn(P->es))) E_ERROR_0;
+        m1 = pj_msfn(sinphi, cosphi, P->es);
+        ml1 = pj_qsfn(sinphi, P->e, P->one_es);
+        if (secant) { /* secant cone */
+            double ml2, m2;
+
+            sinphi = sin(Q->phi2);
+            cosphi = cos(Q->phi2);
+            m2 = pj_msfn(sinphi, cosphi, P->es);
+            ml2 = pj_qsfn(sinphi, P->e, P->one_es);
+            if (ml2 == ml1)
+                return NULL;
+            Q->n = (m1 * m1 - m2 * m2) / (ml2 - ml1);
+        }
+        Q->ec = 1. - .5 * P->one_es * log((1. - P->e) /
+            (1. + P->e)) / P->e;
+        Q->c = m1 * m1 + Q->n * ml1;
+        Q->dd = 1. / Q->n;
+        Q->rho0 = Q->dd * sqrt(Q->c - Q->n * pj_qsfn(sin(P->phi0),
+            P->e, P->one_es));
+    } else {
+        if (secant) Q->n = .5 * (Q->n + sin(Q->phi2));
+        Q->n2 = Q->n + Q->n;
+        Q->c = cosphi * cosphi + Q->n2 * sinphi;
+        Q->dd = 1. / Q->n;
+        Q->rho0 = Q->dd * sqrt(Q->c - Q->n2 * sin(P->phi0));
+    }
+
+    return P;
 }
 
 
@@ -191,8 +193,8 @@ PJ *PROJECTION(aea) {
         return freeup_new (P);
     P->opaque = Q;
 
-	Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
-	Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
+    Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
+    Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
     return setup(P);
 }
 
@@ -203,8 +205,8 @@ PJ *PROJECTION(leac) {
         return freeup_new (P);
     P->opaque = Q;
 
-	Q->phi2 = pj_param(P->ctx, P->params, "rlat_1").f;
-	Q->phi1 = pj_param(P->ctx, P->params, "bsouth").i ? - M_HALFPI: M_HALFPI;
+    Q->phi2 = pj_param(P->ctx, P->params, "rlat_1").f;
+    Q->phi1 = pj_param(P->ctx, P->params, "bsouth").i ? - M_HALFPI: M_HALFPI;
     return setup(P);
 }
 
diff --git a/src/PJ_eck3.c b/src/PJ_eck3.c
index d70838d2..3fe5c49f 100644
--- a/src/PJ_eck3.c
+++ b/src/PJ_eck3.c
@@ -24,9 +24,14 @@ static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
 static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
     LP lp = {0.0,0.0};
     struct pj_opaque *Q = P->opaque;
+    double denominator;
 
     lp.phi = xy.y / Q->C_y;
-    lp.lam = xy.x / (Q->C_x * (Q->A + asqrt(1. - Q->B * lp.phi * lp.phi)));
+    denominator = (Q->C_x * (Q->A + asqrt(1. - Q->B * lp.phi * lp.phi)));
+    if ( denominator == 0.0)
+        lp.lam = HUGE_VAL;
+    else
+        lp.lam = xy.x / denominator;
     return lp;
 }
 
diff --git a/src/PJ_stere.c b/src/PJ_stere.c
index c7fb02e3..f338d16d 100644
--- a/src/PJ_stere.c
+++ b/src/PJ_stere.c
@@ -34,7 +34,7 @@ static double ssfn_ (double phit, double sinphi, double eccen) {
 static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
     XY xy = {0.0,0.0};
     struct pj_opaque *Q = P->opaque;
-    double coslam, sinlam, sinX = 0.0, cosX = 0.0, X, A, sinphi;
+    double coslam, sinlam, sinX = 0.0, cosX = 0.0, X, A = 0.0, sinphi;
 
     coslam = cos (lp.lam);
     sinlam = sin (lp.lam);
@@ -52,8 +52,13 @@ static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
         goto xmul; /* but why not just  xy.x = A * cosX; break;  ? */
 
     case EQUIT:
-        A = Q->akm1 / (1. + cosX * coslam);
-        xy.y = A * sinX;
+        /* avoid zero division */
+        if (1. + cosX * coslam == 0.0) {
+            xy.y = HUGE_VAL;
+        } else {
+            A = Q->akm1 / (1. + cosX * coslam);
+            xy.y = A * sinX;
+        }
 xmul:
         xy.x = A * cosX;
         break;
diff --git a/src/geocent.c b/src/geocent.c
index 66e2c314..e340e73d 100644
--- a/src/geocent.c
+++ b/src/geocent.c
@@ -414,6 +414,11 @@ void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi,
 /*  ellipsoidal (geodetic) height */
         *Height = P*CPHI0+Z*SPHI0-RN*(1.0-gi->Geocent_e2*SPHI0*SPHI0);
 
+        /* avoid zero division */
+        if (RN+*Height==0.0) {
+            *Latitude = 0.0;
+            return;
+        }
         RK = gi->Geocent_e2*RN/(RN+*Height);
         RX = 1.0/sqrt(1.0-RK*(2.0-RK)*ST*ST);
         CPHI = ST*(1.0-RK)*RX;
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index 865726fa..72892ac2 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -74,6 +74,10 @@ int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
             *a = sqrt(*a * b);
             *es = 0.;
         } else if (pj_param(ctx,pl, "bR_h").i) { /* sphere--harmonic mean */
+            if ( (*a + b) == 0.0) {
+                pj_ctx_set_errno(ctx, -20);
+                goto bomb;
+            }
             *a = 2. * *a * b / (*a + b);
             *es = 0.;
         } else if ((i = pj_param(ctx,pl, "tR_lat_a").i) || /* sphere--arith. */
diff --git a/src/pj_gauss.c b/src/pj_gauss.c
index 67a1ab07..70b057f3 100644
--- a/src/pj_gauss.c
+++ b/src/pj_gauss.c
@@ -29,65 +29,68 @@
 #define MAX_ITER 20
 
 struct GAUSS {
-	double C;
-	double K;
-	double e;
-	double ratexp;
+    double C;
+    double K;
+    double e;
+    double ratexp;
 };
 #define EN ((struct GAUSS *)en)
 #define DEL_TOL 1e-14
-	static double
-srat(double esinp, double exp) {
-	return(pow((1.-esinp)/(1.+esinp), exp));
+
+static double srat(double esinp, double exp) {
+    return(pow((1.-esinp)/(1.+esinp), exp));
 }
 
-	void *
-pj_gauss_ini(double e, double phi0, double *chi, double *rc) {
-	double sphi, cphi, es;
-	struct GAUSS *en;
+void *pj_gauss_ini(double e, double phi0, double *chi, double *rc) {
+    double sphi, cphi, es;
+    struct GAUSS *en;
 
-	if ((en = (struct GAUSS *)malloc(sizeof(struct GAUSS))) == NULL)
-		return (NULL);
-	es = e * e;
-	EN->e = e;
-	sphi = sin(phi0);
-	cphi = cos(phi0);  cphi *= cphi;
-	*rc = sqrt(1. - es) / (1. - es * sphi * sphi);
-	EN->C = sqrt(1. + es * cphi * cphi / (1. - es));
-	*chi = asin(sphi / EN->C);
-	EN->ratexp = 0.5 * EN->C * e;
-	EN->K = tan(.5 * *chi + M_FORTPI) / (
-		pow(tan(.5 * phi0 + M_FORTPI), EN->C) *
-		srat(EN->e * sphi, EN->ratexp)  );
-	return ((void *)en);
+    if ((en = (struct GAUSS *)malloc(sizeof(struct GAUSS))) == NULL)
+        return (NULL);
+    es = e * e;
+    EN->e = e;
+    sphi = sin(phi0);
+    cphi = cos(phi0);  cphi *= cphi;
+    *rc = sqrt(1. - es) / (1. - es * sphi * sphi);
+    EN->C = sqrt(1. + es * cphi * cphi / (1. - es));
+    if (en->C == 0.0) {
+        free(en);
+        return NULL;
+    }
+    *chi = asin(sphi / EN->C);
+    EN->ratexp = 0.5 * EN->C * e;
+    EN->K = tan(.5 * *chi + M_FORTPI) / (
+        pow(tan(.5 * phi0 + M_FORTPI), EN->C) *
+        srat(EN->e * sphi, EN->ratexp)  );
+    return ((void *)en);
 }
-	LP
-pj_gauss(projCtx ctx, LP elp, const void *en) {
-	LP slp;
-	(void) ctx;
 
-	slp.phi = 2. * atan( EN->K *
-		pow(tan(.5 * elp.phi + M_FORTPI), EN->C) *
-		srat(EN->e * sin(elp.phi), EN->ratexp) ) - M_HALFPI;
-	slp.lam = EN->C * (elp.lam);
-	return(slp);
+LP pj_gauss(projCtx ctx, LP elp, const void *en) {
+    LP slp;
+    (void) ctx;
+
+    slp.phi = 2. * atan( EN->K *
+        pow(tan(.5 * elp.phi + M_FORTPI), EN->C) *
+        srat(EN->e * sin(elp.phi), EN->ratexp) ) - M_HALFPI;
+    slp.lam = EN->C * (elp.lam);
+    return(slp);
 }
-	LP
-pj_inv_gauss(projCtx ctx, LP slp, const void *en) {
-	LP elp;
-	double num;
-	int i;
 
-	elp.lam = slp.lam / EN->C;
-	num = pow(tan(.5 * slp.phi + M_FORTPI)/EN->K, 1./EN->C);
-	for (i = MAX_ITER; i; --i) {
-		elp.phi = 2. * atan(num * srat(EN->e * sin(slp.phi), -.5 * EN->e))
-			- M_HALFPI;
-		if (fabs(elp.phi - slp.phi) < DEL_TOL) break;
-			slp.phi = elp.phi;
-	}	
-	/* convergence failed */
-	if (!i)
-		pj_ctx_set_errno( ctx, -17 );
-	return (elp);
+LP pj_inv_gauss(projCtx ctx, LP slp, const void *en) {
+    LP elp;
+    double num;
+    int i;
+
+    elp.lam = slp.lam / EN->C;
+    num = pow(tan(.5 * slp.phi + M_FORTPI)/EN->K, 1./EN->C);
+    for (i = MAX_ITER; i; --i) {
+        elp.phi = 2. * atan(num * srat(EN->e * sin(slp.phi), -.5 * EN->e))
+            - M_HALFPI;
+        if (fabs(elp.phi - slp.phi) < DEL_TOL) break;
+            slp.phi = elp.phi;
+    }
+    /* convergence failed */
+    if (!i)
+        pj_ctx_set_errno( ctx, -17 );
+    return (elp);
 }
diff --git a/src/pj_init.c b/src/pj_init.c
index 6568033f..113cc9f4 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -626,6 +626,10 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         PIN->to_meter = pj_strtod(s, &s);
         if (*s == '/') /* ratio number */
             PIN->to_meter /= pj_strtod(++s, 0);
+        if (PIN->to_meter <= 0.0) {
+            pj_ctx_set_errno( ctx, -51);
+            goto bum_call;
+        }
         PIN->fr_meter = 1. / PIN->to_meter;
     } else
         PIN->to_meter = PIN->fr_meter = 1.;
@@ -641,6 +645,10 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         PIN->vto_meter = pj_strtod(s, &s);
         if (*s == '/') /* ratio number */
             PIN->vto_meter /= pj_strtod(++s, 0);
+        if (PIN->vto_meter <= 0.0) {
+            pj_ctx_set_errno( ctx, -51);
+            goto bum_call;
+        }
         PIN->vfr_meter = 1. / PIN->vto_meter;
     } else {
         PIN->vto_meter = PIN->to_meter;
diff --git a/src/pj_qsfn.c b/src/pj_qsfn.c
index ccb12308..c3610140 100644
--- a/src/pj_qsfn.c
+++ b/src/pj_qsfn.c
@@ -3,14 +3,20 @@
 #include <projects.h>
 
 # define EPSILON 1.0e-7
-	double
-pj_qsfn(double sinphi, double e, double one_es) {
-	double con;
 
-	if (e >= EPSILON) {
-		con = e * sinphi;
-		return (one_es * (sinphi / (1. - con * con) -
-		   (.5 / e) * log ((1. - con) / (1. + con))));
-	} else
-		return (sinphi + sinphi);
+double pj_qsfn(double sinphi, double e, double one_es) {
+    double con, div1, div2;
+
+    if (e >= EPSILON) {
+        con = e * sinphi;
+        div1 = 1.0 - con * con;
+        div2 = 1.0 + con;
+
+        /* avoid zero division, fail gracefully */
+        if (div1 == 0.0 || div2 == 0.0)
+            return HUGE_VAL;
+
+        return (one_es * (sinphi / div1 - (.5 / e) * log ((1. - con) / div2 )));
+    } else
+        return (sinphi + sinphi);
 }
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index 8a2a9c4b..36b7de8a 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -56,6 +56,7 @@ pj_err_list[] = {
     "point not within available datum shift grids",                    /* -48 */
     "invalid sweep axis, choose x or y",                               /* -49 */
     "malformed pipeline",                                              /* -50 */
+    "unit conversion factor must be > 0",                              /* -51 */
 };
 
 char *pj_strerrno(int err) {