/*
Copyright (C) 2003, 2004 by J. David Sexton.
Some rights NOT reserved!
This source code should compile with any C compiler that meets the
ISO 9899:1990 ("C89") standard.
This source code is not public domain or Open Source, but I hereby license
you to copy it and give it to others so long as you don't profit by doing
so. In addition, you may include these functions in your own software, so
long as that software is not allowed to be sold or used for profit.
If I sold this, I would guarantee it; but as I don't, I won't.
*/
#include
#include
#include
long double HexTimesX (char *theString, long double theMultiplier);
long double LDFract (long double theX);
long double ATimesLog2OfX (long double theA, long double theX,
long double *theInt);
long double Log2OfX (long double theX, long double *theInt);
long double TwoToTheX (long double theX);
long double Log2OfGamma (long double theX, long double *theInt);
long double IncompleteGammaFactor (long double theA, long double theX);
double IncompleteGammaBySeries (double theA, double theX);
double IncompleteGammaByFraction (double theA, double theX);
double IncompleteGamma (double theA, double theX);
double ChiSquaredPValue (double degOfFreedom, double theChiSquared);
/*
Here begins a small library of math functions. These functions
don't call any functions in , so we have to make sure
the following macros are appropriately defined. HUGE_VAL is part
of the C89 standard. HUGE_VALL and NAN are not part of the C89
standard, but they are in C99. These macros, and other macro
definitions below, depend on definitions in . The
functions also require and .
*/
#ifndef HUGE_VAL
#define HUGE_VAL (DBL_MAX / DBL_MIN)
#endif
#ifndef HUGE_VALL
#define HUGE_VALL (LDBL_MAX / LDBL_MIN)
#endif
#ifndef NAN
#define NAN (0.0L / 0.0L)
#endif
/*
HexTimesX multiplies the hexadecimal floating-point
number specified in theString by theMultiplier.
This is a rather computationally expensive way to
insure maximum accuracy when multiplying by floating-
point literals. It's probably only worth the trouble
for literals that can't be precisely represented by
the long double type.
*/
long double
HexTimesX (char *theString, long double theMultiplier)
{
char *strPtr, theChar;
long double theProduct, placeValue;
size_t numOfDigits, theDigit;
strPtr = theString;
placeValue = theMultiplier;
if (*strPtr == '-')
{
++strPtr;
placeValue = -placeValue;
}
if (placeValue == HUGE_VALL || placeValue == -HUGE_VALL)
return (errno = ERANGE, placeValue);
theProduct = 0.0L;
numOfDigits = theDigit = 0;
while ((theChar = *(strPtr++)))
{
if (theChar == '.')
break;
if ((theChar >= 'a' && theChar <= 'f') ||
(theChar >= 'A' && theChar <= 'F') ||
(theChar >= '0' && theChar <= '9'))
++numOfDigits;
else
break;
}
if (theChar == '.')
{
while ((theChar = *(strPtr++)))
{
if ((theChar >= 'a' && theChar <= 'f') ||
(theChar >= 'A' && theChar <= 'F') ||
(theChar >= '0' && theChar <= '9'))
{
++numOfDigits;
placeValue /= 16.0L;
}
else
break;
}
}
--strPtr;
while (theDigit < numOfDigits)
{
theChar = *(--strPtr);
if (theChar == '.')
continue;
if (theChar >= 'a' && theChar <= 'f')
theChar -= 'a' - 10;
else if (theChar >= 'A' && theChar <= 'F')
theChar -= 'A' - 10;
else
theChar -= '0';
theProduct += placeValue * ((long double) theChar);
placeValue *= 16.0L;
++theDigit;
}
if (theProduct == HUGE_VALL || theProduct == -HUGE_VALL)
errno = ERANGE;
return (theProduct);
}
/*
LDFract returns theX minus the largest
integer not greater than theX.
E.g., LDFract (2.1L) returns 0.1L,
and LDFract (-2.1L) returns 0.9L.
*/
long double
LDFract (long double theX)
{
long double ldX, theIndex;
if (theX == HUGE_VALL || theX == -HUGE_VALL)
return (errno = ERANGE, 0.0L);
ldX = theX;
while (ldX < 0.0L)
{
theIndex = 2.0L;
do
{
ldX += theIndex;
theIndex *= theIndex;
}
while (-ldX >= theIndex);
}
while (ldX >= 2.0L)
{
theIndex = 2.0L;
do
{
ldX -= theIndex;
theIndex *= theIndex;
}
while (ldX >= theIndex);
}
if (ldX >= 1.0L)
--ldX;
return (ldX);
}
/*
ATimesLog2OfX calculates theA times the base 2 logarithm of
theX and returns the fractional part of this product. The
integer part of the product is returned in the variable
pointed to by theInt; this will be the largest integer not
greater than the product. E.g., if the product is 2.1,
ATimesLog2OfX will return 0.1 and 2.0; if the product is -2.1,
it will return 0.9 and -3.0.
*/
long double
ATimesLog2OfX (long double theA, long double theX, long double *theInt)
{
long double fractA, fractB, intA, intB;
fractA = Log2OfX (theX, &intA);
fractA *= theA;
intA *= theA;
fractB = LDFract (fractA);
intB = fractA - fractB;
fractA = LDFract (intA);
intA -= fractA;
if (fractA + fractB >= 1.0L)
{
--fractB;
++intB;
}
fractA += fractB;
intA += intB;
*theInt = intA;
return (fractA);
}
/*
Log2OfX returns the mantissa of the base 2 logarithm of theX.
The integer part of the logarithm is returned in the variable
pointed to by theInt; this will be the largest integer not
greater than the logarithm. E.g., if the logarithm is 2.1,
Log2OfX will return 0.1 and 2.0; if the logarithm is -2.1,
it will return 0.9 and -3.0.
If FLT_RADIX is 2, the mantissa returned should be accurate to
LDBL_MANT_DIG bits. The value returned in the variable pointed
to by theInt should be precise.
*/
long double
Log2OfX (long double theX, long double *theInt)
{
long double theSum, lastSum, pSum, qSum;
long double ldX, ldY, ldZ, theIndex;
if (theX < 0.0L)
{
errno = EDOM;
*theInt = NAN;
return (0.0L);
}
if (theX == 0.0L)
{
errno = ERANGE;
*theInt = -HUGE_VALL;
return (0.0L);
}
if (theX == HUGE_VALL)
{
errno = ERANGE;
*theInt = HUGE_VALL;
return (0.0L);
}
ldX = theX;
theSum = pSum = 0.0L;
ldZ = 1.0L;
while (ldX / ldZ < 1.0L)
{
theIndex = 1.0L;
ldY = 2.0L;
do
{
theSum -= theIndex;
ldZ /= ldY;
theIndex += theIndex;
ldY *= ldY;
}
while (ldZ / ldX >= ldY);
}
while (ldX / ldZ >= 2.0L)
{
theIndex = 1.0L;
ldY = 2.0L;
do
{
theSum += theIndex;
ldZ *= ldY;
theIndex += theIndex;
ldY *= ldY;
}
while (ldX / ldZ >= ldY);
}
ldX = (ldX - ldZ) / ldZ;
if (ldX > 0.0L)
{
theIndex = 1.0L;
do
{
theIndex /= 2.0L;
ldY = ldX * ldX;
ldZ = ldX + ldX;
ldX = ldY + ldZ;
}
while (ldX < 1.0L);
ldX = (ldY + --ldZ) / 2.0L;
qSum = theIndex;
if (ldX > 0.0L)
{
do
{
lastSum = pSum;
do
{
theIndex /= 2.0L;
ldY = ldX * ldX;
ldZ = ldX + ldX;
ldX = ldY + ldZ;
}
while (ldX < 1.0L);
ldX = (ldY + --ldZ) / 2.0L;
pSum += theIndex;
}
while (ldX > 0.0L && pSum != lastSum);
}
pSum += qSum;
}
*theInt = theSum;
return (pSum);
}
/*
TwoToTheX returns 2 to the power of a long double argument.
If FLT_RADIX is 2, the value returned should be accurate to
LDBL_MANT_DIG bits.
*/
long double
TwoToTheX (long double theX)
{
long double ldX, theIndex, theProd, pProd;
long double thePower, prevPower, lastApprox;
if (theX == -HUGE_VALL)
return (errno = ERANGE, 0.0L);
if (theX == HUGE_VALL)
return (errno = ERANGE, HUGE_VALL);
ldX = theX;
theProd = 1.0L;
while (ldX < 0.0L)
{
theIndex = 1.0L;
thePower = 2.0L;
do
{
theProd /= thePower;
ldX += theIndex;
if (theProd == 0.0L)
return (errno = ERANGE, theProd);
theIndex += theIndex;
thePower *= thePower;
}
while (-ldX >= theIndex);
}
while (ldX >= 1.0L)
{
theIndex = 1.0L;
thePower = 2.0L;
do
{
theProd *= thePower;
ldX -= theIndex;
if (theProd == HUGE_VALL)
return (errno = ERANGE, theProd);
theIndex += theIndex;
thePower *= thePower;
}
while (ldX >= theIndex);
}
if (ldX > 0.0L)
{
pProd = 0.0L;
thePower = theIndex = 1.0L;
do
{
prevPower = thePower;
thePower = 0.0L;
do
{
lastApprox = thePower;
thePower += (prevPower - thePower) / (thePower + 1.0L);
thePower /= 2.0L;
}
while (thePower != lastApprox);
theIndex /= 2.0L;
if (ldX >= theIndex)
{
pProd += (pProd * thePower) + thePower;
ldX -= theIndex;
}
}
while (ldX > 0.0L && (thePower + 1.0L) > 1.0L);
theProd += theProd * pProd;
if (theProd == HUGE_VALL)
errno = ERANGE;
}
return (theProd);
}
/*
These macros are used in calls to HexTimesX.
Cantrell's coefficients are rational numbers, but can't be precisely
represented as long doubles. The other constants are irrational.
The "S" at the beginning of S_SQUARE_ROOT_OF_PI and S_BASE_2_LOG_OF_E
is meant to prevent possible conflicts. MANT_OF_BASE_2_LOG_SQRT_2PI
is the mantissa of the base 2 log of the square root of 2pi.
*/
#if LDBL_DIG < 19
#define MANT_OF_BASE_2_LOG_SQRT_2PI "0.536439A4C6EFBAD86E"
#define S_SQUARE_ROOT_OF_PI "1.C5BF891B4EF6AA79C"
#define S_BASE_2_LOG_OF_E "1.71547652B82FE1778"
#define CANTRELL_COEFFICIENT_EIGHT "1.8F27478402B14D893"
#define CANTRELL_COEFFICIENT_SEVEN "0.28E791360B93BC5F5D"
#define CANTRELL_COEFFICIENT_SIX "2.0D5E32BC6D9F5D78E"
#define CANTRELL_COEFFICIENT_FIVE "0.372642BE5EE71DC552"
#define CANTRELL_COEFFICIENT_FOUR "2.FDCDC5670AA34302A"
#define CANTRELL_COEFFICIENT_THREE "0.54C0965CAE8514023D"
#define CANTRELL_COEFFICIENT_TWO "5.734CE827B8D8E9839"
#define CANTRELL_COEFFICIENT_ONE "0.B6679DB440592CF7C4"
#else
#define MANT_OF_BASE_2_LOG_SQRT_2PI "0.536439A4C6EFBAD86D9727840544713A5C"
#define S_SQUARE_ROOT_OF_PI "1.C5BF891B4EF6AA79C3B0520D5DB9383FE"
#define S_BASE_2_LOG_OF_E "1.71547652B82FE1777D0FFDA0D23A7D11D"
#define CANTRELL_COEFFICIENT_EIGHT "1.8F27478402B14D89312C473E2C34ACB43"
#define CANTRELL_COEFFICIENT_SEVEN "0.28E791360B93BC5F5C9D6B186484A0EBEA"
#define CANTRELL_COEFFICIENT_SIX "2.0D5E32BC6D9F5D78DF50633852B8005E1"
#define CANTRELL_COEFFICIENT_FIVE "0.372642BE5EE71DC5524A4434F49AE415FB"
#define CANTRELL_COEFFICIENT_FOUR "2.FDCDC5670AA343029D7681AE3E727CB38"
#define CANTRELL_COEFFICIENT_THREE "0.54C0965CAE8514023D78641D8EE67C0D33"
#define CANTRELL_COEFFICIENT_TWO "5.734CE827B8D8E983955383CB54E551F73"
#define CANTRELL_COEFFICIENT_ONE "0.B6679DB440592CF7C43636F98A326A2641"
#endif
/*
Log2OfGamma returns the mantissa of the base 2 logarithm of the
(complete) Gamma function of theX. The integer part of the logarithm
is returned in the variable pointed to by theInt; this will be the
largest integer not greater than the logarithm. E.g., if the logarithm
is 2.1, Log2OfGamma will return 0.1 and 2.0; if the logarithm is -0.1,
it will return 0.9 and -1.0.
Log2OfGamma uses an approximation developed by David W. Cantrell;
thanks to Hugo Pfoertner for pointing me to it. See:
http://mathforum.org/epigone/sci.math.num-analysis/zonchahee/
*/
long double
Log2OfGamma (long double theX, long double *theInt)
{
long double theValue, ldX, thePochhammer;
long double fractA, intA, fractB, intB;
if (theX < 0.0L)
{
errno = EDOM;
*theInt = NAN;
return (0.0L);
}
if (theX == 0.0L || theX == HUGE_VALL)
{
errno = ERANGE;
*theInt = HUGE_VALL;
return (0.0L);
}
ldX = theX;
/* If it's easy, calculate Gamma the obvious way. */
if (ldX <= 9.0L && LDFract (ldX * 2.0L) == 0.0L)
{
theValue = 1.0L;
while (--ldX > 0.0L)
theValue *= ldX;
if (ldX < 0.0L)
theValue = HexTimesX (S_SQUARE_ROOT_OF_PI, theValue);
return (Log2OfX (theValue, theInt));
}
/* Otherwise, use Cantrell's approximation. */
thePochhammer = 1.0L;
while (ldX <= 9.0L)
thePochhammer *= ldX++;
ldX -= 0.5L;
theValue = HexTimesX (CANTRELL_COEFFICIENT_EIGHT, ldX);
theValue = HexTimesX (CANTRELL_COEFFICIENT_SEVEN, ldX) + (1.0L / theValue);
theValue = HexTimesX (CANTRELL_COEFFICIENT_SIX, ldX) + (1.0L / theValue);
theValue = HexTimesX (CANTRELL_COEFFICIENT_FIVE, ldX) + (1.0L / theValue);
theValue = HexTimesX (CANTRELL_COEFFICIENT_FOUR, ldX) + (1.0L / theValue);
theValue = HexTimesX (CANTRELL_COEFFICIENT_THREE, ldX) + (1.0L / theValue);
theValue = HexTimesX (CANTRELL_COEFFICIENT_TWO, ldX) + (1.0L / theValue);
theValue = HexTimesX (CANTRELL_COEFFICIENT_ONE, ldX) + (1.0L / theValue);
theValue = (24.0L * ldX) + (1.0L / theValue) - 0.5L;
theValue = (1.0L / (1.0L + (1.0L / theValue))) / thePochhammer;
fractA = Log2OfX (theValue, &intA);
fractB = HexTimesX (MANT_OF_BASE_2_LOG_SQRT_2PI, 1.0L);
intB = 1.0L;
if (fractA + fractB >= 1.0L)
{
--fractB;
++intB;
}
fractA += fractB;
intA += intB;
fractB = ATimesLog2OfX (ldX, ldX, &intB);
if (fractA + fractB >= 1.0L)
{
--fractB;
++intB;
}
fractA += fractB;
intA += intB;
intB = HexTimesX (S_BASE_2_LOG_OF_E, ldX);
fractB = LDFract (intB);
intB -= fractB;
if (fractA - fractB < 0.0L)
{
--fractB;
++intB;
}
fractA -= fractB;
intA -= intB;
theValue = fractA;
*theInt = intA;
if (intA == HUGE_VALL)
errno = ERANGE;
return (theValue);
}
/*
IncompleteGammaFactor is called by the incomplete gamma functions.
Where "^" means "to the power of," and "Gamma" stands for an
upper-case Greek letter gamma, this function calculates
(x^a)/(Gamma(a)(e^x)).
*/
long double
IncompleteGammaFactor (long double theA, long double theX)
{
long double fractA, fractB, intA, intB;
fractA = ATimesLog2OfX (theA, theX, &intA);
fractB = Log2OfGamma (theA, &intB);
if (fractA - fractB < 0.0L)
{
--fractB;
++intB;
}
fractA -= fractB;
intA -= intB;
intB = HexTimesX (S_BASE_2_LOG_OF_E, theX);
fractB = LDFract (intB);
intB -= fractB;
if (fractA - fractB < 0.0L)
{
--fractB;
++intB;
}
fractA -= fractB;
intA -= intB;
return (TwoToTheX (fractA) * TwoToTheX (intA));
}
/* Don't call IncompleteGammaBySeries directly; call IncompleteGamma. */
double
IncompleteGammaBySeries (double theA, double theX)
{
long double theSum, lastSum, aValue, theProd;
long double ldA, ldX;
ldA = theA;
ldX = theX;
aValue = ldA;
theSum = theProd = 1.0L / aValue;
do
{
lastSum = theSum;
theProd /= ++aValue;
theProd *= ldX;
theSum += theProd;
}
while (theSum != lastSum);
return (1.0L - (theSum * IncompleteGammaFactor (ldA, ldX)));
}
/* Don't call IncompleteGammaByFraction directly; call IncompleteGamma. */
double
IncompleteGammaByFraction (double theA, double theX)
{
long double aValue, theB, theC, theD, ldA, ldX;
long double theIndex, theProd, lastProd;
ldA = theA;
ldX = theX;
theB = 1.0L + (ldX - ldA);
theD = theProd = 1.0L / theB;
aValue = ldA - 1.0L;
theB += 2.0L;
theC = theB;
theD = 1.0L / (theB + (aValue * theD));
theProd *= theC * theD;
theIndex = 1.0L;
do
{
lastProd = theProd;
++theIndex;
aValue = theIndex * (ldA - theIndex);
theB += 2.0L;
theC = theB + (aValue / theC);
theD = 1.0L / (theB + (aValue * theD));
theProd *= theC * theD;
}
while (theProd != lastProd);
return (theProd * IncompleteGammaFactor (ldA, ldX));
}
/*
IncompleteGamma calculates the regularized incomplete gamma function.
Where "Gamma" stands for an upper-case Greek letter gamma, IncompleteGamma
calculates Gamma(a,x)/Gamma(a). W. H. Press, et al (in Numerical Recipes
in C) write this function as Q(a,x). It is the complement of P(a,x); i.e.
Q(a,x) = 1 - P(a,x). Though the functions here are not from Numerical
Recipes in C, they do implement formulae (6.2.5) and (6.2.7) from that
book. Cf. formulae (2) and (4) in "Computing the incomplete Gamma function
to arbitrary precision" (2003: LNCS 2667) by Serge Winitzki. The continued
fraction expansion [Press, et al (6.2.6) and Winitzki (4)] is due to Legendre.
*/
double
IncompleteGamma (double theA, double theX)
{
if (theA <= 0.0 || theX < 0.0)
return (errno = EDOM, 0.0);
if (theA == HUGE_VAL)
return (errno = ERANGE, theX / HUGE_VAL);
if (theX == HUGE_VAL)
return (0.0);
if (theX == 0.0)
return (1.0);
if (theX < theA)
return (IncompleteGammaBySeries (theA, theX));
return (IncompleteGammaByFraction (theA, theX));
}
/*
ChiSquaredPValue returns the probability of a chi-squared statistic
greater than or equal to theChiSquared.
*/
double
ChiSquaredPValue (double degOfFreedom, double theChiSquared)
{
return (IncompleteGamma (degOfFreedom / 2.0, theChiSquared / 2.0));
}
/* Here ends this small library of math functions. */