Avoid possible zero division in pj_qsfn().

Fixes https://oss-fuzz.com/v2/testcase-detail/5224053675655168?noredirect=1

Credit to OSS-Fuzz

1 files changed, 15 insertions, 9 deletions

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);
 }