sub erfinv {
    
    my $y = @_ ? $_[0] : $_;

    return  0             if $y == 0;
    return  erfcinv(1-$y) if $y >  0.5;
    return -erfcinv(1+$y) if $y < -0.5;

    #
    # Halley's rational 3rd order method:
    #   u <- f(x)/f'(x)
    #   v <- f''(x)/f'(x)
    #   x <- x - u/(1-u*v/2)
    #
    # Here:
    #   f(x) = erf(x) - y
    #   f'(x) = 2/sqrt(pi)*exp(-x*x)
    #   f''(x) = -4/sqrt(pi)*x*exp(-x*x)
    #
    my $x = 0;
    my $dx;
    my $c = .88622692545275801364908374167055;  # sqrt(pi)/2
    my $eps = 5e-20;
    do {
        my $f = erf($x) - $y;
        my $u = $c*$f*exp($x*$x);
        print "u: $u";
       # print "u: $u";
        $dx = -$u/(1+$u*$x);
        $x += $dx;
    } until abs($dx/$x) <= $eps;
    return $x;
}