Math-Cephes
view release on metacpan or search on metacpan
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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:
lib/Math/Cephes.pod view on Meta::CPAN
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:
lib/Math/Cephes.pod view on Meta::CPAN
-
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:
lib/Math/Cephes.pod view on Meta::CPAN
-
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:
lib/Math/Cephes.pod view on Meta::CPAN
| 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
lib/Math/Cephes.pod view on Meta::CPAN
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).
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
$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.
lib/Math/Cephes.pod view on Meta::CPAN
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:
lib/Math/Cephes.pod view on Meta::CPAN
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:
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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:
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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.
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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
lib/Math/Cephes.pod view on Meta::CPAN
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
libmd/bdtr.c view on Meta::CPAN
* 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:
libmd/bdtr.c view on Meta::CPAN
* 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:
libmd/chdtr.c view on Meta::CPAN
* -
* 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 md_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:
*
libmd/chdtr.c view on Meta::CPAN
* -
* 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 md_gamma integral is used, according to the
* formula
*
* y = chdtr( v, x ) = igamc( v/2.0, x/2.0 ).
*
*
* The arguments must both be positive.
*
*
*
* ACCURACY:
*
libmd/expn.c view on Meta::CPAN
* 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
*
*/
libmd/fdtr.c view on Meta::CPAN
* 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).
libmd/fdtr.c view on Meta::CPAN
* 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
libmd/gamma.c view on Meta::CPAN
*
* Returns md_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 md_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
libmd/gamma.c view on Meta::CPAN
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of the absolute
* value of the md_gamma function of the argument.
* The sign (+1 or -1) of the md_gamma function is returned in a
* global (extern) variable named sgngam.
*
* For arguments greater than 13, the logarithm of the md_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:
*
libmd/gamma.c view on Meta::CPAN
0x3fb2,0x4944,0xc9cd,0x3c51,
0x3ff0,0x0000,0x0000,0x0000
};
#define MAXGAM 171.624376956302725
static unsigned short LPI[4] = {
0x3ff2,0x50d0,0x48e7,0xa1bd,
};
#define LOGPI *(double *)LPI
#endif
/* Stirling's formula for the md_gamma function */
#if UNK
static double STIR[5] = {
7.87311395793093628397E-4,
-2.29549961613378126380E-4,
-2.68132617805781232825E-3,
3.47222221605458667310E-3,
8.33333333333482257126E-2,
};
#define MAXSTIR 143.01608
static double SQTPI = 2.50662827463100050242E0;
libmd/gamma.c view on Meta::CPAN
static double stirf();
double lgam();
#endif
#ifdef INFINITIES
extern double INFINITY;
#endif
#ifdef NANS
extern double NAN;
#endif
/* Gamma function computed by Stirling's formula.
* The polynomial STIR is valid for 33 <= x <= 172.
*/
static double stirf(x)
double x;
{
double y, w, v;
w = 1.0/x;
w = 1.0 + w * polevl( w, STIR, 4 );
y = md_exp(x);
libmd/gamma.c view on Meta::CPAN
mtherr( "md_gamma", SING );
return( MAXNUM );
#endif
}
else
return( z/((1.0 + 0.5772156649015329 * x) * x) );
}
/* A[]: Stirling's formula expansion of md_log md_gamma
* B[], C[]: md_log md_gamma function between 2 and 3
*/
#ifdef UNK
static double A[] = {
8.11614167470508450300E-4,
-5.95061904284301438324E-4,
7.93650340457716943945E-4,
-2.77777777730099687205E-3,
8.33333333333331927722E-2
};
libmd/hyperg.c view on Meta::CPAN
* 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:
libmd/hyperg.c view on Meta::CPAN
blowup:
*err = pcanc;
return( sum );
}
/* hy1f1a() */
/* asymptotic formula for hypergeometric function:
*
* ( -a
* -- ( |z|
* | (b) ( -------- 2f0( a, 1+a-b, -1/x )
* ( --
* ( | (b-a)
*
*
* x a-b )
* e |x| )
libmd/hyperg.c view on Meta::CPAN
{
temp = md_gamma(b);
asum *= temp;
acanc *= md_fabs(temp);
}
if( asum != 0.0 )
acanc /= md_fabs(asum);
acanc *= 30.0; /* fudge factor, since error of asymptotic formula
* often seems this much larger than advertised */
adone:
*err = acanc;
return( asum );
}
/* hyp2f0() */
*
*
*
* 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 ) / md_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
libmd/nbdtr.c view on Meta::CPAN
* 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:
libmd/nbdtr.c view on Meta::CPAN
* 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:
libmd/pdtr.c view on Meta::CPAN
* 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
* md_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.
libmd/planck.c view on Meta::CPAN
/* planck.c
*
* 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.
libmd/planck.c view on Meta::CPAN
*
* double lambda, T, y, planckc();
*
* y = planckc( lambda, T );
*
*
*
* DESCRIPTION:
*
* Integral from w to infinity (area under right hand tail)
* of Planck's radiation formula.
*
* The program for large lambda uses an asymptotic series in inverse
* powers of the wavelength.
*
* ACCURACY:
*
* Relative error.
* The domain refers to lambda T / c2.
* arithmetic domain # trials peak rms
* IEEE 0.6, 10 50000 1.1e-15 2.2e-16
libmd/planck.c view on Meta::CPAN
y = ((((y - 1./13305600.)*p + 1./272160.)*p - 1./5040.)*p + 1./60.)*p;
y = y - 0.125*u + 1./3.;
#endif
y = y * planck_c1 * b / (w*w*w);
return y;
}
/* planckd
*
* Planck's black body radiation formula
*
*
*
* SYNOPSIS:
*
* double lambda, T, y, planckd();
*
* y = planckd( lambda, T );
*
*
*
* DESCRIPTION:
*
* Evaluates Planck's radiation formula
* -5
* c1 lambda
* E = ------------------
* c2/(lambda T)
* e - 1
*
*/
double
planckd(w, T)
libmd/polylog.c view on Meta::CPAN
extern double MACHEP;
double
polylog (n, x)
int n;
double x;
{
double h, k, p, s, t, u, xc, z;
int i, j;
/* This recurrence provides formulas for n < 2.
d 1
-- Li (x) = --- Li (x) .
dx n x n-1
*/
if (n == -1)
{
p = 1.0 - x;
libmd/polylog.c view on Meta::CPAN
n n
*/
if (x == -1.0 && n > 1)
{
/* Li_n(1) = zeta(n) */
s = zetac ((double) n) + 1.0;
s = s * (md_powi (2.0, 1 - n) - 1.0);
return s;
}
/* Inversion formula:
* [n/2] n-2r
* n 1 n - md_log (z)
* Li (-z) + (-1) Li (-1/z) = - --- md_log (z) + 2 > ----------- Li (-1)
* n n n! - (n - 2r)! 2r
* r=1
*/
if (x < -1.0 && n > 1)
{
double q, w;
int r;
libmd/psi.c view on Meta::CPAN
* dx
*
* is the logarithmic derivative of the md_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) = md_log(x) - 1/2x - > -------
* - 2k
* k=1 2k x
*
libmd/rgamma.c view on Meta::CPAN
*
* DESCRIPTION:
*
* Returns one divided by the md_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 md_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:
libmd/simpsn.c view on Meta::CPAN
/* simpsn.c */
/* simpsn.c
* Numerical integration of function tabulated
* at equally spaced arguments
*/
/* Coefficients for Cote integration formulas */
/* Note: these numbers were computed using 40-decimal precision. */
#define NCOTE 8
/* 6th order formula */
/*
static double simcon[] =
{
4.88095238095238095E-2,
2.57142857142857142857E-1,
3.2142857142857142857E-2,
3.2380952380952380952E-1,
};
*/
/* 8th order formula */
static double simcon[] =
{
3.488536155202821869E-2,
2.076895943562610229E-1,
-3.27336860670194003527E-2,
3.7022927689594356261E-1,
-1.6014109347442680776E-1,
};
/* 10th order formula */
/*
static double simcon[] =
{
2.68341483619261397039E-2,
1.77535941424830313719E-1,
-8.1043570626903960237E-2,
4.5494628827962161295E-1,
-4.3515512265512265512E-1,
7.1376463043129709796E-1,
};
*/
/* simpsn.c 2 */
/* 20th order formula */
/*
static double simcon[] =
{
1.182527324903160319E-2,
1.14137717644606974987E-1,
-2.36478370511426964E-1,
1.20618689348187566E+0,
-3.7710317267153304677E+0,
1.03367982199398011435E+1,
-2.270881584397951229796E+1,
libmd/spence.c view on Meta::CPAN
*
* x
* -
* | | md_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
libmd/zeta.c view on Meta::CPAN
*
*
*
* 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
libmd/zeta.c view on Meta::CPAN
#ifdef ANSIPROT
extern double md_fabs ( double );
extern double md_pow ( double, double );
extern double md_floor ( double );
#else
double md_fabs(), md_pow(), md_floor();
#endif
extern double MAXNUM, MACHEP;
/* Expansion coefficients
* for Euler-Maclaurin summation formula
* (2k)! / B2k
* where B2k are Bernoulli numbers
*/
static double A[] = {
12.0,
-720.0,
30240.0,
-1209600.0,
47900160.0,
-1.8924375803183791606e9, /*1.307674368e12/691*/
libmd/zeta.c view on Meta::CPAN
if(q == md_floor(q))
{
mtherr( "zeta", SING );
retinf:
return( MAXNUM );
}
if( x != md_floor(x) )
goto domerr; /* because q^-x not defined */
}
/* Euler-Maclaurin summation formula */
/*
if( x < 25.0 )
*/
{
/* Permit negative q but continue sum until n+q > +9 .
* This case should be handled by a reflection formula.
* If q<0 and x is an integer, there is a relation to
* the polygamma function.
*/
s = md_pow( q, -x );
a = q;
i = 0;
b = 0.0;
while( (i < 9) || (a <= 9.0) )
{
i += 1;
libmd/zetac.c view on Meta::CPAN
* 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 > md_log2(MAXNUM).
*
* An overflow error may occur for large negative x, due to the
* md_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
*
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.
',
'nbdtri' => 'nbdtri: Functional inverse of negative binomial distribution
SYNOPSIS:
-
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.
',
'zetac' => 'zetac: Riemann zeta function
SYNOPSIS:
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.
',
'ellpj' => 'ellpj: Jacobian Elliptic Functions
SYNOPSIS:
# double u, m, sn, cn, dn, phi;
# int ellpj();
($flag, $sn, $cn, $dn, $phi) = ellpj( $u, $m );
-
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.
',
'beta' => 'beta: Beta function
SYNOPSIS:
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.
',
'chdtri' => 'chdtri: Inverse of complemented Chi-square distribution
SYNOPSIS:
# double df, x, y, chdtri();
$x = chdtri( $df, $y );
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
# 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.
',
'plancki' => 'plancki: Integral of Planck 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 the Planck 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.
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.
',
'fabs' => 'fabs: Absolute value
SYNOPSIS:
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.
',
'rgamma' => 'rgamma: Reciprocal gamma function
SYNOPSIS:
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.
',
'shichi' => 'shichi: Hyperbolic sine and cosine integrals
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.
',
'log2' => 'log2: Base 2 logarithm
| 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.
',
'dawsn' => 'dawsn: Dawson\'s Integral
SYNOPSIS:
# double x, y, dawsn();
$y = dawsn( $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
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.
',
'incbi' => 'incbi: Inverse of imcomplete beta integral
SYNOPSIS:
# double a, b, x, y, incbi();
$x = incbi( $a, $b, $y );
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.
',
'gdtr' => 'gdtr: Gamma distribution function
SYNOPSIS:
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) ).
',
'bdtri' => 'bdtri: Inverse binomial distribution
SYNOPSIS:
# int k, n;
# double p, y, bdtri();
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.
',
'igam' => 'igam: Incomplete gamma integral
SYNOPSIS:
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.
',
'cosh' => 'cosh: Hyperbolic cosine
SYNOPSIS:
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.
',
'iv' => 'iv: Modified Bessel function of noninteger order
SYNOPSIS:
$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.
'
);
}
( run in 0.787 second using v1.01-cache-2.11-cpan-26ccb49234f )