Add upper limit to number of iterations for all (documented as such) computations using Newton-Raphson. Credit to OSS Fuzz

5 files changed, 29 insertions, 6 deletions

diff --git a/src/PJ_aitoff.c b/src/PJ_aitoff.c
index faab9aba..77d3f1c6 100644
--- a/src/PJ_aitoff.c
+++ b/src/PJ_aitoff.c
@@ -141,7 +141,11 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
 	/* if too far from given values of x,y, repeat with better approximation of phi,lam */
 	} while (((fabs(xy.x-x) > EPSILON) || (fabs(xy.y-y) > EPSILON)) && (round++ < MAXROUND));
 
-	if (iter == MAXITER && round == MAXROUND) fprintf(stderr, "Warning: Accuracy of 1e-12 not reached. Last increments: dlat=%e and dlon=%e\n", dp, dl);
+	if (iter == MAXITER && round == MAXROUND)
+        {
+            pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );
+            /* fprintf(stderr, "Warning: Accuracy of 1e-12 not reached. Last increments: dlat=%e and dlon=%e\n", dp, dl); */
+        }
 
 	return lp;
 }
diff --git a/src/PJ_natearth.c b/src/PJ_natearth.c
index 22bf340d..8d0dae08 100644
--- a/src/PJ_natearth.c
+++ b/src/PJ_natearth.c
@@ -34,6 +34,8 @@ PROJ_HEAD(natearth, "Natural Earth") "\n\tPCyl., Sph.";
 #define C4 (11 * B4)
 #define EPS 1e-11
 #define MAX_Y (0.8707 * 0.52 * M_PI)
+/* Not sure at all of the appropriate number for MAX_ITER... */
+#define MAX_ITER 100
 
 
 static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
@@ -52,6 +54,7 @@ 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};
     double yc, tol, y2, y4, f, fder;
+    int i;
     (void) P;
 
     /* make sure y is inside valid range */
@@ -63,7 +66,7 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
 
     /* latitude */
     yc = xy.y;
-        for (;;) { /* Newton-Raphson */
+    for (i = MAX_ITER; i ; --i) { /* Newton-Raphson */
         y2 = yc * yc;
         y4 = y2 * y2;
         f = (yc * (B0 + y2 * (B1 + y4 * (B2 + B3 * y2 + B4 * y4)))) - xy.y;
@@ -73,6 +76,8 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
             break;
         }
     }
+    if( i == 0 )
+        pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );
     lp.phi = yc;
 
     /* longitude */
diff --git a/src/PJ_natearth2.c b/src/PJ_natearth2.c
index 9cd93fdf..aab9f9ac 100644
--- a/src/PJ_natearth2.c
+++ b/src/PJ_natearth2.c
@@ -26,6 +26,8 @@ PROJ_HEAD(natearth2, "Natural Earth 2") "\n\tPCyl., Sph.";
 #define C3 (13 * B3)
 #define EPS 1e-11
 #define MAX_Y (0.84719 * 0.535117535153096 * M_PI)
+/* Not sure at all of the appropriate number for MAX_ITER... */
+#define MAX_ITER 100
 
 
 static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
@@ -46,6 +48,7 @@ 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};
     double yc, tol, y2, y4, y6, f, fder;
+    int i;
     (void) P;
 
     /* make sure y is inside valid range */
@@ -57,7 +60,7 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
 
     /* latitude */
     yc = xy.y;
-    for (;;) { /* Newton-Raphson */
+    for (i = MAX_ITER; i ; --i) { /* Newton-Raphson */
         y2 = yc * yc;
         y4 = y2 * y2;
         f = (yc * (B0 + y4 * y4 * (B1 + B2 * y2 + B3 * y4))) - xy.y;
@@ -67,6 +70,8 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
             break;
         }
     }
+    if( i == 0 )
+        pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );
     lp.phi = yc;
 
     /* longitude */
diff --git a/src/PJ_patterson.c b/src/PJ_patterson.c
index 50ca0981..952b63ea 100644
--- a/src/PJ_patterson.c
+++ b/src/PJ_patterson.c
@@ -53,6 +53,8 @@ PROJ_HEAD(patterson, "Patterson Cylindrical") "\n\tCyl.";
 #define C4 (9.0 * K4)
 #define EPS11 1.0e-11
 #define MAX_Y 1.790857183
+/* Not sure at all of the appropriate number for MAX_ITER... */
+#define MAX_ITER 100
 
 
 static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
@@ -71,6 +73,7 @@ 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};
     double yc, tol, y2, f, fder;
+    int i;
     (void) P;
 
     yc = xy.y;
@@ -82,7 +85,7 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
         xy.y = -MAX_Y;
     }
 
-    for (;;) { /* Newton-Raphson */
+    for (i = MAX_ITER; i ; --i) { /* Newton-Raphson */
         y2 = yc * yc;
         f = (yc * (K1 + y2 * y2 * (K2 + y2 * (K3 + K4 * y2)))) - xy.y;
         fder = C1 + y2 * y2 * (C2 + y2 * (C3 + C4 * y2));
@@ -91,6 +94,8 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
             break;
         }
     }
+    if( i == 0 )
+        pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );
     lp.phi = yc;
 
     /* longitude */
diff --git a/src/PJ_robin.c b/src/PJ_robin.c
index 93330072..4eab552a 100644
--- a/src/PJ_robin.c
+++ b/src/PJ_robin.c
@@ -70,7 +70,8 @@ static const struct COEFS Y[] = {
 #define NODES   18
 #define ONEEPS  1.000001
 #define EPS 1e-8
-
+/* Not sure at all of the appropriate number for MAX_ITER... */
+#define MAX_ITER 100
 
 static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
     XY xy = {0.0,0.0};
@@ -96,6 +97,7 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
     int i;
     double t, t1;
     struct COEFS T;
+    int iters;
 
     lp.lam = xy.x / FXC;
     lp.phi = fabs(xy.y / FYC);
@@ -120,11 +122,13 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
         t = 5. * (lp.phi - T.c0)/(Y[i+1].c0 - T.c0);
         /* make into root */
         T.c0 = (float)(T.c0 - lp.phi);
-        for (;;) { /* Newton-Raphson reduction */
+        for (iters = MAX_ITER; iters ; --iters) { /* Newton-Raphson */
             t -= t1 = V(T,t) / DV(T,t);
             if (fabs(t1) < EPS)
                 break;
         }
+        if( iters == 0 )
+            pj_ctx_set_errno( P->ctx, PJD_ERR_NON_CONVERGENT );
         lp.phi = (5 * i + t) * DEG_TO_RAD;
         if (xy.y < 0.) lp.phi = -lp.phi;
         lp.lam /= V(X[i], t);