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;
}

Milliard

why start from base 2? because j % 1 will neer have a remainder. – Cheery

@Marc B: Variable p counts the number of times that the value of variable j is divisible by the value of variable i, and it does so only for prime values of variable i (due to the j /= i). The loop starts from 2 because this is the smallest prime factor. – Tiliaceous

@MarcB That's not really a good explanation, since 1 is in fact a divisor, so why not consider it? The reason is that we only want prime divisors – Commend

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.

Commend

I think the code in the question does unnecessary things. We simply need to find the minimum of the highest power of 2 and the highest power of 5 in the factorial of the number to find the number of trailing zeros in the factorial. – Hurl

@Hurl The code works for general bases. Of course what you describe is exactly what that algorithm would do for base = 10 = 2^1 * 5^1 – Commend

Answer 2

Basically, it finds prime factors of base, where i is prime and i^p is a factor of base, and then figures out how many i factors will exist in n!, divides that by p, and tracks the minimum number of that result over all the prime factors of base.

So to answer the questions in the code:

1. why do we start from 2
   • because 2 is the smallest prime
2. what is p?
   • p is the 'power' for the prime factor i of the base (so i^p is a factor of base)
3. What is the maths behind this loop
   • the loop is counting (c) the number of factors of i in n!

Nara