Math-Cephes
view release on metacpan or search on metacpan
lib/Math/Cephes.pod view on Meta::CPAN
# double a, x, y, igamc();
$y = igamc( $a, $x );
DESCRIPTION:
The function is defined by
igamc(a,x) = 1 - igam(a,x)
inf.
-
1 | | -t a-1
= ----- | e t dt.
- | |
| (a) -
x
In this implementation both arguments must be positive.
The integral is evaluated by either a power series or
continued fraction expansion, depending on the relative
values of a and x.
ACCURACY:
Tested at random a, x.
a x Relative error:
arithmetic domain domain # trials peak rms
IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
=item I<igami>: Inverse of complemented imcomplete gamma integral
SYNOPSIS:
# double a, x, p, igami();
$x = igami( $a, $p );
DESCRIPTION:
Given p, the function finds x such that
igamc( a, x ) = p.
It is valid in the right-hand tail of the distribution, p < 0.5.
Starting with the approximate value
3
x = a t
where
t = 1 - d - ndtri(p) sqrt(d)
and
d = 1/9a,
the routine performs up to 10 Newton iterations to find the
root of igamc(a,x) - p = 0.
ACCURACY:
Tested at random a, p in the intervals indicated.
a p Relative error:
arithmetic domain domain # trials peak rms
IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15
IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15
IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14
=item I<incbet>: Incomplete beta integral
SYNOPSIS:
# double a, b, x, y, incbet();
$y = incbet( $a, $b, $x );
DESCRIPTION:
Returns incomplete beta integral of the arguments, evaluated
from zero to x. The function is defined as
x
- -
| (a+b) | | a-1 b-1
----------- | t (1-t) dt.
- - | |
| (a) | (b) -
0
The domain of definition is 0 <= x <= 1. In this
implementation a and b are restricted to positive values.
The integral from x to 1 may be obtained by the symmetry
relation
1 - incbet( a, b, x ) = incbet( b, a, 1-x ).
The integral is evaluated by a continued fraction expansion
or, when b*x is small, by a power series.
ACCURACY:
Tested at uniformly distributed random points (a,b,x) with a and b
in "domain" and x between 0 and 1.
Relative error
arithmetic domain # trials peak rms
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 );
lib/Math/Cephes.pod view on Meta::CPAN
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;
# double *err;
($value, $err) = onef2( $a, $b, $c, $x)
ACCURACY:
( run in 0.540 second using v1.01-cache-2.11-cpan-71847e10f99 )