Computing n! modulo pᵏ for small p

Hi everyone!

There is an article on cp-algorithms about how to compute modulo prime number . Today I was wondering, how to do it if we need to compute modulo rather than just . Turns out one can do it efficiently if is small.

Motivational example

Xiaoxu Guo Contest 3 — Binomial Coefficient. Given , find modulo .

There are also some projecteuler problems that require it, including with several queries of distinct and .

Task formulation and outline of result

To clarify the task a bit, our ultimate goal here is to be able to compute e.g. binomial coefficients modulo . Thus, what we really need is to represent , where , and then report and modulo . This is sufficient to also compute modulo .

We will show that, assuming that polynomial multiplication of size requires operations, and assuming that arithmetic operations modulo take , we can find and in , where . It requires pre-computation that takes memory.

Algorithm with polynomials

Great thanks to Endagorion for sharing this algorithm with me!

Let , then we can represent as

So, there is a part contributed by numbers divisible by , which can be reduced to computing , and a part contributed by numbers not divisible by . Let , then the later can be represented as follows:

To compute it quickly, we can define a family of polynomial for such that

so the value of the factorial would be represented as

and expanding it then rewrites into

which simplifies as

Now, what would it take us to use this setup? First of all, notice that can be computed from one another:

Note that for shorter and more consistent implementation, this recurrent formula also mostly works for if we put , but for we should go up to instead of . We should also note that for , we only care for coefficients up to , as the larger ones are divisible by .

This allows to compute in for all for a given . Over all from to it sums up to time and memory. Then, evaluating all the polynomials for any specific requires operations, where .

As it requires some manipulations with polynomials, I implemented it in Python with sympy just as a proof of concept:

Algorithm with -adic logarithms and exponents

The algorithm above requires some heavy machinery on polynomials. I also found out an alternative approach that allows to compute and for any given divisible by in with precomputation, where . It doesn't rely on polynomial operations, except for Lagrange interpolation.

Let , where , then we can represent its factorial as

We can further rewrite it as

Let's learn how to compute the product

as with it we can represent the factorial as

We can take a -adic logarithm (please refer to this article for definitions and properties) from the product to get

We can precompute sums of for each up to and up to in , and we can find the sum of for from to in via Lagrange interpolation. Therefore, with precomputation, we can compute in , which then allows to find in , where . I implemented it in Python, as a proof of concept, but I didn't bother with faster sums over yet:

Note: Due to inherent properties of -adic logarithms and exponents, the algorithm only works with , and it might return instead with . This is because with when has remainder modulo . It is accounted for in the code by tracking the number of integers below that have remainder modulo .

Some other approaches

There is also a blog by prabowo, based apparently on another blog by Min_25, and a scholarly article by Andrew Granville. I'm not sure whether they talk about the same algorithms as described here.