MAINT: special: avoid iteration in Wright Omega for large negatives

Just return `exp(x)`, which is accurate to double precision for
roughly `x < -50`.

scipy/special/_precompute/wrightomega.py

@@ -19,11 +19,25 @@ def wrightomega_series_error(x):
     return abs(series - desired) / desired
 
 
+def wrightomega_exp_error(x):
+    exponential_approx = mpmath.exp(x)
+    desired = mpmath_wrightomega(x)
+    return abs(exponential_approx - desired) / desired
+
+
 def main():
     desired_error = 2 * np.finfo(float).eps
+    print('Series Error')
     for x in [1e5, 1e10, 1e15, 1e20]:
-        error = wrightomega_series_error(x)
-        print(error, error < desired_error)
+        with mpmath.workdps(100):
+            error = wrightomega_series_error(x)
+        print(x, error, error < desired_error)
+
+    print('Exp error')
+    for x in [-10, -25, -50, -100, -200, -400, -700, -740]:
+        with mpmath.workdps(100):
+            error = wrightomega_exp_error(x)
+        print(x, error, error < desired_error)
 
 
 if __name__ == '__main__':

scipy/special/tests/test_mpmath.py

@@ -1860,7 +1860,8 @@ def mpmath_wrightomega_real(x):
             sc.wrightomega,
             mpmath_wrightomega_real,
             [Arg()],
-            rtol=1e-14,
+            rtol=1e-15,
+            atol=0,
             nan_ok=False,
         )
 

scipy/special/tests/test_wrightomega.py

@@ -90,6 +90,28 @@ def test_wrightomega_real_series_crossover():
     )
 
 
+def test_wrightomega_exp_approximation_crossover():
+    desired_error = 2 * np.finfo(float).eps
+    crossover = -50
+    x_before_crossover = np.nextafter(crossover, np.inf)
+    x_after_crossover = np.nextafter(crossover, -np.inf)
+    # Computed using Mpmath
+    desired_before_crossover = 1.9287498479639314876e-22
+    desired_after_crossover = 1.9287498479639040784e-22
+    assert_allclose(
+        sc.wrightomega(x_before_crossover),
+        desired_before_crossover,
+        atol=0,
+        rtol=desired_error,
+    )
+    assert_allclose(
+        sc.wrightomega(x_after_crossover),
+        desired_after_crossover,
+        atol=0,
+        rtol=desired_error,
+    )
+
+
 def test_wrightomega_real_versus_complex():
     x = np.linspace(-500, 500, 1001)
     assert_allclose(

scipy/special/wright.cc

@@ -391,31 +391,39 @@ wright::wrightomega_real(double x)
 	}
     }
 
-  /* Split int three distinct intervals (-inf,-2), [-2,1), [1,inf) */
-  if (x < -2.0)
+  if (x < -50.0)
     {
-      /* exponential is approx < 1.3e-1 accurate */
+      /*
+       * Skip the iterative scheme because it exp(x) is already
+       * accurate to double precision.
+       */
       w = exp(x);
       if (w == 0.0)
         {
-          /* Skip the iterative scheme because it computes log(w) */
           sf_error("wrightomega", SF_ERROR_UNDERFLOW, "underflow in exponential series");
-          return w;
         }
+      return w;
     }
-  else if (x < 1)
-    {
-      /* on [-2,1) approx < 1.5e-1 accurate */
-      w = exp(2.0*(x-1.0)/3.0);
-    }
-  else if (x > 1e20)
+  if (x > 1e20)
     {
       /*
        * Skip the iterative scheme because the result is just x to
        * double precision
        */
       return x;
     }
+
+  /* Split int three distinct intervals (-inf,-2), [-2,1), [1,inf) */
+  if (x < -2.0)
+    {
+      /* exponential is approx < 1.3e-1 accurate */
+      w = exp(x);
+    }
+  else if (x < 1)
+    {
+      /* on [-2,1) approx < 1.5e-1 accurate */
+      w = exp(2.0*(x-1.0)/3.0);
+    }
   else
     {
       /* infinite series with 2 terms approx <1.7e-1 accurate */