[MATH-456] Erf.erf - handle infinities and large values

• Type: Bug
• Status: Closed
• Priority: Minor
• Resolution: Fixed
• Affects Version/s: 2.2
• Fix Version/s: 2.2
• Assignee: Unassigned
• Reporter: Josh Milthorpe
• Created: 09/Dec/10 11:46
• Updated: 23/Mar/11 20:29
• Resolved: 01/Jan/11 01:24

Description

Erf.erf(double) crashes when presented with infinities or large values, as follows:

org.apache.commons.math.ConvergenceException: Continued fraction diverged to NaN for value ∞
	at org.apache.commons.math.util.ContinuedFraction.evaluate(ContinuedFraction.java:186)
	at org.apache.commons.math.special.Gamma.regularizedGammaQ(Gamma.java:266)
	at org.apache.commons.math.special.Gamma.regularizedGammaP(Gamma.java:173)
	at org.apache.commons.math.special.Erf.erf(Erf.java:56)
	at TestInfErf.main(TestInfErf.java:9)

The following code demonstrates this crash:

import org.apache.commons.math.MathException;
import org.apache.commons.math.special.Erf;

public class TestInfErf {
    public static void main(String[] args) {
        try {
            System.out.println("erf(Inf) = " + Erf.erf(Double.POSITIVE_INFINITY));
            System.out.println("erf(-Inf) = " + Erf.erf(Double.NEGATIVE_INFINITY));
            System.out.println("erf(Huge) = " + Erf.erf(1e300));
        } catch (MathException e) { 
            e.printStackTrace(); 
        }
    }
}

At double precision, erf(x) = 1.0 for x > 6.0 and erf(x) = -1.0 for x < -6.0. Therefore Erf.java could be patched as follows:

Index: src/main/java/org/apache/commons/math/special/Erf.java
===================================================================
--- src/main/java/org/apache/commons/math/special/Erf.java	(revision 1043888)
+++ src/main/java/org/apache/commons/math/special/Erf.java	(working copy)
@@ -48,6 +48,12 @@
      * @throws MathException if the algorithm fails to converge.
      */
     public static double erf(double x) throws MathException {
+        // at double precision, erf(x) = (+/-)1.0 for |x| > 6.0
+        if (x > 6.0) {
+            return 1.0;
+        } else if (x <-6.0) {
+            return -1.0;
+        }
         double ret = Gamma.regularizedGammaP(0.5, x * x, 1.0e-15, 10000);
         if (x < 0) {
             ret = -ret;

Issue Links

• is blocked by: MATH-364 Make Erf more precise in the tails by providing erfc
• is related to: MATH-301 Erf(z) should return 1.0 for z 'large' but fails with a MaxIterationsExceededException for z > 26.0.