package ij.util;

public class IJMath {

	static double A = 0.05857895078654250866288;
	static double B = -0.00626245895772819579;
	static double C = -0.00299946450696036814;
	static double D = 0.289389696496082416;
	static double E = 0.0539962589851632982;
	static double F = 0.00508516909930653109;
	static double G = 0.000215969713046142876;
	static double H = -0.000225663858340491571;
	static double I = -3.06833213472529049e-7;

	/* This approximation of the error function erf has maximum absolute and relative errors below 1e-12
	 * Except for low and high x, it is based on the type erf(x) = sgn(x) * sqrt(1-exp(-x^2*f(x)),
	 * with the function f(x) having a smooth transition from 4/pi = 1.273... at x=0 to 1 at high x */
	public static double erf(double x) {
		double x2 = sqr(x);
		if (x2 < 1e-8)
			return (2/Math.sqrt(Math.PI))*x*(1+x2*(-1./3.+x2*(1./10.)));  // Taylor series for low x
		double erf = x2 > 36 ? 1 :  //the polynomials go crazy for some large x2; erf(6) is 1 - 2e-17, less than ulp from 1
				Math.sqrt(1 - Math.exp(-sqr(x*((2/Math.sqrt(Math.PI) - A) + A *
				(1 + x2*x2*(B + x2*(C + x2*(H + x2*I)))) /
				(1 + x2*(D + x2*(E + x2*(F + x2*G))))
				))));
		return x>0 ? erf : -erf; //could be Math.copySign in Java 1.6 & up
	}

	static double sqr(double x) {
		return x*x;
	}
	
}