package hep.aida.ref.pdf;

/**
 *
 * @author turri
 */
public class MathUtils {


    public static double erf(double d) {
        double d1 = (0.0D / 0.0D);
        double d2 = Math.abs(d);
        if(d2 <= 1.0D) {
            if(d2 <= erf_sqeps)
                d1 = (2D * d) / 1.7724538509055161D;
            else
                d1 = d * (1.0D + csevl(2D * d * d - 1.0D, erfc_coef));
        } else
            if(d2 <= erf_xbig)
                d1 = sign(1.0D - erfc(d2), d);
            else
                d1 = sign(1.0D, d);
        return d1;
    }


    public static double erfc(double d) {
        double d1 = (0.0D / 0.0D);
        if(d <= erfc_xsml)
            d1 = 2D;
        else
            if(d <= erfc_xmax) {
                double d2 = Math.abs(d);
                if(d2 <= 1.0D) {
                    if(d2 < erf_sqeps)
                        d1 = 1.0D - (2D * d) / 1.7724538509055161D;
                    else
                        d1 = 1.0D - d * (1.0D + csevl(2D * d * d - 1.0D, erfc_coef));
                } else {
                    d2 *= d2;
                    if(d2 <= 4D)
                        d1 = (Math.exp(-d2) / Math.abs(d)) * (0.5D + csevl((8D / d2 - 5D) / 3D, erfc2_coef));
                    else
                        d1 = (Math.exp(-d2) / Math.abs(d)) * (0.5D + csevl(8D / d2 - 1.0D, erfcc_coef));
                    if(d < 0.0D)
                        d1 = 2D - d1;
                }
            } else {
                d1 = 0.0D;
            }
        return d1;
    }

    private static double csevl(double d, double ad[]) {
        int i = ad.length;
        double d2 = 0.0D;
        double d1 = 0.0D;
        double d3 = 0.0D;
        double d4 = 2D * d;
        for(int j = i - 1; j >= 0; j--) {
            d3 = d2;
            d2 = d1;
            d1 = (d4 * d2 - d3) + ad[j];
        }

        return 0.5D * (d1 - d3);
    }

    public static double sign(double d, double d1) {
        double d2 = d >= 0.0D ? d : -d;
        return d1 >= 0.0D ? d2 : -d2;
    }

    private static final double erfc_coef[] = {
        -0.049046121234691806D, -0.14226120510371365D, 0.010035582187599796D, -0.00057687646997674853D,
	2.741993125219606E-005D, -1.1043175507344507E-006D, 3.8488755420345036E-008D, -1.1808582533875466E-009D,
	3.2334215826050907E-011D, -7.9910159470045487E-013D,
        1.7990725113961456E-014D, -3.7186354878186928E-016D, 7.1035990037142532E-018D
    };
    private static final double erfc2_coef[] = {
        -0.069601346602309502D, -0.041101339362620892D, 0.0039144958666896268D,
	-0.00049063956505489791D, 7.1574790013770361E-005D, -1.1530716341312328E-005D,
	1.9946705902019974E-006D, -3.6426664715992229E-007D, 6.9443726100050124E-008D,
	-1.3712209021043659E-008D, 2.7883896610071373E-009D, -5.8141647243311614E-010D,
	1.2389204917527532E-010D, -2.6906391453067435E-011D, 5.9426143508479114E-012D,
	-1.3323867357581197E-012D, 3.0280468061771323E-013D, -6.9666488149410327E-014D,
	1.620854541053923E-014D, -3.8099344652504917E-015D, 9.0404878159788311E-016D,
	-2.1640061950896072E-016D, 5.2221022339958551E-017D, -1.2697296023645554E-017D,
	3.1091455042761977E-018D
    };
    private static final double erfcc_coef[] = {
        0.071517931020292483D, -0.026532434337606717D, 0.0017111539779208558D,
	-0.00016375166345851787D, 1.9871293500552038E-005D, -2.8437124127665552E-006D,
	4.6061613089631305E-007D, -8.227753025879209E-008D, 1.5921418727709012E-008D,
	-3.2950713622528431E-009D, 7.2234397604005558E-010D, -1.6648558133987297E-010D,
	4.0103925882376649E-011D, -1.0048162144257311E-011D, 2.6082759133003339E-012D,
	-6.9911105604040245E-013D, 1.9294923332617072E-013D, -5.4701311887543309E-014D,
	1.5896633097626975E-014D, -4.7268939801975551E-015D, 1.4358733767849847E-015D,
	-4.4495105618173579E-016D, 1.4048108847682335E-016D, -4.5138183877642106E-017D,
	1.4745215410451331E-017D, -4.8926214069457765E-018D, 1.6476121414106467E-018D,
	-5.6268171763294081E-019D, 1.9474433822320786E-019D
    };
    private static final double erf_xbig = Math.sqrt(-Math.log(1.9678190753608168E-016D));
    private static final double erf_sqeps = Math.sqrt(2.2204460492503E-016D);
    private static final double erfc_xsml = -Math.sqrt(-Math.log(1.9678190753608168E-016D));
    private static double erfc_xmax;


}