double  Do_logfac (long n)  /* pass by value ANSI way */

/*  Given positive integer "n", the function calculates log_2 n! */
{ /*  Begin function Do_logfac() */
  static double ln2fac [16] =
  { 0.00000,  0.00000, 1.00000, 2.585,
    4.585,    6.907,   9.4915,  12.3,
    15.3,     18.47,   21.79,   25.25,
    28.83,    32.53,   36.34,   40.25
  };
  long i;
  double n_flt, result;

  if (n < 0)
  {
     fprintf (stderr, "Bad logfac input %9d\n", n);
     return(ln2fac[0]);
  }
  else /* valid "n => 0" */
  {
     if (n < 16)  /* small "n"? */
     {     i = n;
          return (ln2fac[i]);
     }  /* Table look-up path was taken */
     else
     {
          /* Use Stirling's Formula (one of several versions)
           *  Most popular version of Stirling's formula:
           *  n!  =(approx) Sqrt[2 Pi] n^{n+0.5) e^{-n}.
           * Value "e" is the base of the natural logarithm.
           * Notation a^{b} means:  value "a" raised to the power "b";
           * i.e., ^{b} denotes "b" as the exponent.
           * Above version of Stirling's Formula, modified to log_2  */

        n_flt = (double) n; /* type converstion to floating point */
        result = (n_flt + 0.5) * log (n_flt);
        result = (result - n_flt) + 0.9189;  /*  asymptotic fudge factor
  */       result = result * 1.442695;    /* convert to log base 2 */
        return (result);
     } /* End Stirling Formula */
  } /* End clause where n > 0 */
} /*  End Function Do_logfac() */
 /********************************************************/

 /*  Use of log factorial function (log_2 n!) for Enumerative CodeLength
 enl

 The enumerative code length applies to a data file (or data sequence)
 D comprised of a succession (or sequence) of data values (typical
 value"i").  In general, the values range from 0 to N-1.  For data
 bytes, the values range from 0 to 255.  The number of occurences of
 each value "i" is that value's frequency "hst[i]".  The sum comprised
 of each value's frequency hst[i] is the total number of values (bytes)
 in the data sequence D.  A histogram is a vector arranged as a
 sequence of frequencies hst[i] corresponding to the numerical ordering
 of the values themselves.  For example, given the 8-bit binary integer
 data values, the ordering is 0, 1, ..., i, ..., N-1.

 The total count obtained by summing the N histogram values may be
 viewed as the "length" L (number of elements in) data sequence D.
 Thus the sum of the 256 numbers of the histogram vector, since the
 numbers themselves are *occurrences*, must equal L, the total
 occurrences of values in the data.  The values themselves (0, 1, ...,
 255) are accounted for by their position within the 256-element
 histogram vector.  Thus value 0 is represented as the first vector
 postion, and the number of 0s found in the data sequence is found
 there.

 As an example, a data file of 10,000 8-bit bytes has a histogram on
 the ordered values 0, 1, ..., 255 since these are the numberical
 values assigned to 8-bit (binary) bytes 0000000, 00000001, ...,
 11111111.  The sum of the histogram frequencies in this example should
 be 10,000:  otherwise a mistake must have occurred.

 The relative frequency of value "i" is calculated as the frequecy of
 "i" divided by the total number of values in the sequence.

 Given:  Histogram hst[i] of N values 0, 1, ..., i, ..., N-1.

 The "elements" of vector hst[] represent the counts or occurrences for
 each value i (values i range from 0 through N-1).  Element 0 of the
 vector 'hst' represents the number of occurences of value 0 in the
 data from which the histogram was created.  So hst[0] is the frequency
 of the first, and hst[N-1] is the frequency of the last value of the
 histogram for the data sequence.

 Variable "totct" is the sum formed by adding each count:
                  totct =  sum (from i =0 to i= N-1) of "hst[i]"
 */

 /* ****************************************************** */
 double  Sum_logfac_hst(long hst_ary[MAXval])

 /* -log2 of [(N-1)! x n1! x ... x nN!) ] / [(N + Totcount)!], N is nr
 of symbols in the alphabet, and n1 is count of symbol 1, etc. */
 /*  Generate a bit-per-pixel array per context: bpp_ary */
 {
         long     subtot;
         double  Big_logfac, Small_logfac;
         long     i, j, k;
         double  sum = 0.0;
         for (i=0; i <=  0; i++)  { /* [i] for future use in contexts */
           subtot = 0;
           Small_logfac = 0.0;
           for (k=0; k<=(MAXval-1); k++) {
                   subtot += hst_ary/*[i]*/[k];
                   Small_logfac += Do_logfac(hst_ary/*[i]*/[k]);
           }
           /* Next: adjust per alphabet size (MAXval) */
           Small_logfac += Do_logfac(MAXval-1);
           Big_logfac = Do_logfac(subtot + (MAXval-1));
           sum += (Big_logfac - Small_logfac);
                 }
         return (sum);
 }