Trailing zeroes of falling factorials

29 Apr 2017

def exp_of_b_in(b, n):
  if n < b:
    return 0
  elif n % b != 0:
    return 0
  else:
    return 1 + exp_of_b_in(b, n / b)

print [exp_of_b_in(2, i) for i in range(1, 20)]

def falling_factorial(x, n):
  if n == 0:
    return 1
  else:
    return x * falling_factorial(x - 1, n - 1)

print [[falling_factorial(x, n) for n in range(2, 4)] for x in range(3, 6)]

def trailing_zeroes_of_falling_factorial(x, n):
  return exp_of_b_in(10, falling_factorial(x, n))

print [[trailing_zeroes_of_falling_factorial(x, n) for n in range(2, 4)] for x in range(3, 6)]

Is there a systematic way to transform trailing_zeroes_of_falling_factorial to avoid creating very large intermediate values?

Let's try a standard refactoring, recursion into iteration, on the recursive function exp_of_b_in. For safety, let's take the existing code as a specification, and use it to generate some clamp tests. To create a clamp test, you just generate some characteristic output of the thing, and then copy and paste the output back into the source code, and assert that the same input should generate the same output, even after refactoring.

I'm not certain this will work, but I have to try something.

def exp_of_b_in(b, n):
  if n < b:
    return 0
  elif n % b != 0:
    return 0
  else:
    return 1 + exp_of_b_in(b, n / b)

assert [exp_of_b_in(2, i) for i in range(1, 100)] == [0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5, 0, 1, 0]
assert [exp_of_b_in(3, i) for i in range(1, 100)] == [0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 3, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 3, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 4, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 2]

Getting it to work took some trial and error (I forgot that the opposite of < is >= not >), but the clamp tests caught my bugs, and they are extensive enough that I am now convinced that this is right.

def exp_of_b_in(b, n):
  accumulator = 0
  while n >= b and n % b == 0:
    accumulator += 1
    n = n / b
  return accumulator

assert [exp_of_b_in(2, i) for i in range(1, 100)] == [0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5, 0, 1, 0]
assert [exp_of_b_in(3, i) for i in range(1, 100)] == [0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 3, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 3, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 4, 0, 0, 1, 0, 0, 1, 0, 0, 2, 0, 0, 1, 0, 0, 1, 0, 0, 2]

That worked. What if we did the same transformation to the falling_factorial function?

def falling_factorial(x, n):
  if n == 0:
    return 1
  else:
    return x * falling_factorial(x - 1, n - 1)

print [[falling_factorial(x, n) for n in range(1, 10)] for x in range(1, 20)]
assert [[falling_factorial(x, n) for n in range(1, 10)] for x in range(1, 20)] == [
  [1, 0, 0, 0, 0, 0, 0, 0, 0],
  [2, 2, 0, 0, 0, 0, 0, 0, 0],
  [3, 6, 6, 0, 0, 0, 0, 0, 0],
  [4, 12, 24, 24, 0, 0, 0, 0, 0]
  [5, 20, 60, 120, 120, 0, 0, 0, 0],
  [6, 30, 120, 360, 720, 720, 0, 0, 0],
  [7, 42, 210, 840, 2520, 5040, 5040, 0, 0],
  [8, 56, 336, 1680, 6720, 20160, 40320, 40320, 0],
  [9, 72, 504, 3024, 15120, 60480, 181440, 362880, 362880],
  [10, 90, 720, 5040, 30240, 151200, 604800, 1814400, 3628800],
  [11, 110, 990, 7920, 55440, 332640, 1663200, 6652800, 19958400],
  [12, 132, 1320, 11880, 95040, 665280, 3991680, 19958400, 79833600],
  [13, 156, 1716, 17160, 154440, 1235520, 8648640, 51891840, 259459200],
  [14, 182, 2184, 24024, 240240, 2162160, 17297280, 121080960, 726485760],
  [15, 210, 2730, 32760, 360360, 3603600, 32432400, 259459200, 1816214400],
  [16, 240, 3360, 43680, 524160, 5765760, 57657600, 518918400, 4151347200],
  [17, 272, 4080, 57120, 742560, 8910720, 98017920, 980179200, 8821612800],
  [18, 306, 4896, 73440, 1028160, 13366080, 160392960, 1764322560, 17643225600],
  [19, 342, 5814, 93024, 1395360, 19535040, 253955520, 3047466240, 33522128640]]

OK, so we get something like this:

def falling_factorial(x, n):
  accumulator = 1
  while n > 0:
    accumulator *= x
    x -= 1
    n -= 1
  return accumulator

What next? Well, we could inline the two refactored functions into the trailing_zeroes_of_falling_factorial function.

First we get some clamp tests around it:

def exp_of_b_in(b, n):
  accumulator = 0
  while n >= b and n % b == 0:
    accumulator += 1
    n = n / b
  return accumulator

def falling_factorial(x, n):
  accumulator = 1
  while n > 0:
    accumulator *= x
    x -= 1
    n -= 1
  return accumulator

def trailing_zeroes_of_falling_factorial(x, n):
  return exp_of_b_in(10, falling_factorial(x, n))

print [[trailing_zeroes_of_falling_factorial(x, n) for n in range(1, 20)] for x in range(1, 20)]
assert [[trailing_zeroes_of_falling_factorial(x, n) for n in range(1, 20)] for x in range(1, 20)] == [
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0],
  [0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0, 0, 0, 0],
  [0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0, 0, 0],
  [0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0, 0],
  [0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 0],
  [0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3]]

Then we might get something like this (renaming apart the two accumulators):

def trailing_zeroes_of_falling_factorial(x, n):
  falling_factorial_accumulator = 1
  while n > 0:
    falling_factorial_accumulator *= x
    x -= 1
    n -= 1
  exp_of_b_in_accumulator = 0
  while falling_factorial_accumulator >= 10 and falling_factorial_accumulator % 10 == 0:
    exp_of_b_in_accumulator += 1
    falling_factorial_accumulator /= 10
  return exp_of_b_in_accumulator

So now I want to fuse the two loops. I know that loop fusion is a thing, but I'm not very practiced at it, I don't know all the rules, and so I will just have to muddle through (with clamp tests to make sure I don't mess anything up).

My idea was would like to do something like this:

def trailing_zeroes_of_falling_factorial(x, n):
  falling_factorial_accumulator = 1
  exp_of_b_in_accumulator = 0
  while n > 0:
    falling_factorial_accumulator *= x
    x -= 1
    n -= 1
    while falling_factorial_accumulator > 0 and falling_factorial_accumulator % 10 == 0:
      exp_of_b_in_accumulator += 1
      falling_factorial_accumulator /= 10
  return exp_of_b_in_accumulator

My reasoning was that if we know we're going to divide by ten out eventually, and we want to avoid big intermediate results, maybe we could divide by ten immediately? But my clamp tests alerted me that I had messed something up. After a bit of debugging, I figured out that there is a case falling_factorial_accumulator doesn't just keep getting bigger and bigger - when we multiply by zero.

So, special casing that possibility out of the way, we end up with this:

def trailing_zeroes_of_falling_factorial(x, n):
  falling_factorial_accumulator = 1
  exp_of_b_in_accumulator = 0
  if n > x:
    return 0
  while n > 0:
    falling_factorial_accumulator *= x
    x -= 1
    n -= 1
    while falling_factorial_accumulator > 0 and falling_factorial_accumulator % 10 == 0:
      exp_of_b_in_accumulator += 1
      falling_factorial_accumulator /= 10
  return exp_of_b_in_accumulator

That while loop looks annoying. Do we really need that first part of the test?

def trailing_zeroes_of_falling_factorial(x, n):
  falling_factorial_accumulator = 1
  exp_of_b_in_accumulator = 0
  if n > x:
    return 0
  while n > 0:
    falling_factorial_accumulator *= x
    x -= 1
    n -= 1
    while falling_factorial_accumulator % 10 == 0:
      exp_of_b_in_accumulator += 1
      falling_factorial_accumulator /= 10
  return exp_of_b_in_accumulator

Tried it, apparently not. What about moving those x and n changes away, so that we can get the multiplication closer to the modulus?

def trailing_zeroes_of_falling_factorial(x, n):
  falling_factorial_accumulator = 1
  exp_of_b_in_accumulator = 0
  if n > x:
    return 0
  while n > 0:
    falling_factorial_accumulator *= x
    while falling_factorial_accumulator % 10 == 0:
      exp_of_b_in_accumulator += 1
      falling_factorial_accumulator /= 10
    x -= 1
    n -= 1
  return exp_of_b_in_accumulator

What if we renamed this long falling_factorial_accumulator just f? And while we're at it, the other accumulator a?

def trailing_zeroes_of_falling_factorial(x, n):
  f = 1
  a = 0
  if n > x:
    return 0
  while n > 0:
    f *= x
    while f % 10 == 0:
      a += 1
      f /= 10
    x -= 1
    n -= 1
  return a

What can we say about multiplication followed by testing whether something divides 10? Intuitively, the output is round if the input (of the multiplication f * x) is round. For example, 100 * 100 (round times round) is 10000 (very round).

This situation reminds me of an idea I've heard of: "all numbers can be expressed as products of prime powers", or "the fundamental theorem of arithmetic". So f can be expressed as 2^x*3^y*5^z*... and x can too. And then when we multiply them together, the answer will be zero mod 10 if x (the exponent of 2 in f * x), is greater than zero AND z (the exponent of 5 in f * x) is greater than zero.

So if we just tracked the exponent of 2 in f, and the exponent of 5 in f, then we could sortof directly track the roundness of f, rather than the value of f, which would cut down on the size of those intermediate numbers.

def trailing_zeroes_of_falling_factorial(x, n):
  f2 = 0
  f5 = 0
  a = 0
  if n > x:
    return 0
  while n > 0:
    f2 += exp_of_b_in(2, x)
    f5 += exp_of_b_in(5, x)
    while f2 > 0 and f5 > 0:
      a += 1
      f2 -= 1
      f5 -= 1
    x -= 1
    n -= 1
  return a

Now that loop, counting out 2s and 5s into the accumulator a looks a little silly. Why go step by step, when we know that the number of times round the loop will be whichever is smaller, f2 or f5?

def trailing_zeroes_of_falling_factorial(x, n):
  if n > x:
    return 0
  f2 = 0
  f5 = 0
  while n > 0:
    f2 += exp_of_b_in(2, x)
    f5 += exp_of_b_in(5, x)
    x -= 1
    n -= 1
  m = min(f2, f5)
  a = m
  f2 -= m
  f5 -= m
  return a

We don't really need to compute those last updates to a, f2, and f5. Eliminating dead code is such a simple refactoring, that it often doesn't get mentioned as a refactoring.

def trailing_zeroes_of_falling_factorial(x, n):
  if n > x:
    return 0
  f2 = 0
  f5 = 0
  while n > 0:
    f2 += exp_of_b_in(2, x)
    f5 += exp_of_b_in(5, x)
    x -= 1
    n -= 1
  return min(f2, f5)

There's some Python syntactic sugar that might express this more tersely.

def trailing_zeroes_of_falling_factorial(x, n):
  if n > x:
    return 0
  f2 = sum(exp_of_b_in(2, x) for x in range(x, x - n, -1))
  f5 = sum(exp_of_b_in(5, x) for x in range(x, x - n, -1))
  return min(f2, f5)

That's a bit elegant, but maybe we could save a few characters by flipping the direction of the range? This is not exactly because we care about saving characters, and more that we want to make changes, and sometimes we use the excuse of saving characters to motivate making changes.

def trailing_zeroes_of_falling_factorial(x, n):
  if n > x:
    return 0
  f2 = sum(exp_of_b_in(2, x) for x in range(x - n + 1, x + 1))
  f5 = sum(exp_of_b_in(5, x) for x in range(x - n + 1, x + 1))
  return min(f2, f5)

Getting the bounds on those ranges was annoying. Maybe using the python list comprehension wasn't such a good idea.

def exp_of_b_in(b, n):
  accumulator = 0
  while n >= b and n % b == 0:
    accumulator += 1
    n = n / b
  return accumulator

def trailing_zeroes_of_falling_factorial(x, n):
  if n > x:
    return 0
  f2 = 0
  for y in range(x - n + 1, x + 1):
    f2 += exp_of_b_in(2, y)
  f5 = 0
  for z in range(x - n + 1, x + 1):
    f5 += exp_of_b_in(5, z)
  return min(f2, f5)

Should we inline these two calls to the helper function?

def trailing_zeroes_of_falling_factorial(x, n):
  if n > x:
    return 0
  f2 = 0
  for y in range(x - n + 1, x + 1):
    while y >= 2 and y % 2 == 0:
      f2 += 1
      y /= 2
  f5 = 0
  for z in range(x - n + 1, x + 1):
    while z >= 5 and z % 5 == 0:
      f5 += 1
      z /= 5
  return min(f2, f5)

Huh. I think this is might be "going too far". Let's go back a couple refactorings, and call it good?

def exp_of_b_in(b, n):
  accumulator = 0
  while n >= b and n % b == 0:
    accumulator += 1
    n = n / b
  return accumulator

def trailing_zeroes_of_falling_factorial(x, n):
  if n > x:
    return 0
  f2 = sum(exp_of_b_in(2, x) for x in range(x - n + 1, x + 1))
  f5 = sum(exp_of_b_in(5, x) for x in range(x - n + 1, x + 1))
  return min(f2, f5)