В некоторых случаях необходимо считать по некоторому простому модулю p сложные формулы, которые в том числе могут содержать факториалы. Здесь мы рассмотрим случай, когда модуль p сравнительно мал. Понятно, что эта задача имеет смысл только в том случае, когда факториалы входят и в числитель, и в знаменатель дробей. Действительно, факториал p! и все последующие обращаются в ноль по модулю p, однако в дробях все множители, содержащие p, могут сократиться, и полученное выражение уже будет отлично от нуля по модулю p.

Таким образом, формально задача такая. Требуется вычислить n! по простому модулю p, при этом не учитывая все кратные p множители, входящие в факториал. Научившись эффективно вычислять такой факториал, мы сможем быстро вычислять значение различных комбинаторных формул (например, Биномиальные коэффициенты).

Выпишем этот "модифицированный" факториал в явном виде:

n!{%p} = 1 * 2 * 3 * ... * (p - 2) * (p - 1) * {1}p * (p + 1) * (p + 2) * ... * (2p - 1) * {2}2p * (2p + 1) * ... * (p^2 - 1) * {1}p^2 * (p^2 + 1) * ... * n = 1 * 2 * 3 * ... * (p - 2) * (p - 1) * {1}p * 1 * 2 * ... * (p - 1) * {2}2p * 1 * 2 * ... * (p - 1) * {1}p^2 * 1 * 2 * ... * (n % p) * (mod p)

При такой записи видно, что "модифицированный" факториал распадается на несколько блоков длины p, которые все одинаковы, за исключением последнего элемента:
n!{%p} = {1 * 2 * ... * (p - 2) * (p - 1) * 1} * {1 * 2 * (p - 1) * 2} * ... * {1 * 2 * (p - 1) * 1} * ... * {1 * 2 * ... * n%p } (mod p)

Общую часть блоков посчитать легко — это просто (p-1)! mod p, которую можно посчитать программно или по теореме Вильсона (Wilson) сразу найти (p-1)! mod p = p - 1. Чтобы перемножить эти общие части всех блоков, надо найденную величину возвести в степень по модулю p, что можно сделать за O(log n) операций (см. статью Бинарное_возведение_в_степень); впрочем, можно заметить, что мы фактически возводим минус единицу в какую-то степень, а потому результатом всегда будет либо 1, либо p-1, в зависимости от чётности показателя. Значение в последнем, неполном блоке тоже можно посчитать отдельно за O(p). Остались только последние элементы блоков, рассмотрим их внимательнее:

n%p = {... * 1} * {... * 2} * {... * 3} * ... * {... * (p - 1)} * {... * 2} * {... * 1} * {... * 1}

И мы снова пришли к "модифицированному" факториалу, но уже меньшей размерности (столько, сколько было полных блоков, а их было [ n / p ]. Таким образом, вычисление "модифицированного" факториала n!{%p} мы свели за O(p) операций к вычислению уже (n/p)!{%p}. Раскрывая эту рекуррентную зависимость, мы получаем, что глубина рекурсии будет O (log_p(n)), итого асимптотика алгоритма получается O(p log_p(n)).

Понятно, что при реализации не обязательно использовать рекурсию в явном виде: поскольку рекурсия хвостовая, её легко развернуть в цикл:

int factmod (int n, int p) {
	int res = 1;
	while (n > 1) {
		res = (res * ((n/p) % 2 ? p-1 : 1)) % p;
		for (int i=2; i<=n%p; ++i)
			res = (res * i) % p;
		n /= p;
	}
	return res % p;
}

Эта реализация работает за O(p log_p n).