home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / s / snip1292.zip / FACTORYL.C < prev    next >
C/C++ Source or Header  |  1992-12-26  |  3KB  |  104 lines
  1. /*
  2. **  FACTORYL.C
  3. **
  4. **  Original Copyright 1992 by Bob Stout as part of
  5. **  the MicroFirm Function Library (MFL)
  6. **
  7. **  This subset version is hereby donated to the public domain.
  8. */
  9.  
  10. #include <math.h>
  11. #include <float.h>
  12. #include <assert.h>
  13.  
  14. #define dfrac(x) ((x)-dround(x))
  15.  
  16. /* Use #defines for Permutations and Combinations     */
  17.  
  18. #define log10P(n,r) (log10factorial(n)-log10factorial((n)-(r)))
  19. #define log10C(n,r) (log10P((n),(r))-log10factorial(r))
  20.  
  21. #ifndef PI
  22.  #define PI 3.14159265358979323846
  23. #endif
  24.  
  25. #define SQRT2PI sqrt(2 * PI)
  26. #define ONESIXTH (1.0/6.0)
  27.  
  28. /*
  29. **  DROUND.C - Rounds a double to the nearest whole number
  30. **  public domain by Ross Cottrell
  31. */
  32.  
  33. double dround(double x)
  34. {
  35.       assert(1 == FLT_ROUNDS);
  36.       x += 1.0 / DBL_EPSILON;
  37.       return x - 1.0 / DBL_EPSILON;
  38. }
  39.  
  40. /*
  41. **  log10factorial()
  42. **
  43. **  Returns the logarithm (base 10) of the factorial of a given number.
  44. **  The logarithm is returned since this allows working wil extremely
  45. **  large values which would otherwise overflow the F.P. range.
  46. **
  47. **  Parameters: 1 - Number whose factorial to return.
  48. **
  49. **  Returns: log10() of the passed value, -1.0 if error
  50. **
  51. */
  52.  
  53. double log10factorial(double N)
  54. {
  55.       if (N < 40)                         /* Small, explicitly compute  */
  56.       {
  57.             int i;
  58.             double f;
  59.  
  60.             for (i = 1, f = 1.0; i <= (int)N; ++i)
  61.                   f *= i;
  62.             return log10(f);
  63.       }
  64.       else                                /* Large, use approximation   */
  65.       {
  66.             return log10(SQRT2PI)+((N + 0.5) *
  67.                   (log10(sqrt(N * N + N + ONESIXTH) / exp(1))));
  68.       }
  69. }
  70.  
  71. #ifdef TEST
  72.  
  73. #include <stdio.h>
  74. #include <stdlib.h>
  75.  
  76. void main(int argc, char *argv[])
  77. {
  78.       double f, lf;
  79.       char *dummy;
  80.  
  81.       while (--argc)
  82.       {
  83.             f  = strtod((const char *)(*(++argv)), &dummy);
  84.             lf = log10factorial(f);
  85.             if (171.0 > f)
  86.                   printf("%.14g! = %.14g\n", f, pow(10.0, lf));
  87.             else
  88.             {
  89.                   printf("%.14g! = %.14ge+%ld\n", f,
  90.                         pow(10.0, dfrac(lf)), (long)dround(lf));
  91.             }
  92.       }
  93.       lf = log10C(1000000,750000);
  94.       printf("\nJust to dazzle with you with big numbers:\n"
  95.             "C(1000000,750000) = %.14ge+%ld\n",
  96.             pow(10.0, dfrac(lf)), (long)dround(lf));
  97.       lf = log10P(1000000,750000);
  98.       printf("\n...once more:\n"
  99.             "P(1000000,750000) = %.14ge+%ld\n",
  100.             pow(10.0, dfrac(lf)), (long)dround(lf));
  101. }
  102.  
  103. #endif /* TEST */
  104.