Math-Cephes
view release on metacpan or search on metacpan
lib/Math/Cephes.pod view on Meta::CPAN
ndtri: Inverse of Normal distribution function
pdtr: Poisson distribution
pdtrc: Complemented poisson distribution
pdtri: Inverse Poisson distribution
stdtr: Student's t distribution
stdtri: Functional inverse of Student's t distribution
=item use Math::Cephes qw(:gammas);
imports
fac: Factorial function
gamma: Gamma function
igam: Incomplete gamma integral
igamc: Complemented incomplete gamma integral
igami: Inverse of complemented imcomplete gamma integral
psi: Psi (digamma) function
rgamma: Reciprocal gamma function
=item use Math::Cephes qw(:betas);
imports
beta: Beta function
incbet: Incomplete beta integral
incbi: Inverse of imcomplete beta integral
lbeta: Natural logarithm of |beta|
=item use Math::Cephes qw(:elliptics);
imports
ellie: Incomplete elliptic integral of the second kind
ellik: Incomplete elliptic integral of the first kind
ellpe: Complete elliptic integral of the second kind
ellpj: Jacobian Elliptic Functions
ellpk: Complete elliptic integral of the first kind
=item use Math::Cephes qw(:hypergeometrics);
imports
hyp2f0: Gauss hypergeometric function F
hyp2f1: Gauss hypergeometric function F
hyperg: Confluent hypergeometric function
onef2: Hypergeometric function 1F2
threef0: Hypergeometric function 3F0
=item use Math::Cephes qw(:misc);
imports
airy: Airy function
bernum: Bernoulli numbers
dawsn: Dawson's Integral
ei: Exponential integral
erf: Error function
erfc: Complementary error function
expn: Exponential integral En
fresnl: Fresnel integral
plancki: Integral of Planck's black body radiation formula
polylog: Polylogarithm function
shichi: Hyperbolic sine and cosine integrals
sici: Sine and cosine integrals
simpson: Simpson's rule to find an integral
spence: Dilogarithm
struve: Struve function
vecang: angle between two vectors
zeta: Riemann zeta function of two arguments
zetac: Riemann zeta function
=item use Math::Cephes qw(:fract);
imports
new_fract: create a new fraction object
radd: add two fractions
rmul: multiply two fractions
rsub: subtracttwo fractions
rdiv: divide two fractions
euclid: finds the greatest common divisor
=back
=head1 FUNCTIONS
A description of the various functions available follows.
=over 4
=item I<acosh>: Inverse hyperbolic cosine
SYNOPSIS:
# double x, y, acosh();
$y = acosh( $x );
DESCRIPTION:
Returns inverse hyperbolic cosine of argument.
If 1 <= x < 1.5, a rational approximation
sqrt(z) * P(z)/Q(z)
where z = x-1, is used. Otherwise,
acosh(x) = log( x + sqrt( (x-1)(x+1) ).
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC 1,3 30000 4.2e-17 1.1e-17
IEEE 1,3 30000 4.6e-16 8.7e-17
ERROR MESSAGES:
message condition value returned
acosh domain |x| < 1 NAN
lib/Math/Cephes.pod view on Meta::CPAN
DESCRIPTION:
Returns radian angle whose tangent is y/x.
Define compile time symbol ANSIC = 1 for ANSI standard,
range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
0 to 2PI, args (x,y).
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
IEEE -10, 10 10^6 2.5e-16 6.9e-17
See atan.c.
=item I<atanh>: Inverse hyperbolic tangent
SYNOPSIS:
# double x, y, atanh();
$y = atanh( $x );
DESCRIPTION:
Returns inverse hyperbolic tangent of argument in the range
MINLOG to MAXLOG.
If |x| < 0.5, the rational form x + x**3 P(x)/Q(x) is
employed. Otherwise,
atanh(x) = 0.5 * log( (1+x)/(1-x) ).
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC -1,1 50000 2.4e-17 6.4e-18
IEEE -1,1 30000 1.9e-16 5.2e-17
=item I<bdtr>: Binomial distribution
SYNOPSIS:
# int k, n;
# double p, y, bdtr();
$y = bdtr( $k, $n, $p );
DESCRIPTION:
Returns the sum of the terms 0 through k of the Binomial
probability density:
k
-- ( n ) j n-j
> ( ) p (1-p)
-- ( j )
j=0
The terms are not summed directly; instead the incomplete
beta integral is employed, according to the formula
y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
The arguments must be positive, with p ranging from 0 to 1.
ACCURACY:
Tested at random points (a,b,p), with p between 0 and 1.
a,b Relative error:
arithmetic domain # trials peak rms
For p between 0.001 and 1:
IEEE 0,100 100000 4.3e-15 2.6e-16
See also incbet.c.
ERROR MESSAGES:
message condition value returned
bdtr domain k < 0 0.0
n < k
x < 0, x > 1
=item I<bdtrc>: Complemented binomial distribution
SYNOPSIS:
# int k, n;
# double p, y, bdtrc();
$y = bdtrc( $k, $n, $p );
DESCRIPTION:
Returns the sum of the terms k+1 through n of the Binomial
probability density:
n
-- ( n ) j n-j
> ( ) p (1-p)
-- ( j )
j=k+1
The terms are not summed directly; instead the incomplete
beta integral is employed, according to the formula
y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
The arguments must be positive, with p ranging from 0 to 1.
ACCURACY:
Tested at random points (a,b,p).
a,b Relative error:
arithmetic domain # trials peak rms
For p between 0.001 and 1:
IEEE 0,100 100000 6.7e-15 8.2e-16
For p between 0 and .001:
IEEE 0,100 100000 1.5e-13 2.7e-15
ERROR MESSAGES:
message condition value returned
bdtrc domain x<0, x>1, n<k 0.0
=item I<bdtri>: Inverse binomial distribution
SYNOPSIS:
# int k, n;
# double p, y, bdtri();
$p = bdtr( $k, $n, $y );
DESCRIPTION:
Finds the event probability p such that the sum of the
terms 0 through k of the Binomial probability density
is equal to the given cumulative probability y.
This is accomplished using the inverse beta integral
function and the relation
1 - p = incbi( n-k, k+1, y ).
ACCURACY:
Tested at random points (a,b,p).
a,b Relative error:
arithmetic domain # trials peak rms
For p between 0.001 and 1:
IEEE 0,100 100000 2.3e-14 6.4e-16
IEEE 0,10000 100000 6.6e-12 1.2e-13
For p between 10^-6 and 0.001:
IEEE 0,100 100000 2.0e-12 1.3e-14
IEEE 0,10000 100000 1.5e-12 3.2e-14
See also incbi.c.
ERROR MESSAGES:
message condition value returned
bdtri domain k < 0, n <= k 0.0
x < 0, x > 1
lib/Math/Cephes.pod view on Meta::CPAN
The complemented function is
1 - P(1-x) = incbet( b, a, x );
ACCURACY:
See incbet.c.
=item I<cbrt>: Cube root
SYNOPSIS:
# double x, y, cbrt();
$y = cbrt( $x );
DESCRIPTION:
Returns the cube root of the argument, which may be negative.
Range reduction involves determining the power of 2 of
the argument. A polynomial of degree 2 applied to the
mantissa, and multiplication by the cube root of 1, 2, or 4
approximates the root to within about 0.1%. Then Newton's
iteration is used three times to converge to an accurate
result.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC -10,10 200000 1.8e-17 6.2e-18
IEEE 0,1e308 30000 1.5e-16 5.0e-17
=item I<chdtr>: Chi-square distribution
SYNOPSIS:
# double v, x, y, chdtr();
$y = chdtr( $v, $x );
DESCRIPTION:
Returns the area under the left hand tail (from 0 to x)
of the Chi square probability density function with
v degrees of freedom.
inf.
-
1 | | v/2-1 -t/2
P( x | v ) = ----------- | t e dt
v/2 - | |
2 | (v/2) -
x
where x is the Chi-square variable.
The incomplete gamma integral is used, according to the
formula
y = chdtr( v, x ) = igam( v/2.0, x/2.0 ).
The arguments must both be positive.
ACCURACY:
See igam().
ERROR MESSAGES:
message condition value returned
chdtr domain x < 0 or v < 1 0.0
=item I<chdtrc>: Complemented Chi-square distribution
SYNOPSIS:
# double v, x, y, chdtrc();
$y = chdtrc( $v, $x );
DESCRIPTION:
Returns the area under the right hand tail (from x to
infinity) of the Chi square probability density function
with v degrees of freedom:
inf.
-
1 | | v/2-1 -t/2
P( x | v ) = ----------- | t e dt
v/2 - | |
2 | (v/2) -
x
where x is the Chi-square variable.
The incomplete gamma integral is used, according to the
formula
y = chdtrc( v, x ) = igamc( v/2.0, x/2.0 ).
The arguments must both be positive.
ACCURACY:
See igamc().
ERROR MESSAGES:
message condition value returned
chdtrc domain x < 0 or v < 1 0.0
=item I<chdtri>: Inverse of complemented Chi-square distribution
SYNOPSIS:
# double df, x, y, chdtri();
$x = chdtri( $df, $y );
DESCRIPTION:
Finds the Chi-square argument x such that the integral
from x to infinity of the Chi-square density is equal
to the given cumulative probability y.
This is accomplished using the inverse gamma integral
function and the relation
x/2 = igami( df/2, y );
ACCURACY:
See igami.c.
ERROR MESSAGES:
message condition value returned
chdtri domain y < 0 or y > 1 0.0
v < 1
=item I<clog>: Complex natural logarithm
SYNOPSIS:
# void clog();
# cmplx z, w;
$z = new_cmplx(2, 3); # $z = 2 + 3 i
$w = new_cmplx();
clog($z, $w );
print $w->{r}, ' ', $w->{i}; # prints real and imaginary parts of $w
DESCRIPTION:
Returns complex logarithm to the base e (2.718...) of
the complex argument x.
lib/Math/Cephes.pod view on Meta::CPAN
exp underflow x < -MAXL2 0.0
exp overflow x > MAXL2 MAXNUM
For DEC arithmetic, MAXL2 = 127.
For IEEE arithmetic, MAXL2 = 1024.
=item I<ei>: Exponential integral
SYNOPSIS:
#double x, y, ei();
$y = ei( $x );
DESCRIPTION:
x
- t
| | e
Ei(x) = -|- --- dt .
| | t
-
-inf
Not defined for x <= 0.
See also expn.c.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
IEEE 0,100 50000 8.6e-16 1.3e-16
=item I<expn>: Exponential integral En
SYNOPSIS:
# int n;
# double x, y, expn();
$y = expn( $n, $x );
DESCRIPTION:
Evaluates the exponential integral
inf.
-
| | -xt
| e
E (x) = | ---- dt.
n | n
| | t
-
1
Both n and x must be nonnegative.
The routine employs either a power series, a continued
fraction, or an asymptotic formula depending on the
relative values of n and x.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC 0, 30 5000 2.0e-16 4.6e-17
IEEE 0, 30 10000 1.7e-15 3.6e-16
=item I<fabs>: Absolute value
SYNOPSIS:
# double x, y;
$y = fabs( $x );
DESCRIPTION:
Returns the absolute value of the argument.
=item I<fac>: Factorial function
SYNOPSIS:
# double y, fac();
# int i;
$y = fac( $i );
DESCRIPTION:
Returns factorial of i = 1 * 2 * 3 * ... * i.
fac(0) = 1.0.
Due to machine arithmetic bounds the largest value of
i accepted is 33 in DEC arithmetic or 170 in IEEE
arithmetic. Greater values, or negative ones,
produce an error message and return MAXNUM.
ACCURACY:
For i < 34 the values are simply tabulated, and have
full machine accuracy. If i > 55, fac(i) = gamma(i+1);
see gamma.c.
Relative error:
arithmetic domain peak
IEEE 0, 170 1.4e-15
DEC 0, 33 1.4e-17
=item I<fdtr>: F distribution
SYNOPSIS:
# int df1, df2;
# double x, y, fdtr();
$y = fdtr( $df1, $df2, $x );
DESCRIPTION:
Returns the area from zero to x under the F density
function (also known as Snedcor's density or the
variance ratio density). This is the density
of x = (u1/df1)/(u2/df2), where u1 and u2 are random
variables having Chi square distributions with df1
and df2 degrees of freedom, respectively.
The incomplete beta integral is used, according to the
formula
P(x) = incbet( df1/2, df2/2, df1*x/(df2 + df1*x) ).
The arguments a and b are greater than zero, and x is
nonnegative.
ACCURACY:
Tested at random points (a,b,x).
x a,b Relative error:
arithmetic domain domain # trials peak rms
IEEE 0,1 0,100 100000 9.8e-15 1.7e-15
IEEE 1,5 0,100 100000 6.5e-15 3.5e-16
IEEE 0,1 1,10000 100000 2.2e-11 3.3e-12
IEEE 1,5 1,10000 100000 1.1e-11 1.7e-13
See also incbet.c.
ERROR MESSAGES:
message condition value returned
fdtr domain a<0, b<0, x<0 0.0
=item I<fdtrc>: Complemented F distribution
SYNOPSIS:
# int df1, df2;
# double x, y, fdtrc();
$y = fdtrc( $df1, $df2, $x );
DESCRIPTION:
Returns the area from x to infinity under the F density
function (also known as Snedcor's density or the
variance ratio density).
inf.
-
1 | | a-1 b-1
1-P(x) = ------ | t (1-t) dt
B(a,b) | |
-
x
The incomplete beta integral is used, according to the
formula
P(x) = incbet( df2/2, df1/2, df2/(df2 + df1*x) ).
ACCURACY:
Tested at random points (a,b,x) in the indicated intervals.
x a,b Relative error:
arithmetic domain domain # trials peak rms
IEEE 0,1 1,100 100000 3.7e-14 5.9e-16
IEEE 1,5 1,100 100000 8.0e-15 1.6e-15
IEEE 0,1 1,10000 100000 1.8e-11 3.5e-13
IEEE 1,5 1,10000 100000 2.0e-11 3.0e-12
See also incbet.c.
ERROR MESSAGES:
message condition value returned
fdtrc domain a<0, b<0, x<0 0.0
=item I<fdtri>: Inverse of complemented F distribution
SYNOPSIS:
# int df1, df2;
# double x, p, fdtri();
$x = fdtri( $df1, $df2, $p );
DESCRIPTION:
Finds the F density argument x such that the integral
from x to infinity of the F density is equal to the
given probability p.
This is accomplished using the inverse beta integral
function and the relations
z = incbi( df2/2, df1/2, p )
x = df2 (1-z) / (df1 z).
Note: the following relations hold for the inverse of
the uncomplemented F distribution:
z = incbi( df1/2, df2/2, p )
x = df2 z / (df1 (1-z)).
ACCURACY:
Tested at random points (a,b,p).
a,b Relative error:
arithmetic domain # trials peak rms
For p between .001 and 1:
IEEE 1,100 100000 8.3e-15 4.7e-16
IEEE 1,10000 100000 2.1e-11 1.4e-13
For p between 10^-6 and 10^-3:
IEEE 1,100 50000 1.3e-12 8.4e-15
IEEE 1,10000 50000 3.0e-12 4.8e-14
See also fdtrc.c.
lib/Math/Cephes.pod view on Meta::CPAN
($flag, $S, $C) = fresnl( $x );
DESCRIPTION:
Evaluates the Fresnel integrals
x
-
| |
C(x) = | cos(pi/2 t**2) dt,
| |
-
0
x
-
| |
S(x) = | sin(pi/2 t**2) dt.
| |
-
0
The integrals are evaluated by a power series for x < 1.
For x >= 1 auxiliary functions f(x) and g(x) are employed
such that
C(x) = 0.5 + f(x) sin( pi/2 x**2 ) - g(x) cos( pi/2 x**2 )
S(x) = 0.5 - f(x) cos( pi/2 x**2 ) - g(x) sin( pi/2 x**2 )
ACCURACY:
Relative error.
Arithmetic function domain # trials peak rms
IEEE S(x) 0, 10 10000 2.0e-15 3.2e-16
IEEE C(x) 0, 10 10000 1.8e-15 3.3e-16
DEC S(x) 0, 10 6000 2.2e-16 3.9e-17
DEC C(x) 0, 10 5000 2.3e-16 3.9e-17
=item I<gamma>: Gamma function
SYNOPSIS:
# double x, y, gamma();
# extern int sgngam;
$y = gamma( $x );
DESCRIPTION:
Returns gamma function of the argument. The result is
correctly signed, and the sign (+1 or -1) is also
returned in a global (extern) variable named sgngam.
This variable is also filled in by the logarithmic gamma
function lgam().
Arguments |x| <= 34 are reduced by recurrence and the function
approximated by a rational function of degree 6/7 in the
interval (2,3). Large arguments are handled by Stirling's
formula. Large negative arguments are made positive using
a reflection formula.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC -34, 34 10000 1.3e-16 2.5e-17
IEEE -170,-33 20000 2.3e-15 3.3e-16
IEEE -33, 33 20000 9.4e-16 2.2e-16
IEEE 33, 171.6 20000 2.3e-15 3.2e-16
Error for arguments outside the test range will be larger
owing to error amplification by the exponential function.
=item I<lgam>: Natural logarithm of gamma function
SYNOPSIS:
# double x, y, lgam();
# extern int sgngam;
$y = lgam( $x );
DESCRIPTION:
Returns the base e (2.718...) logarithm of the absolute
value of the gamma function of the argument.
The sign (+1 or -1) of the gamma function is returned in a
global (extern) variable named sgngam.
For arguments greater than 13, the logarithm of the gamma
function is approximated by the logarithmic version of
Stirling's formula using a polynomial approximation of
degree 4. Arguments between -33 and +33 are reduced by
recurrence to the interval [2,3] of a rational approximation.
The cosecant reflection formula is employed for arguments
less than -33.
Arguments greater than MAXLGM return MAXNUM and an error
message. MAXLGM = 2.035093e36 for DEC
arithmetic or 2.556348e305 for IEEE arithmetic.
ACCURACY:
arithmetic domain # trials peak rms
DEC 0, 3 7000 5.2e-17 1.3e-17
DEC 2.718, 2.035e36 5000 3.9e-17 9.9e-18
IEEE 0, 3 28000 5.4e-16 1.1e-16
IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
The error criterion was relative when the function magnitude
was greater than one but absolute when it was less than one.
The following test used the relative error criterion, though
at certain points the relative error could be much higher than
indicated.
IEEE -200, -4 10000 4.8e-16 1.3e-16
=item I<gdtr>: Gamma distribution function
SYNOPSIS:
# double a, b, x, y, gdtr();
$y = gdtr( $a, $b, $x );
DESCRIPTION:
Returns the integral from zero to x of the gamma probability
density function:
x
b -
a | | b-1 -at
y = ----- | t e dt
- | |
| (b) -
0
The incomplete gamma integral is used, according to the
relation
y = igam( b, ax ).
ACCURACY:
See igam().
ERROR MESSAGES:
message condition value returned
gdtr domain x < 0 0.0
=item I<gdtrc>: Complemented gamma distribution function
SYNOPSIS:
lib/Math/Cephes.pod view on Meta::CPAN
hyp2f1( a, b, c, x ) = F ( a, b; c; x )
2 1
inf.
- a(a+1)...(a+k) b(b+1)...(b+k) k+1
= 1 + > ----------------------------- x .
- c(c+1)...(c+k) (k+1)!
k = 0
Cases addressed are
Tests and escapes for negative integer a, b, or c
Linear transformation if c - a or c - b negative integer
Special case c = a or c = b
Linear transformation for x near +1
Transformation for x < -0.5
Psi function expansion if x > 0.5 and c - a - b integer
Conditionally, a recurrence on c to make c-a-b > 0
|x| > 1 is rejected.
The parameters a, b, c are considered to be integer
valued if they are within 1.0e-14 of the nearest integer
(1.0e-13 for IEEE arithmetic).
ACCURACY:
Relative error (-1 < x < 1):
arithmetic domain # trials peak rms
IEEE -1,7 230000 1.2e-11 5.2e-14
Several special cases also tested with a, b, c in
the range -7 to 7.
ERROR MESSAGES:
A "partial loss of precision" message is printed if
the internally estimated relative error exceeds 1^-12.
A "singularity" message is printed on overflow or
in cases not addressed (such as x < -1).
=item I<hyperg>: Confluent hypergeometric function
SYNOPSIS:
# double a, b, x, y, hyperg();
$y = hyperg( $a, $b, $x );
DESCRIPTION:
Computes the confluent hypergeometric function
1 2
a x a(a+1) x
F ( a,b;x ) = 1 + ---- + --------- + ...
1 1 b 1! b(b+1) 2!
Many higher transcendental functions are special cases of
this power series.
As is evident from the formula, b must not be a negative
integer or zero unless a is an integer with 0 >= a > b.
The routine attempts both a direct summation of the series
and an asymptotic expansion. In each case error due to
roundoff, cancellation, and nonconvergence is estimated.
The result with smaller estimated error is returned.
ACCURACY:
Tested at random points (a, b, x), all three variables
ranging from 0 to 30.
Relative error:
arithmetic domain # trials peak rms
DEC 0,30 2000 1.2e-15 1.3e-16
IEEE 0,30 30000 1.8e-14 1.1e-15
Larger errors can be observed when b is near a negative
integer or zero. Certain combinations of arguments yield
serious cancellation error in the power series summation
and also are not in the region of near convergence of the
asymptotic series. An error message is printed if the
self-estimated relative error is greater than 1.0e-12.
=item I<i0>: Modified Bessel function of order zero
SYNOPSIS:
# double x, y, i0();
$y = i0( $x );
DESCRIPTION:
Returns modified Bessel function of order zero of the
argument.
The function is defined as i0(x) = j0( ix ).
The range is partitioned into the two intervals [0,8] and
(8, infinity). Chebyshev polynomial expansions are employed
in each interval.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC 0,30 6000 8.2e-17 1.9e-17
IEEE 0,30 30000 5.8e-16 1.4e-16
=item I<i0e>: Modified Bessel function of order zero, exponentially scaled
SYNOPSIS:
# double x, y, i0e();
$y = i0e( $x );
DESCRIPTION:
Returns exponentially scaled modified Bessel function
lib/Math/Cephes.pod view on Meta::CPAN
IEEE 0,5 10000 6.9e-15 4.5e-16
IEEE 0,85 250000 2.2e-13 1.7e-14
IEEE 0,1000 30000 5.3e-12 6.3e-13
IEEE 0,10000 250000 9.3e-11 7.1e-12
IEEE 0,100000 10000 8.7e-10 4.8e-11
Outputs smaller than the IEEE gradual underflow threshold
were excluded from these statistics.
ERROR MESSAGES:
message condition value returned
incbet domain x<0, x>1 0.0
incbet underflow 0.0
=item I<incbi>: Inverse of imcomplete beta integral
SYNOPSIS:
# double a, b, x, y, incbi();
$x = incbi( $a, $b, $y );
DESCRIPTION:
Given y, the function finds x such that
incbet( a, b, x ) = y .
The routine performs interval halving or Newton iterations to find the
root of incbet(a,b,x) - y = 0.
ACCURACY:
Relative error:
x a,b
arithmetic domain domain # trials peak rms
IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13
IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15
IEEE 0,1 0,5 50000 1.1e-12 5.5e-15
VAX 0,1 .5,100 25000 3.5e-14 1.1e-15
With a and b constrained to half-integer or integer values:
IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13
IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16
With a = .5, b constrained to half-integer or integer values:
IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11
=item I<iv>: Modified Bessel function of noninteger order
SYNOPSIS:
# double v, x, y, iv();
$y = iv( $v, $x );
DESCRIPTION:
Returns modified Bessel function of order v of the
argument. If x is negative, v must be integer valued.
The function is defined as Iv(x) = Jv( ix ). It is
here computed in terms of the confluent hypergeometric
function, according to the formula
v -x
Iv(x) = (x/2) e hyperg( v+0.5, 2v+1, 2x ) / gamma(v+1)
If v is a negative integer, then v is replaced by -v.
ACCURACY:
Tested at random points (v, x), with v between 0 and
30, x between 0 and 28.
Relative error:
arithmetic domain # trials peak rms
DEC 0,30 2000 3.1e-15 5.4e-16
IEEE 0,30 10000 1.7e-14 2.7e-15
Accuracy is diminished if v is near a negative integer.
See also hyperg.c.
=item I<j0>: Bessel function of order zero
SYNOPSIS:
# double x, y, j0();
$y = j0( $x );
DESCRIPTION:
Returns Bessel function of order zero of the argument.
The domain is divided into the intervals [0, 5] and
(5, infinity). In the first interval the following rational
approximation is used:
2 2
(w - r ) (w - r ) P (w) / Q (w)
1 2 3 8
2
where w = x and the two r's are zeros of the function.
In the second interval, the Hankel asymptotic expansion
is employed with two rational functions of degree 6/6
and 7/7.
ACCURACY:
Absolute error:
arithmetic domain # trials peak rms
DEC 0, 30 10000 4.4e-17 6.3e-18
IEEE 0, 30 60000 4.2e-16 1.1e-16
=item I<y0>: Bessel function of the second kind, order zero
SYNOPSIS:
# double x, y, y0();
$y = y0( $x );
lib/Math/Cephes.pod view on Meta::CPAN
char *fctnam;
# int code;
# int mtherr();
mtherr( $fctnam, $code );
DESCRIPTION:
This routine may be called to report one of the following
error conditions (in the include file mconf.h).
Mnemonic Value Significance
DOMAIN 1 argument domain error
SING 2 function singularity
OVERFLOW 3 overflow range error
UNDERFLOW 4 underflow range error
TLOSS 5 total loss of precision
PLOSS 6 partial loss of precision
EDOM 33 Unix domain error code
ERANGE 34 Unix range error code
The default version of the file prints the function name,
passed to it by the pointer fctnam, followed by the
error condition. The display is directed to the standard
output device. The routine then returns to the calling
program. Users may wish to modify the program to abort by
calling exit() under severe error conditions such as domain
errors.
Since all error conditions pass control to this function,
the display may be easily changed, eliminated, or directed
to an error logging device.
SEE ALSO: mconf.h
=item I<nbdtr>: Negative binomial distribution
SYNOPSIS:
# int k, n;
# double p, y, nbdtr();
$y = nbdtr( $k, $n, $p );
DESCRIPTION:
Returns the sum of the terms 0 through k of the negative
binomial distribution:
k
-- ( n+j-1 ) n j
> ( ) p (1-p)
-- ( j )
j=0
In a sequence of Bernoulli trials, this is the probability
that k or fewer failures precede the nth success.
The terms are not computed individually; instead the incomplete
beta integral is employed, according to the formula
y = nbdtr( k, n, p ) = incbet( n, k+1, p ).
The arguments must be positive, with p ranging from 0 to 1.
ACCURACY:
Tested at random points (a,b,p), with p between 0 and 1.
a,b Relative error:
arithmetic domain # trials peak rms
IEEE 0,100 100000 1.7e-13 8.8e-15
See also incbet.c.
=item I<nbdtrc>: Complemented negative binomial distribution
SYNOPSIS:
# int k, n;
# double p, y, nbdtrc();
$y = nbdtrc( $k, $n, $p );
DESCRIPTION:
Returns the sum of the terms k+1 to infinity of the negative
binomial distribution:
inf
-- ( n+j-1 ) n j
> ( ) p (1-p)
-- ( j )
j=k+1
The terms are not computed individually; instead the incomplete
beta integral is employed, according to the formula
y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
The arguments must be positive, with p ranging from 0 to 1.
ACCURACY:
Tested at random points (a,b,p), with p between 0 and 1.
a,b Relative error:
arithmetic domain # trials peak rms
IEEE 0,100 100000 1.7e-13 8.8e-15
See also incbet.c.
=item I<nbdtrc>: Complemented negative binomial distribution
SYNOPSIS:
# int k, n;
# double p, y, nbdtrc();
$y = nbdtrc( $k, $n, $p );
DESCRIPTION:
Returns the sum of the terms k+1 to infinity of the negative
binomial distribution:
inf
-- ( n+j-1 ) n j
> ( ) p (1-p)
-- ( j )
j=k+1
The terms are not computed individually; instead the incomplete
beta integral is employed, according to the formula
y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
The arguments must be positive, with p ranging from 0 to 1.
ACCURACY:
See incbet.c.
=item I<nbdtri>: Functional inverse of negative binomial distribution
SYNOPSIS:
# int k, n;
# double p, y, nbdtri();
$p = nbdtri( $k, $n, $y );
DESCRIPTION:
Finds the argument p such that nbdtr(k,n,p) is equal to y.
ACCURACY:
Tested at random points (a,b,y), with y between 0 and 1.
a,b Relative error:
arithmetic domain # trials peak rms
IEEE 0,100 100000 1.5e-14 8.5e-16
See also incbi.c.
=item I<ndtr>: Normal distribution function
SYNOPSIS:
# double x, y, ndtr();
$y = ndtr( $x );
DESCRIPTION:
Returns the area under the Gaussian probability density
function, integrated from minus infinity to x:
x
-
1 | | 2
ndtr(x) = --------- | exp( - t /2 ) dt
sqrt(2pi) | |
-
-inf.
= ( 1 + erf(z) ) / 2
where z = x/sqrt(2). Computation is via the functions
erf and erfc.
ACCURACY:
Relative error:
lib/Math/Cephes.pod view on Meta::CPAN
IEEE 3e-308, 0.135 50000 4.6e-16 9.8e-17
ERROR MESSAGES:
message condition value returned
ndtri domain x <= 0 -MAXNUM
ndtri domain x >= 1 MAXNUM
=item I<pdtr>: Poisson distribution
SYNOPSIS:
# int k;
# double m, y, pdtr();
$y = pdtr( $k, $m );
DESCRIPTION:
Returns the sum of the first k terms of the Poisson
distribution:
k j
-- -m m
> e --
-- j!
j=0
The terms are not summed directly; instead the incomplete
gamma integral is employed, according to the relation
y = pdtr( k, m ) = igamc( k+1, m ).
The arguments must both be positive.
ACCURACY:
See igamc().
=item I<pdtrc>: Complemented poisson distribution
SYNOPSIS:
# int k;
# double m, y, pdtrc();
$y = pdtrc( $k, $m );
DESCRIPTION:
Returns the sum of the terms k+1 to infinity of the Poisson
distribution:
inf. j
-- -m m
> e --
-- j!
j=k+1
The terms are not summed directly; instead the incomplete
gamma integral is employed, according to the formula
y = pdtrc( k, m ) = igam( k+1, m ).
The arguments must both be positive.
ACCURACY:
See igam.c.
=item I<pdtri>: Inverse Poisson distribution
SYNOPSIS:
# int k;
# double m, y, pdtr();
$m = pdtri( $k, $y );
DESCRIPTION:
Finds the Poisson variable x such that the integral
from 0 to x of the Poisson density is equal to the
given probability y.
This is accomplished using the inverse gamma integral
function and the relation
m = igami( k+1, y ).
ACCURACY:
See igami.c.
ERROR MESSAGES:
message condition value returned
pdtri domain y < 0 or y >= 1 0.0
k < 0
=item I<pow>: Power function
SYNOPSIS:
# double x, y, z, pow();
$z = pow( $x, $y );
DESCRIPTION:
Computes x raised to the yth power. Analytically,
x**y = exp( y log(x) ).
Following Cody and Waite, this program uses a lookup table
of 2**-i/16 and pseudo extended precision arithmetic to
obtain an extra three bits of accuracy in both the logarithm
and the exponential.
ACCURACY:
lib/Math/Cephes.pod view on Meta::CPAN
-26 < y < 26, y uniformly distributed.
IEEE 0,8700 30000 1.5e-14 2.1e-15
0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
ERROR MESSAGES:
message condition value returned
pow overflow x**y > MAXNUM INFINITY
pow underflow x**y < 1/MAXNUM 0.0
pow domain x<0 and y noninteger 0.0
=item I<powi>: Real raised to integer power
SYNOPSIS:
# double x, y, powi();
# int n;
$y = powi( $x, $n );
DESCRIPTION:
Returns argument x raised to the nth power.
The routine efficiently decomposes n as a sum of powers of
two. The desired power is a product of two-to-the-kth
powers of x. Thus to compute the 32767 power of x requires
28 multiplications instead of 32767 multiplications.
ACCURACY:
Relative error:
arithmetic x domain n domain # trials peak rms
DEC .04,26 -26,26 100000 2.7e-16 4.3e-17
IEEE .04,26 -26,26 50000 2.0e-15 3.8e-16
IEEE 1,2 -1022,1023 50000 8.6e-14 1.6e-14
Returns MAXNUM on overflow, zero on underflow.
=item I<psi>: Psi (digamma) function
SYNOPSIS:
# double x, y, psi();
$y = psi( $x );
DESCRIPTION:
d -
psi(x) = -- ln | (x)
dx
is the logarithmic derivative of the gamma function.
For integer x,
n-1
-
psi(n) = -EUL + > 1/k.
-
k=1
This formula is used for 0 < n <= 10. If x is negative, it
is transformed to a positive argument by the reflection
formula psi(1-x) = psi(x) + pi cot(pi x).
For general positive x, the argument is made greater than 10
using the recurrence psi(x+1) = psi(x) + 1/x.
Then the following asymptotic expansion is applied:
inf. B
- 2k
psi(x) = log(x) - 1/2x - > -------
- 2k
k=1 2k x
where the B2k are Bernoulli numbers.
ACCURACY:
Relative error (except absolute when |psi| < 1):
arithmetic domain # trials peak rms
DEC 0,30 2500 1.7e-16 2.0e-17
IEEE 0,30 30000 1.3e-15 1.4e-16
IEEE -30,0 40000 1.5e-15 2.2e-16
ERROR MESSAGES:
message condition value returned
psi singularity x integer <=0 MAXNUM
=item I<rgamma>: Reciprocal gamma function
SYNOPSIS:
# double x, y, rgamma();
$y = rgamma( $x );
DESCRIPTION:
Returns one divided by the gamma function of the argument.
The function is approximated by a Chebyshev expansion in
the interval [0,1]. Range reduction is by recurrence
for arguments between -34.034 and +34.84425627277176174.
1/MAXNUM is returned for positive arguments outside this
range. For arguments less than -34.034 the cosecant
reflection formula is applied; lograrithms are employed
to avoid unnecessary overflow.
The reciprocal gamma function has no singularities,
but overflow and underflow may occur for large arguments.
These conditions return either MAXNUM or 1/MAXNUM with
appropriate sign.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC -30,+30 4000 1.2e-16 1.8e-17
IEEE -30,+30 30000 1.1e-15 2.0e-16
For arguments less than -34.034 the peak error is on the
order of 5e-15 (DEC), excepting overflow or underflow.
=item I<round>: Round double to nearest or even integer valued double
SYNOPSIS:
# double x, y, round();
$y = round( $x );
DESCRIPTION:
Returns the nearest integer to x as a double precision
floating point result. If x ends in 0.5 exactly, the
nearest even integer is chosen.
ACCURACY:
If x is greater than 1/(2*MACHEP), its closest machine
representation is already an integer, so rounding does
not change it.
=item I<shichi>: Hyperbolic sine and cosine integrals
SYNOPSIS:
# double x, Chi, Shi, shichi();
($flag, $Shi, $Chi) = shichi( $x );
DESCRIPTION:
Approximates the integrals
x
-
| | cosh t - 1
Chi(x) = eul + ln x + | ----------- dt,
| | t
-
0
x
-
| | sinh t
Shi(x) = | ------ dt
lib/Math/Cephes.pod view on Meta::CPAN
Two polynomial approximating functions are employed.
Between 0 and pi/4 the cosine is approximated by
1 - x**2 P(x**2).
Between pi/4 and pi/2 the sine is represented as
x + x**3 P(x**2).
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC +-1000 3400 3.5e-17 9.1e-18
IEEE +-1000 30000 2.1e-16 5.7e-17
See also sin().
=item I<sinh>: Hyperbolic sine
SYNOPSIS:
# double x, y, sinh();
$y = sinh( $x );
DESCRIPTION:
Returns hyperbolic sine of argument in the range MINLOG to
MAXLOG.
The range is partitioned into two segments. If |x| <= 1, a
rational function of the form x + x**3 P(x)/Q(x) is employed.
Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC +- 88 50000 4.0e-17 7.7e-18
IEEE +-MAXLOG 30000 2.6e-16 5.7e-17
=item I<spence>: Dilogarithm
SYNOPSIS:
# double x, y, spence();
$y = spence( $x );
DESCRIPTION:
Computes the integral
x
-
| | log t
spence(x) = - | ----- dt
| | t - 1
-
1
for x >= 0. A rational approximation gives the integral in
the interval (0.5, 1.5). Transformation formulas for 1/x
and 1-x are employed outside the basic expansion range.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
IEEE 0,4 30000 3.9e-15 5.4e-16
DEC 0,4 3000 2.5e-16 4.5e-17
=item I<sqrt>: Square root
SYNOPSIS:
# double x, y, sqrt();
$y = sqrt( $x );
DESCRIPTION:
Returns the square root of x.
Range reduction involves isolating the power of two of the
argument and using a polynomial approximation to obtain
a rough value for the square root. Then Heron's iteration
is used three times to converge to an accurate value.
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
DEC 0, 10 60000 2.1e-17 7.9e-18
IEEE 0,1.7e308 30000 1.7e-16 6.3e-17
ERROR MESSAGES:
message condition value returned
sqrt domain x < 0 0.0
=item I<stdtr>: Student's t distribution
SYNOPSIS:
# double t, stdtr();
short k;
$y = stdtr( $k, $t );
DESCRIPTION:
Computes the integral from minus infinity to t of the Student
t distribution with integer k > 0 degrees of freedom:
t
-
| |
- | 2 -(k+1)/2
| ( (k+1)/2 ) | ( x )
---------------------- | ( 1 + --- ) dx
- | ( k )
sqrt( k pi ) | ( k/2 ) |
lib/Math/Cephes.pod view on Meta::CPAN
Relation to incomplete beta integral:
1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
where
z = k/(k + t**2).
For t < -2, this is the method of computation. For higher t,
a direct method is derived from integration by parts.
Since the function is symmetric about t=0, the area under the
right tail of the density is found by calling the function
with -t instead of t.
ACCURACY:
Tested at random 1 <= k <= 25. The "domain" refers to t.
Relative error:
arithmetic domain # trials peak rms
IEEE -100,-2 50000 5.9e-15 1.4e-15
IEEE -2,100 500000 2.7e-15 4.9e-17
=item I<stdtri>: Functional inverse of Student's t distribution
SYNOPSIS:
# double p, t, stdtri();
# int k;
$t = stdtri( $k, $p );
DESCRIPTION:
Given probability p, finds the argument t such that stdtr(k,t)
is equal to p.
ACCURACY:
Tested at random 1 <= k <= 100. The "domain" refers to p:
Relative error:
arithmetic domain # trials peak rms
IEEE .001,.999 25000 5.7e-15 8.0e-16
IEEE 10^-6,.001 25000 2.0e-12 2.9e-14
=item I<struve>: Struve function
SYNOPSIS:
# double v, x, y, struve();
$y = struve( $v, $x );
DESCRIPTION:
Computes the Struve function Hv(x) of order v, argument x.
Negative x is rejected unless v is an integer.
ACCURACY:
Not accurately characterized, but spot checked against tables.
=item I<plancki>: Integral of Planck's black body radiation formula
SYNOPSIS:
# double lambda, T, y, plancki()
$y = plancki( $lambda, $T );
DESCRIPTION:
Evaluates the definite integral, from wavelength 0 to lambda,
of Planck's radiation formula
-5
c1 lambda
E = ------------------
c2/(lambda T)
e - 1
Physical constants c1 = 3.7417749e-16 and c2 = 0.01438769 are built in
to the function program. They are scaled to provide a result
in watts per square meter. Argument T represents temperature in degrees
Kelvin; lambda is wavelength in meters.
The integral is expressed in closed form, in terms of polylogarithms
(see polylog.c).
The total area under the curve is
(-1/8) (42 zeta(4) - 12 pi^2 zeta(2) + pi^4 ) c1 (T/c2)^4
= (pi^4 / 15) c1 (T/c2)^4
= 5.6705032e-8 T^4
where sigma = 5.6705032e-8 W m^2 K^-4 is the Stefan-Boltzmann constant.
ACCURACY:
The left tail of the function experiences some relative error
amplification in computing the dominant term exp(-c2/(lambda T)).
For the right-hand tail see planckc, below.
Relative error.
The domain refers to lambda T / c2.
arithmetic domain # trials peak rms
IEEE 0.1, 10 50000 7.1e-15 5.4e-16
=item I<polylog>: polylogarithm function
SYNOPSIS:
# double x, y, polylog();
# int n;
$y = polylog( $n, $x );
The polylogarithm of order n is defined by the series
inf k
- x
Li (x) = > --- .
n - n
k=1 k
For x = 1,
inf
- 1
Li (1) = > --- = Riemann zeta function (n) .
n - n
k=1 k
When n = 2, the function is the dilogarithm, related to Spence's integral:
x 1-x
- -
| | -ln(1-t) | | ln t
Li (x) = | -------- dt = | ------ dt = spence(1-x) .
2 | | t | | 1 - t
- -
0 1
ACCURACY:
Relative error:
arithmetic domain n # trials peak rms
IEEE 0, 1 2 50000 6.2e-16 8.0e-17
IEEE 0, 1 3 100000 2.5e-16 6.6e-17
IEEE 0, 1 4 30000 1.7e-16 4.9e-17
IEEE 0, 1 5 30000 5.1e-16 7.8e-17
=item I<bernum>: Bernoulli numbers
SYNOPSIS:
($num, $den) = bernum( $n);
($num_array, $den_array) = bernum();
DESCRIPTION:
This calculates the Bernoulli numbers, up to 30th order.
If called with an integer argument, the numerator and denominator
of that Bernoulli number is returned; if called with no argument,
two array references representing the numerator and denominators
of the first 30 Bernoulli numbers are returned.
=item I<simpson>: Simpson's rule to find an integral
SYNOPSIS:
$result = simpson(\&fun, $a, $b, $abs_err, $rel_err, $nmax);
sub fun {
my $x = shift;
return cos($x)*exp($x);
}
DESCRIPTION:
This evaluates the area under the graph of a function,
represented in a subroutine, from $a to $b, using an 8-point
Newton-Cotes formula. The routine divides up the interval into
equal segments, evaluates the integral, then compares that
to the result with double the number of segments. If the two
results agree, to within an absolute error $abs_err or a
relative error $rel_err, the result is returned; otherwise,
the number of segments is doubled again, and the results
compared. This continues until the desired accuracy is attained,
or until the maximum number of iterations $nmax is reached.
=item I<vecang>: angle between two vectors
SYNOPSIS:
# double p[3], q[3], vecang();
$y = vecang( $p, $q );
DESCRIPTION:
For two vectors p, q, the angle A between them is given by
p.q / (|p| |q|) = cos A .
where "." represents inner product, "|x|" the length of vector x.
If the angle is small, an expression in sin A is preferred.
Set r = q - p. Then
p.q = p.p + p.r ,
|p|^2 = p.p ,
|q|^2 = p.p + 2 p.r + r.r ,
p.p^2 + 2 p.p p.r + p.r^2
cos^2 A = ----------------------------
p.p (p.p + 2 p.r + r.r)
p.p + 2 p.r + p.r^2 / p.p
= --------------------------- ,
p.p + 2 p.r + r.r
sin^2 A = 1 - cos^2 A
r.r - p.r^2 / p.p
= --------------------
p.p + 2 p.r + r.r
= (r.r - p.r^2 / p.p) / q.q .
ACCURACY:
Relative error:
arithmetic domain # trials peak rms
IEEE -1, 1 10^6 1.7e-16 4.2e-17
=item I<onef2>: Hypergeometric function 1F2
SYNOPSIS:
# double a, b, c, x, value;
lib/Math/Cephes.pod view on Meta::CPAN
$y = expm1( $x );
# cosm1(x) = cos(x) - 1
$y = cosm1( $x );
=item I<yn>: Bessel function of second kind of integer order
SYNOPSIS:
# double x, y, yn();
# int n;
$y = yn( $n, $x );
DESCRIPTION:
Returns Bessel function of order n, where n is a
(possibly negative) integer.
The function is evaluated by forward recurrence on
n, starting with values computed by the routines
y0() and y1().
If n = 0 or 1 the routine for y0 or y1 is called
directly.
ACCURACY:
Absolute error, except relative
when y > 1:
arithmetic domain # trials peak rms
DEC 0, 30 2200 2.9e-16 5.3e-17
IEEE 0, 30 30000 3.4e-15 4.3e-16
ERROR MESSAGES:
message condition value returned
yn singularity x = 0 MAXNUM
yn overflow MAXNUM
Spot checked against tables for x, n between 0 and 100.
=item I<zeta>: Riemann zeta function of two arguments
SYNOPSIS:
# double x, q, y, zeta();
$y = zeta( $x, $q );
DESCRIPTION:
inf.
- -x
zeta(x,q) = > (k+q)
-
k=0
where x > 1 and q is not a negative integer or zero.
The Euler-Maclaurin summation formula is used to obtain
the expansion
n
- -x
zeta(x,q) = > (k+q)
-
k=1
1-x inf. B x(x+1)...(x+2j)
(n+q) 1 - 2j
+ --------- - ------- + > --------------------
x-1 x - x+2j+1
2(n+q) j=1 (2j)! (n+q)
where the B2j are Bernoulli numbers. Note that (see zetac.c)
zeta(x,1) = zetac(x) + 1.
ACCURACY:
REFERENCE:
Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals,
Series, and Products, p. 1073; Academic Press, 1980.
=item I<zetac>: Riemann zeta function
SYNOPSIS:
# double x, y, zetac();
$y = zetac( $x );
DESCRIPTION:
inf.
- -x
zetac(x) = > k , x > 1,
-
k=2
is related to the Riemann zeta function by
Riemann zeta(x) = zetac(x) + 1.
Extension of the function definition for x < 1 is implemented.
Zero is returned for x > log2(MAXNUM).
An overflow error may occur for large negative x, due to the
gamma function in the reflection formula.
ACCURACY:
Tabulated values have full machine accuracy.
Relative error:
arithmetic domain # trials peak rms
IEEE 1,50 10000 9.8e-16 1.3e-16
DEC 1,50 2000 1.1e-16 1.9e-17
=back
=head1 TODO
=over 4
=item * Include more operating systems when generating mconf.h.
=back
=head1 MAINTAINER
Shlomi Fish, L<http://www.shlomifish.org/>,
L<https://metacpan.org/author/SHLOMIF> .
=head1 BUGS
Please report any on the rt.cpan.org interface:
L<https://rt.cpan.org/Dist/Display.html?Queue=Math-Cephes>
=head1 VERSION CONTROL
This distribution is maintained in this GitHub repository:
L<https://github.com/shlomif/Math-Cephes>.
=head1 SEE ALSO
For interfaces to programs which can do symbolic manipulation,
see L<PDL>, L<Math::Pari>, and L<Math::ematica>.
For a command line interface to the routines of I<Math::Cephes>,
see the included C<pmath> script. For a
different interface to the fraction and complex number routines,
see L<Math::Cephes::Fraction> and L<Math::Cephes::Complex>.
For an interface to some polynomial routines, see
L<Math::Cephes::Polynomial>, and for some matrix routines,
see L<Math::Cephes::Matrix>.
=head1 COPYRIGHT
The C code for the Cephes Math Library is
Copyright 1984, 1987, 1989, 2002 by Stephen L. Moshier,
and is available at L<http://www.netlib.org/cephes/>.
Direct inquiries to 30 Frost Street, Cambridge, MA 02140.
The file arrays.c included here to handle passing arrays
into and out of C routines comes from the PGPLOT module
of Karl Glazebrook <kgb@zzoepp.aao.gov.au>.
The perl interface is copyright 2000, 2002 by Randy Kobes.
( run in 1.820 second using v1.01-cache-2.11-cpan-39bf76dae61 )