How does this code find the number of trailing zeros from any base number factorial?

Question

The code below works perfectly but I would like someone to explain to me the mathematics behind it. Basically, how does it work?

#include <stdio.h> 
#include <stdlib.h>  /* atoi */

#define min(x, y) (((x) < (y)) ? (x) : (y))

int main(int argc, char* argv[])
{
   const int base = 16;
   int n,i,j,p,c,noz,k;

   n = 7;  /* 7! = decimal 5040 or 0x13B0  - 1 trailing zero */  
   noz = n;
   j = base;
   /* Why do we start from 2 */
   for (i=2; i <= base; i++)
   {
      if (j % i == 0)
      {   
         p = 0;  /* What is p? */
         while (j % i== 0)
         {
            p++;
            j /= i;
         }
         c = 0;
         k = n;
         /* What is the maths behind this while loop? */
         while (k/i > 0)
         {
            c += k/i;
            k /= i;
         }
         noz = min(noz, c/p);
      }
   }
   printf("%d! has %d trailing zeros\n", n, noz);

   return 0;
}

Answer 1

Note that the problem is equivalent to finding the highest power of base which divides n!.

If the base was prime (let's call it p), we could use a theorem from number theory to compute the highest power of p that divides n!:

Let's extract the part of your code that does this into a function:

int maximum_power_of_p_in_fac(int p, int n) {
    int mu = 0;
    while (n/p > 0) {
        mu += n/p;
        n /= p;
    }
    return mu;
}

Now what happens if base is a prime power? Let's say we have base = p^q. Then if μ is the highest power of p which divides n! and r = floor(μ/q), we have

(p^q)^r = p^(qr) divides p^μ divides n!

and

(p^q)^(r+1) = p^(q(r+1)) >= p^(μ+1) does not divide n!

So r is the maximum power of p^q in n!. Let's write a function for that as well:

int maximum_power_of_pq_in_fac(int p, int q, int n) {
    return maximum_power_of_p_in_fac(p, n) / q;
}

So what if base is general number? Let's say

base = p₁^q₁ p₂^q₂ ... p_m^q_m

(this is the unique prime factorization of base). Then we just solve the problem for all p_i^q_i and take the minimum of those:

int maximum_power_of_base_in_fac(int base, int n) {
    int res = infinity;
    for every factor p^q in the prime factorization in base:
       res = min(res, maximum_power_of_pq_in_fac(p,q,n));
    return res;
}

How to factorize base? Well we can just use trial division, like your example code does. We start by checking whether 2 is a prime factor. If it is, we compute maximum_power_of_pq_in_fac and divide base by 2 until it is no longer divisible by 2. Then we proceed with the next candidate factor:

void factorize(int base) {
    for (int p = 2; p <= base; ++p) {
        if (base % p == 0) { // if base is divisible by p, p is a prime factor
            int q = 0;
            while (base % p == 0) { // compute q and get rid of all the p factors
                q++;
                base /= p;
            }
            // do something with factor p^q
        }
        // proceed with next candidate divisor
    }
}

By examining the code carefully, you will find that it contains all the above elements, only put together in a single loop, which is a bit confusing.

UPDATE: In case you are interested, the algorithm you presented has complexity O(base * log n). You can easily make it O(sqrt(base) * log n) by adapting the prime factorization routine slightly:

void factorize(int base) {
    for (int p = 2; p*p <= base; ++p) { // only check prime factors up to sqrt(base)
        // ... same as before
    }
    if (base) {
        // there is exactly one prime factor > sqrt(base). 
        // It certainly has multiplicity 1.

        // process prime factor base^1
    }   
}

And of course you can use any other more sophisticated prime factorization algorithm if you want to speed things up even more.