summaryrefslogtreecommitdiff
path: root/lib/Math/SpecFun/Erf.pm
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Math/SpecFun/Erf.pm')
-rw-r--r--lib/Math/SpecFun/Erf.pm565
1 files changed, 0 insertions, 565 deletions
diff --git a/lib/Math/SpecFun/Erf.pm b/lib/Math/SpecFun/Erf.pm
deleted file mode 100644
index c913792..0000000
--- a/lib/Math/SpecFun/Erf.pm
+++ /dev/null
@@ -1,565 +0,0 @@
-#
-# 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;
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-02 01:32:39 +0000