Algorithm to find smallest N such that N! is divisible by a prime raised to a power

Question

Is there an efficient algorithm to compute the smallest integer N such that N! is divisible by p^k where p is a relatively small prime number and k, a very large integer. In other words,

factorial(N) mod p^k == 0

If, given N and p, I wanted to find how many times p divides into N!, I would use the well-known formula

k = Sum(floor(N/p^i) for i=1,2,...

I've done brute force searches for small values of k but that approach breaks down very quickly as k increases and there doesn't appear to be a pattern that I can extrapolate to larger values.

Edited 6/13/2011

Using suggestions proposed by Fiver and Hammar, I used a quasi-binary search to solve the problem but not quite in the manner they suggested. Using a truncated version of the second formula above, I computed an upper bound on N as the product of k and p (using just the first term). I used 1 as the lower bound. Using the classic binary search algorithm, I computed the midpoint between these two values and calculated what k would be using this midpoint value as N in the second formula, this time with all the terms being used.

If the computed k was too small, I adjusted the lower bound and repeated. Too big, I first tested to see if k computed at midpoint-1 was smaller than the desired k. If so, midpoint was returned as the closest N. Otherwise, I adjusted the highpoint and repeated.

If the computed k were equal, I tested whether the value at midpoint-1 was equal to the value at midpoint. If so, I adjusted the highpoint to be the midpoint and repeated. If midpoint-1 was less than the desired k, the midpoint was returned as the desired answer.

Even with very large values for k (10 or more digits), this approach works O(n log(n)) speeds.

Answer 1

OK this is kind of fun.

Define f(i) = (p^i - 1) / (p - 1)

Write k in a kind of funny "base" where the value of position i is this f(i).

You do this from most-significant to least-significant digit. So first, find the largest j such that f(j) <= k. Then compute the quotient and remainder of k / f(j). Store the quotient as q_j and the remainder as r. Now compute the quotient and remainder of r / f(j-1). Store the quotient as q_{j-1} and the remainder as r again. Now compute the quotient and remainder of r / f(j-2). And so on.

This generates a sequence q_j, q_{j-1}, q_{j-2}, ..., q_1. (Note that the sequence ends at 1, not 0.) Then compute q_j*p^j + q_{j-1}*p^(j-1) + ... q_1*p. That's your N.

Example: k = 9, p = 3. So f(i) = (3^i - 1) / 2. f(1) = 1, f(2) = 4, f(3) = 13. So the largest j with f(j) <= 9 is i = 2 with f(2) = 4. Take the quotient and remainder of 9 / 4. That's a quotient of 2 (which is the digit in our 2's place) and remainder of 1.

For that remainder of 1, find the quotient and remainder of 1 / f(1). Quotient is 1, remainder is zero, so we are done.

So q_2 = 2, q_1 = 1. 2*3^2 + 1*3^1 = 21, which is the right N.

I have an explanation on paper for why this works, but I am not sure how to communicate it in text... Note that f(i) answers the question, "how many factors of p are there in (p^i)!". Once you find the largest i,j such that j*f(i) is less than k, and realize what you are really doing is finding the largest j*p^i less than N, the rest kind of falls out of the wash. In our p=3 example, for instance, we get 4 p's contributed by the product of 1-9, 4 more contributed by the product of 10-18, and one more contributed by 21. Those first two are just multiples of p^2; f(2) = 4 is telling us that each multiple of p^2 contributes 4 more p's to the product.

[update]

Code always helps to clarify. Save the following perl script as foo.pl and run it as foo.pl <p> <k>. Note that ** is Perl's exponentiation operator, bdiv computes a quotient and remainder for BigInts (unlimited-precision integers), and use bigint tells Perl to use BigInts everywhere.

#!/usr/bin/env perl

use warnings;
use strict;
use bigint;

@ARGV == 2
    or die "Usage: $0 <p> <k>\n";

my ($p, $k) = map { Math::BigInt->new($_) } @ARGV;

sub f {
    my $i = shift;
    return ($p ** $i - 1) / ($p - 1);
}

my $j = 0;
while (f($j) <= $k) {
    $j++;
}
$j--;

my $N = 0;

my $r = $k;
while ($r > 0) {
    my $val = f($j);
    my ($q, $new_r) = $r->bdiv($val);
    $N += $q * ($p ** $j);
    $r = $new_r;
    $j--;
}

print "Result: $N\n";

exit 0;

Answer 2

Using the formula you mentioned, the sequence of k values given fixed p and N = 1,2... is non-decreasing. This means you can use a variant of binary search to find N given the desired k.

• Start with N = 1, and calculate k.
• Double N until k is greater or equal than your desired k to get an upper bound.
• Do a binary search on the remaining interval to find your k.