Powers of Primes Dividing Factorials

It's often useful to know the highest power of a given prime p that
divides n! for a given positive integer n.  This enables us to 
determine the powers of primes dividing multinomial coefficients,
and so on.  It turns out that there is a nice formula, apparently
first discussed by Legendre, involving the sum of the digits of n
when n is expressed in the base p.

Writing the number n in the base p we have

      n = n0 + n1 p + n2 p^2 + n3 p^3 + ... + nk p^k

The number of powers of p dividing n! is just the sum of the powers
of p dividing each of the numbers n, n-1, n-2, ..., 1.  Of these
numbers, those that are divisible by at least one power of p are
all the multiples of p less than n-n0, of which there are exactly
(n - n0) / p.  Of these numbers, some are divisible by TWO powers 
of p, namely, all the multiples of p^2 less than n - n0 - n1 p, of
which there are exactly (n - n0 - n1 p) / p^2.  Continuing on in this
way we find that the total number of powers of p dividing n! is the
sum

   n - n0     n - n0 - n1 p     n - n0 - n1 p - n2 p^2
   ------  +  -------------  +  ----------------------  +  ...
      p            p^2                    p^3


                         n - n0 - n1 p - n2 p^2 - ... - n_k p^k
                      +  --------------------------------------
                                          p^(k+1)

Placing all these fractions on a common denominator of p^(k+1), we
arrive at

 (p^k+...+p+1)n - (p^k+...+p+1)n0 - (p^k + ... + p)n1 + ... - p^k nk
 ------------------------------------------------------------------
                               p^(k+1)

Noting that the coefficients of n and n0 in the numerator are 
the partial geometric sum [p^k + p^(k-1) + ... + p + 1] and the 
coefficients of n1,n2,... are truncated copies of this, we see 
that the expression can be simplified if we multiply through by 
p-1.  This leads us to set the above sum equal to X/(p-1) and 
for X.  If we let j denote k+1, this gives

       (p^j - 1)n - (p^j - 1)n0 - (p^j - p)n1 - ... - (p^j - p^k)nk
 X  =  ------------------------------------------------------------
                                   p^j

                                       n - (n0 + n1 p + ... + nk p^k)
       =  n - (n0 + n1 + ... + nk)  -  ------------------------------
                                                       p^j

The expression inside the parentheses in the numerator of the right-
hand fraction is simply n, so the fraction vanishes and we are left
with X = n - s_p(n) where s_p(k) denotes the sum of the base-p digits 
of k.  Therefore, if we let f_p(k) denote the number of powers of p
dividing k, we have proven Legendre's result that

                             n - s_p(n)
                f_p(n!)  =  ------------
                               p - 1

Obviously we can then express the powers of p dividing the binomial
coefficient C(m+n,n) = (m+n)!/((m!)(n!)) as (omitting the subscripts)
              _         _
             |  / m+n \  |     s(m) + s(n) - s(m+n)
            f| (       ) |  =  --------------------
             |_ \  n  / _|              p - 1

This is interesting because it shows that the number of powers of p
dividing C(m+n,n) is a measure of the extent to which the operations
of addition and s() are not distributive.  For each of the infinitely 
many primes p that do not divide C(m+n,n) we have

                   s(m+n) = s(m) + s(n)

but for a prime p that divides C(m+n,n) this equality doesn't hold, 
and the extent to which it doesn't hold equals (p-1) times the number
of powers of p that divide C(m+n,n).  Of course, the above identity
also allows us to express the binomial coefficients as the product
of powers of primes

       / m+n \             [s_p(m)+s_p(n)-s_p(m+n)]/(p-1)
      (       )  =  PROD  p
       \  n  /

where the product is evaluated over all primes.

We can also easily verify that, for any fixed prime p, consecutive
values of s() are related to consecutive values of f() as follows

                     /      1          if  f(k) < = f(k-1)
  s(k) - s(k-1)  =  (
                     \ 1 - (p-1) f(k)  if  f(k) > f(k-1)

This is clear simply by examining the sequence of base-p expressions,
such as (in the base 5)

                 ...
                 032
                 033   s=s+1
                 034   s=s+1
                 040   s=s+1-1(4)
                 041   s=s+1
                 042   s=s+1
                 043   s=s+1
                 044   s=s+1
                 100   s=s+1-2(4)
                 101   s=s+1
                 102   s=s+1
                    etc.

For the steps in between multiples of p the value of f(k) is zero
and the sum of the digits increases by 1.  When going from a non-
multiple of p to a multiple of p^t, we see that f(k) increases by t
and s(k) changes by 1 - (p-1) f(k).