Category : C Source Code
Archive   : SNIP1292.ZIP
Filename : FACTORYL.C

 
Output of file : FACTORYL.C contained in archive : SNIP1292.ZIP
/*
** FACTORYL.C
**
** Original Copyright 1992 by Bob Stout as part of
** the MicroFirm Function Library (MFL)
**
** This subset version is hereby donated to the public domain.
*/

#include
#include
#include

#define dfrac(x) ((x)-dround(x))

/* Use #defines for Permutations and Combinations */

#define log10P(n,r) (log10factorial(n)-log10factorial((n)-(r)))
#define log10C(n,r) (log10P((n),(r))-log10factorial(r))

#ifndef PI
#define PI 3.14159265358979323846
#endif

#define SQRT2PI sqrt(2 * PI)
#define ONESIXTH (1.0/6.0)

/*
** DROUND.C - Rounds a double to the nearest whole number
** public domain by Ross Cottrell
*/

double dround(double x)
{
assert(1 == FLT_ROUNDS);
x += 1.0 / DBL_EPSILON;
return x - 1.0 / DBL_EPSILON;
}

/*
** log10factorial()
**
** Returns the logarithm (base 10) of the factorial of a given number.
** The logarithm is returned since this allows working wil extremely
** large values which would otherwise overflow the F.P. range.
**
** Parameters: 1 - Number whose factorial to return.
**
** Returns: log10() of the passed value, -1.0 if error
**
*/

double log10factorial(double N)
{
if (N < 40) /* Small, explicitly compute */
{
int i;
double f;

for (i = 1, f = 1.0; i <= (int)N; ++i)
f *= i;
return log10(f);
}
else /* Large, use approximation */
{
return log10(SQRT2PI)+((N + 0.5) *
(log10(sqrt(N * N + N + ONESIXTH) / exp(1))));
}
}

#ifdef TEST

#include
#include

void main(int argc, char *argv[])
{
double f, lf;
char *dummy;

while (--argc)
{
f = strtod((const char *)(*(++argv)), &dummy);
lf = log10factorial(f);
if (171.0 > f)
printf("%.14g! = %.14g\n", f, pow(10.0, lf));
else
{
printf("%.14g! = %.14ge+%ld\n", f,
pow(10.0, dfrac(lf)), (long)dround(lf));
}
}
lf = log10C(1000000,750000);
printf("\nJust to dazzle with you with big numbers:\n"
"C(1000000,750000) = %.14ge+%ld\n",
pow(10.0, dfrac(lf)), (long)dround(lf));
lf = log10P(1000000,750000);
printf("\n...once more:\n"
"P(1000000,750000) = %.14ge+%ld\n",
pow(10.0, dfrac(lf)), (long)dround(lf));
}

#endif /* TEST */


  3 Responses to “Category : C Source Code
Archive   : SNIP1292.ZIP
Filename : FACTORYL.C

  1. Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!

  2. This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.

  3. But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/