home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / snip9707.zip / FACTORYL.C < prev    next >
C/C++ Source or Header  |  1997-07-05  |  3KB  |  105 lines
  1. /* +++Date last modified: 05-Jul-1997 */
  2.  
  3. /*
  4. **  FACTORYL.C
  5. **
  6. **  Original Copyright 1992 by Bob Stout as part of
  7. **  the MicroFirm Function Library (MFL)
  8. **
  9. **  The user is granted a free limited license to use this source file
  10. **  to create royalty-free programs, subject to the terms of the
  11. **  license restrictions specified in the LICENSE.MFL file.
  12. **
  13. **  Uses DBLROUND.C, also in SNIPPETS
  14. */
  15.  
  16. #include <float.h>
  17. #include "snipmath.h"
  18.  
  19. #define dfrac(x) ((x)-dround(x))
  20.  
  21. #define SQRT2PI sqrt(2 * PI)
  22. #define ONESIXTH (1.0/6.0)
  23.  
  24. /*
  25. **  log10factorial()
  26. **
  27. **  Returns the logarithm (base 10) of the factorial of a given number.
  28. **  The logarithm is returned since this allows working with extremely
  29. **  large values which would otherwise overflow the F.P. range.
  30. **
  31. **  Parameters: 1 - Number whose factorial to return.
  32. **
  33. **  Returns: log10() of the passed value, -1.0 if error
  34. **
  35. **  Limitations: Cannot return 0! since log(0) is undefined.
  36. */
  37.  
  38. double log10factorial(double N)
  39. {
  40.       double dummy;
  41.  
  42.       if ((N < 1.0) || (0.0 != modf(N, &dummy)))
  43.             return -1.0;
  44.       if (N < 40)                         /* Small, explicitly compute  */
  45.       {
  46.             int i;
  47.             double f;
  48.  
  49.             if (0.0 == N)
  50.                   return N;
  51.             for (i = 1, f = 1.0; i <= (int)N; ++i)
  52.                   f *= i;
  53.             return log10(f);
  54.       }
  55.       else                                /* Large, use approximation   */
  56.       {
  57.             return log10(SQRT2PI)+((N + 0.5) *
  58.                   (log10(sqrt(N * N + N + ONESIXTH) / exp(1.0))));
  59.       }
  60. }
  61.  
  62. #ifdef TEST
  63.  
  64. #include <stdio.h>
  65. #include <stdlib.h>
  66.  
  67. main(int argc, char *argv[])
  68. {
  69.       double f, lf;
  70.       char *dummy;
  71.  
  72.       while (--argc)
  73.       {
  74.             f  = strtod((const char *)(*(++argv)), &dummy);
  75.             if (0.0 == f)
  76.             {
  77.                   puts("0! = 0");
  78.                   continue;
  79.             }
  80.             if (-1.0 == (lf = log10factorial(f)))
  81.             {
  82.                   printf(">>> ERROR: %g! is not a valid expression\n", f);
  83.                   continue;
  84.             }
  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(1000000L,750000L);
  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(1000000L,750000L);
  98.       printf("\n...once more:\n"
  99.             "P(1000000,750000) = %.14ge+%ld\n",
  100.             pow(10.0, dfrac(lf)), (long)dround(lf));
  101.       return EXIT_SUCCESS;
  102. }
  103.  
  104. #endif /* TEST */
  105.