#
# Author:      Peter J. Acklam
# Time-stamp:  2000-11-29 23:04:53
# E-mail:      pjacklam@online.no
# URL:         http://home.online.no/~pjacklam

=head1 NAME

Math::SpecFun::Erf - error and scaled and unscaled complementary error
functions and their inverses

=head1 SYNOPSIS

    use Math::SpecFun::Erf qw(erf erfc erfcx erfinv erfcinv erfcxinv);

Imports all the routines explicitly.  Use a subset of the list for the
routines you want.

    use Math::SpecFun::Erf qw(:all);

Imports all the routines, as well.

=head1 DESCRIPTION

This module implements the error function, C<erf>, and its inverse
C<erfinv>, the complementary error function, C<erfc>, and its inverse
C<erfcinv>, and the scaled complementary error function, C<erfcx>, and its
inverse C<erfcxinv>.

For references and details about the algorithms, see the comments inside
this module.

=head1 FUNCTIONS

=over 8

=item erf EXPR

=item erf

Returns the error function evaluated at EXPR.  If EXPR is omitted, C<$_> is
used.  The error function is

    erf(x) = 2/sqrt(PI) * integral from 0 to x of exp(-t*t) dt

=item erfinv EXPR

=item erfinv

Returns the inverse of the error function evaluated at EXPR.  If EXPR is
omitted, C<$_> is used.

=item erfc EXPR

=item erfc

Returns the complementary error function evaluated at EXPR.  If EXPR is
omitted, C<$_> is used.  The complementary error function is

    erfc(x) = 2/sqrt(PI) * integral from x to infinity of exp(-t*t) dt
            = 1 - erf(x)

Here is a function returning the lower tail probability of the standard
normal distribution function

    use Math::SpecFun::Erf qw(erfc);

    sub ltpnorm ($) {
        erfc( - $_[0] / sqrt(2) )/2;
    }

=item erfcinv EXPR

=item erfcinv

Returns the inverse complementary error function evaluated at EXPR.  If EXPR
is omitted, C<$_> is used.

Here is a function returning the lower tail quantile of the standard normal
distribution function

    use Math::SpecFun::Erf qw(erfcinv);

    sub ltqnorm ($) {
        -sqrt(2) * erfcinv( 2 * $_[0] );
    }

=item erfcx EXPR

=item erfcx

Returns the scaled complementary error function evaluated at EXPR.  If EXPR
is omitted, C<$_> is used.  The scaled complementary error function is

    erfcx(x) = exp(x*x) * erfc(x)

=item erfcxinv EXPR

=item erfcxinv

Returns the inverse scaled complementary error function evaluated at EXPR.
If EXPR is omitted, C<$_> is used.

=back

=head1 HISTORY

=over 4

=item Version 0.03

Added the inverse functions.

=item Version 0.02

Minor code tweaking.

=item Version 0.01

First release.

=back

=head1 AUTHOR

Perl translation by Peter J. Acklam E<lt>pjacklam@online.noE<gt>

FORTRAN code by W. J. Cody, Argonne National Laboratory, March 19, 1990.
FORTRAN code can be found at http://www.netlib.org/specfun/erf

=head1 COPYRIGHT

Copyright (c) 1999-2000 Peter J. Acklam.  All rights reserved.
This program is free software; you can redistribute it and/or
modify it under the same terms as Perl itself.

=cut

package Math::SpecFun::Erf;
require 5.000;
require Exporter;

use strict;
use vars qw($VERSION @ISA @EXPORT_OK %EXPORT_TAGS);

$VERSION = '0.02';
@ISA = qw(Exporter);
@EXPORT_OK = qw(erf erfc erfcx erfinv erfcinv erfcxinv);
%EXPORT_TAGS = ( all => [ @EXPORT_OK ] );

########################################################################
## Internal functions.
########################################################################

sub calerf {
    my ($arg, $result, $jint) = @_;
    local $[ = 1;
#------------------------------------------------------------------
#
# This packet evaluates  erf(x),  erfc(x),  and  exp(x*x)*erfc(x)
#   for a real argument  x.  It contains three FUNCTION type
#   subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX),
#   and one SUBROUTINE type subprogram, CALERF.  The calling
#   statements for the primary entries are:
#
#                   Y=ERF(X)     (or   Y=DERF(X)),
#
#                   Y=ERFC(X)    (or   Y=DERFC(X)),
#   and
#                   Y=ERFCX(X)   (or   Y=DERFCX(X)).
#
#   The routine  CALERF  is intended for internal packet use only,
#   all computations within the packet being concentrated in this
#   routine.  The function subprograms invoke  CALERF  with the
#   statement
#
#          CALL CALERF(ARG,RESULT,JINT)
#
#   where the parameter usage is as follows
#
#      Function                     Parameters for CALERF
#       call              ARG                  Result          JINT
#
#     ERF(ARG)      ANY REAL ARGUMENT         ERF(ARG)          0
#     ERFC(ARG)     ABS(ARG) < XBIG           ERFC(ARG)         1
#     ERFCX(ARG)    XNEG < ARG < XMAX         ERFCX(ARG)        2
#
#   The main computation evaluates near-minimax approximations
#   from "Rational Chebyshev approximations for the error function"
#   by W. J. Cody, Math. Comp., 1969, PP. 631-638.  This
#   transportable program uses rational functions that theoretically
#   approximate  erf(x)  and  erfc(x)  to at least 18 significant
#   decimal digits.  The accuracy achieved depends on the arithmetic
#   system, the compiler, the intrinsic functions, and proper
#   selection of the machine-dependent constants.
#
#*******************************************************************
#*******************************************************************
#
# Explanation of machine-dependent constants
#
#   XMIN   = the smallest positive floating-point number.
#   XINF   = the largest positive finite floating-point number.
#   XNEG   = the largest negative argument acceptable to ERFCX;
#            the negative of the solution to the equation
#            2*exp(x*x) = XINF.
#   XSMALL = argument below which erf(x) may be represented by
#            2*x/sqrt(pi)  and above which  x*x  will not underflow.
#            A conservative value is the largest machine number X
#            such that   1.0 + X = 1.0   to machine precision.
#   XBIG   = largest argument acceptable to ERFC;  solution to
#            the equation:  W(x) * (1-0.5/x**2) = XMIN,  where
#            W(x) = exp(-x*x)/[x*sqrt(pi)].
#   XHUGE  = argument above which  1.0 - 1/(2*x*x) = 1.0  to
#            machine precision.  A conservative value is
#            1/[2*sqrt(XSMALL)]
#   XMAX   = largest acceptable argument to ERFCX; the minimum
#            of XINF and 1/[sqrt(pi)*XMIN].
#
#   Approximate values for some important machines are:
#
#                          XMIN       XINF        XNEG     XSMALL
#
#  CDC 7600      (S.P.)  3.13E-294   1.26E+322   -27.220  7.11E-15
#  CRAY-1        (S.P.)  4.58E-2467  5.45E+2465  -75.345  7.11E-15
#  IEEE (IBM/XT,
#    SUN, etc.)  (S.P.)  1.18E-38    3.40E+38     -9.382  5.96E-8
#  IEEE (IBM/XT,
#    SUN, etc.)  (D.P.)  2.23D-308   1.79D+308   -26.628  1.11D-16
#  IBM 195       (D.P.)  5.40D-79    7.23E+75    -13.190  1.39D-17
#  UNIVAC 1108   (D.P.)  2.78D-309   8.98D+307   -26.615  1.73D-18
#  VAX D-Format  (D.P.)  2.94D-39    1.70D+38     -9.345  1.39D-17
#  VAX G-Format  (D.P.)  5.56D-309   8.98D+307   -26.615  1.11D-16
#
#
#                          XBIG       XHUGE       XMAX
#
#  CDC 7600      (S.P.)  25.922      8.39E+6     1.80X+293
#  CRAY-1        (S.P.)  75.326      8.39E+6     5.45E+2465
#  IEEE (IBM/XT,
#    SUN, etc.)  (S.P.)   9.194      2.90E+3     4.79E+37
#  IEEE (IBM/XT,
#    SUN, etc.)  (D.P.)  26.543      6.71D+7     2.53D+307
#  IBM 195       (D.P.)  13.306      1.90D+8     7.23E+75
#  UNIVAC 1108   (D.P.)  26.582      5.37D+8     8.98D+307
#  VAX D-Format  (D.P.)   9.269      1.90D+8     1.70D+38
#  VAX G-Format  (D.P.)  26.569      6.71D+7     8.98D+307
#
#*******************************************************************
#*******************************************************************
#
# Error returns
#
#  The program returns  ERFC = 0      for  ARG >= XBIG;
#
#                       ERFCX = XINF  for  ARG < XNEG;
#      and
#                       ERFCX = 0     for  ARG >= XMAX.
#
#
# Intrinsic functions required are:
#
#     ABS, AINT, EXP
#
#
#  Author: W. J. Cody
#          Mathematics and Computer Science Division
#          Argonne National Laboratory
#          Argonne, IL 60439
#
#  Latest modification: March 19, 1990
#
# Translation to Perl by Peter J. Acklam, December 3, 1999
#
#------------------------------------------------------------------
    my ($i);
    my ($x, $del, $xden, $xnum, $y, $ysq);
#------------------------------------------------------------------
#  Mathematical constants
#------------------------------------------------------------------
    my ($four, $one, $half, $two, $zero) = (4, 1, 0.5, 2, 0);
    my $sqrpi = 5.6418958354775628695e-1;
    my $thresh = 0.46875;
    my $sixten = 16;
#------------------------------------------------------------------
#  Machine-dependent constants
#------------------------------------------------------------------
    my ($xinf, $xneg, $xsmall) = (1.79e308, -26.628, 1.11e-16);
    my ($xbig, $xhuge, $xmax) = (26.543, 6.71e7, 2.53e307);
#------------------------------------------------------------------
#  Coefficients for approximation to  erf  in first interval
#------------------------------------------------------------------
    my @a = (3.16112374387056560e00, 1.13864154151050156e02,
             3.77485237685302021e02, 3.20937758913846947e03,
             1.85777706184603153e-1);
    my @b = (2.36012909523441209e01, 2.44024637934444173e02,
             1.28261652607737228e03, 2.84423683343917062e03);
#------------------------------------------------------------------
#  Coefficients for approximation to  erfc  in second interval
#------------------------------------------------------------------
    my @c = (5.64188496988670089e-1, 8.88314979438837594e0,
             6.61191906371416295e01, 2.98635138197400131e02,
             8.81952221241769090e02, 1.71204761263407058e03,
             2.05107837782607147e03, 1.23033935479799725e03,
             2.15311535474403846e-8);
    my @d = (1.57449261107098347e01, 1.17693950891312499e02,
             5.37181101862009858e02, 1.62138957456669019e03,
             3.29079923573345963e03, 4.36261909014324716e03,
             3.43936767414372164e03, 1.23033935480374942e03);
#------------------------------------------------------------------
#  Coefficients for approximation to  erfc  in third interval
#------------------------------------------------------------------
    my @p = (3.05326634961232344e-1, 3.60344899949804439e-1,
             1.25781726111229246e-1, 1.60837851487422766e-2,
             6.58749161529837803e-4, 1.63153871373020978e-2);
    my @q = (2.56852019228982242e00, 1.87295284992346047e00,
             5.27905102951428412e-1, 6.05183413124413191e-2,
             2.33520497626869185e-3);
#------------------------------------------------------------------
    $x = $arg;
    $y = abs($x);
    if ($y <= $thresh) {
#------------------------------------------------------------------
#  Evaluate  erf  for  |X| <= 0.46875
#------------------------------------------------------------------
        $ysq = $zero;
        if ($y > $xsmall) { $ysq = $y * $y }
        $xnum = $a[5]*$ysq;
        $xden = $ysq;
        for (my $i = 1 ; $i <= 3 ; ++$i) {
            $xnum = ($xnum + $a[$i]) * $ysq;
            $xden = ($xden + $b[$i]) * $ysq;
        }
        $$result = $x * ($xnum + $a[4]) / ($xden + $b[4]);
        if ($jint != 0) { $$result = $one - $$result }
        if ($jint == 2) { $$result = exp($ysq) * $$result }
        goto x800;
#------------------------------------------------------------------
#  Evaluate  erfc  for 0.46875 <= |X| <= 4.0
#------------------------------------------------------------------
    } elsif ($y <= $four) {
        $xnum = $c[9]*$y;
        $xden = $y;
        for (my $i = 1 ; $i <= 7 ; ++$i) {
            $xnum = ($xnum + $c[$i]) * $y;
            $xden = ($xden + $d[$i]) * $y;
        }
        $$result = ($xnum + $c[8]) / ($xden + $d[8]);
        if ($jint != 2) {
            $ysq = int($y*$sixten)/$sixten;
            $del = ($y-$ysq)*($y+$ysq);
            $$result = exp(-$ysq*$ysq) * exp(-$del) * $$result;
        }
#------------------------------------------------------------------
#  Evaluate  erfc  for |X| > 4.0
#------------------------------------------------------------------
    } else {
        $$result = $zero;
        if ($y >= $xbig) {
            if (($jint != 2) || ($y >= $xmax)) { goto x300 }
            if ($y >= $xhuge) {
                $$result = $sqrpi / $y;
                goto x300;
            }
        }
        $ysq = $one / ($y * $y);
        $xnum = $p[6]*$ysq;
        $xden = $ysq;
        for (my $i = 1 ; $i <= 4 ; ++$i) {
            $xnum = ($xnum + $p[$i]) * $ysq;
            $xden = ($xden + $q[$i]) * $ysq;
        }
        $$result = $ysq *($xnum + $p[5]) / ($xden + $q[5]);
        $$result = ($sqrpi -  $$result) / $y;
        if ($jint != 2) {
            $ysq = int($y*$sixten)/$sixten;
            $del = ($y-$ysq)*($y+$ysq);
            $$result = exp(-$ysq*$ysq) * exp(-$del) * $$result;
        }
    }
#------------------------------------------------------------------
#  Fix up for negative argument, erf, etc.
#------------------------------------------------------------------
  x300:
    if ($jint == 0) {
        $$result = ($half - $$result) + $half;
        if ($x < $zero) { $$result = -$$result }
    } elsif ($jint == 1) {
        if ($x < $zero) { $$result = $two - $$result }
    } else {
        if ($x < $zero) {
            if ($x < $xneg) {
                $$result = $xinf;
            } else {
                $ysq = int($x*$sixten)/$sixten;
                $del = ($x-$ysq)*($x+$ysq);
                $y = exp($ysq*$ysq) * exp($del);
                $$result = ($y+$y) - $$result;
            }
        }
    }
  x800:
    return 1;
#---------- Last card of CALERF ----------
}

sub erf {
    my $x = @_ ? $_[0] : $_;
#--------------------------------------------------------------------
#
# This subprogram computes approximate values for erf(x).
#   (see comments heading CALERF).
#
#   Author/date: W. J. Cody, January 8, 1985
#
# Translation to Perl by Peter J. Acklam, December 3, 1999
#
#--------------------------------------------------------------------
    my $result;
    my $jint = 0;
    calerf($x, \$result, $jint);
    my $erf = $result;
    return $erf;
#---------- Last card of ERF ----------
}

########################################################################
## User functions.
########################################################################

sub erfc {
    my $x = @_ ? $_[0] : $_;
#--------------------------------------------------------------------
#
# This subprogram computes approximate values for erfc(x).
#   (see comments heading CALERF).
#
#   Author/date: W. J. Cody, January 8, 1985
#
# Translation to Perl by Peter J. Acklam, December 3, 1999
#
#--------------------------------------------------------------------
    my ($result);
    my $jint = 1;
    calerf($x, \$result, $jint);
    my $erfc = $result;
    return $erfc;
#---------- Last card of ERFC ----------
}

sub erfcx {
    my $x = @_ ? $_[0] : $_;
#------------------------------------------------------------------
#
# This subprogram computes approximate values for exp(x*x) * erfc(x).
#   (see comments heading CALERF).
#
#   Author/date: W. J. Cody, March 30, 1987
#
# Translation to Perl by Peter J. Acklam, December 3, 1999
#
#------------------------------------------------------------------
    my ($result);
    my $jint = 2;
    calerf($x, \$result, $jint);
    my $erfcx = $result;
    return $erfcx;
#---------- Last card of ERFCX ----------
}

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-15;
    do {
        my $f = erf($x) - $y;
        my $u = $c*$f*exp($x*$x);
        $dx = -$u/(1+$u*$x);
        $x += $dx;
    } until abs($dx/$x) <= $eps;
    return $x;
}

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

    return 0 if $y == 1;

    my $flag = $y > 1;
    $y = 2 - $y if $flag;

    #
    # 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) = erfc(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-15;
    do {
        my $u = $c*(erfcx($x) - $y*exp($x*$x));
        $dx = -$u/(1+$u*$x);
        $x += $dx;
    } until abs($dx/$x) <= $eps;

    return $flag ? -$x : $x;
}

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

    return 0 if $y == 1;

    #
    # Halley's 3rd order method:
    #   u <- f(x)/f'(x)
    #   v <- f''(x)/f'(x)
    #   x <- x - u/(1-u*v/2)
    #
    # Here:
    #   f(x) = erfcx(x) - y
    #   f'(x) = 2*(x*erfcx(x)-1/sqrt(pi));
    #   f''(x) = (2+4*x*x)*erfcx(x) - 4*x/sqrt(pi);
    #
    my $x = 0;
    my $dx;
    my $c = .56418958354775628694807945156079;  # 1/sqrt(pi)
    my $d = 2.2567583341910251477923178062432;  # 4/sqrt(pi)
    my $eps = 5e-15;
    do {
        my $f = erfcx($x) - $y;
        my $df = 2*($x*erfcx($x)-$c);
        my $ddf = (2+4*$x*$x)*erfcx($x) - $x*$d;
        my $u = $f/$df;
        my $v = $ddf/$df;
        $dx = -$u/(1-$u*$v/2);
        $x += $dx;
    } until abs($dx/$x) <= $eps;
    return $x;
}