diff --git a/applications/test/incGamma/Test-incGamma.C b/applications/test/incGamma/Test-incGamma.C
index ebaf50bfd9..71f3bf4be4 100644
--- a/applications/test/incGamma/Test-incGamma.C
+++ b/applications/test/incGamma/Test-incGamma.C
@@ -51,6 +51,8 @@ scalar grade(const scalar f)
 
 int main(int argc, char *argv[])
 {
+    argList args(argc, argv);
+
     const scalarField as({0.25, 0.5, 1, 2, 4});
 
     static const scalar nXs = 10000;
@@ -92,7 +94,15 @@ int main(int argc, char *argv[])
         {
             const scalar P = incGammaRatio_P(as[ai], xs[xi]);
             const scalar PStar =
-                incGammaRatio_P(as[ai], invIncGammaRatio_P(as[ai], P));
+                incGammaRatio_P
+                (
+                    as[ai],
+                    invIncGammaRatio_P
+                    (
+                        as[ai],
+                        min(max(P, small), 1 - small)
+                    )
+                );
             scalarPs[ai][xi] = P;
             scalarPErrors[ai][xi] = mag(P - PStar);
         }
diff --git a/src/OpenFOAM/primitives/Scalar/scalar/invIncGamma.C b/src/OpenFOAM/primitives/Scalar/scalar/invIncGamma.C
index 104ee97bbd..2213100709 100644
--- a/src/OpenFOAM/primitives/Scalar/scalar/invIncGamma.C
+++ b/src/OpenFOAM/primitives/Scalar/scalar/invIncGamma.C
@@ -349,8 +349,8 @@ Foam::scalar Foam::invIncGammaRatio_P(const scalar a, const scalar P)
 
     for (label iter = 0; iter < nIter; ++ iter)
     {
-        const scalar dP = incGammaRatio_P(a, x) - P;
-        const scalar t = dP/R(a, x);
+        const scalar dP = incGammaRatio_P(a, max(x, vSmall)) - P;
+        const scalar t = dP/max(R(a, x), vSmall);
         const scalar w = (a - 1 - x)/2;
 
         bool halley = mag(t) < 0.1 && mag(w*t) < 0.1;