Вычисление биномиальных коэффициентов… всё-таки программно

Ранее, в трёх статьях была затронута тема вычисления биномиальных коэффициентов.

Расчет биномиальных коэффициентов на Си (С++)
Расчет биномиальных коэффициентов с использованием Фурье-преобразований
Вычисление биномиальных коэффициентов… вручную

В последней статье автор продемонстрировал чисто математическую оптимизацию вычисления биномиальных коэффициентов. Оказалось, что для их вычисления и компьютер-то не особенно нужен. Хватит бумаги, ручки, калькулятора и пива. Однако, если одной из вышеперечисленных вещей под рукой не окажется, но в области досягаемости будет простаивающий компьютер, то можно сделать то же самое и на компьютере.

Алгоритм

В комментариях к статье автору попеняли за отсутствие формализованного алгоритма. После недолгих размышлений нарисовалось следующее.

При вычислении биномиального коэффициента после приведения факториальной формулы к виду:

C(n, k) = П(n, k+1) / П(1, n-k)

Получаем дробь, которую необходимо сократить. Как проще всего сократить дробь в целых числах? Разложить числитель и знаменатель на простые множители. Дело облегчается тем, что в данном случае, и числитель и знаменатель являются произведениями небольших чисел. Т.е. на простые множители можно раскладывать каждый из сомножителей по отдельности, накапливая результат в виде списка простых чисел с соответствующими показателями степени.

Затем, проходим по списку, полученному для знаменателя, берём каждое число и из показателя степени для этого числа в числителе вычитаем показатель степени для этого числа в знаменателе. Исходя из факта, что биномиальный коэффициент — число целое, знаменатель должен сократиться полностью, т.е стать равным 1.

Последний шаг: проходим по списку, полученному для числителя на предыдущем шаге, и перемножаем все числа, возведённые в степень, соответствующую этому числу.

Реализация

Хранить пары «простое число — показатель степени» было решено в контейнере map, где простое число является ключом, а показатель степени — значением.

Для факторизации чисел был взят самый примитивный алгоритм. Для оптимизации дальнейшего использования, в функцию факторизации по ссылке передаётся внешний контейнер.

#include <iostream>
#include <map>

using namespace std;

typedef unsigned long long ullong;

// Функция разложения числа на простые множители.
// В качестве аргументов получает число и
// ссылку на контейнер map, в который записывает разложение в виде:
// ключ - простое число, значение - степень этого числа в разложении.

template <typename integral_type>
void factorize(integral_type num, map<integral_type, unsigned int>& res) {

	integral_type j = 2;
	while (num / integral_type(2) >= j) {
		if (num % j == 0) {
			res[j]++;
			num /= j;
			j = integral_type(2);
		}
		else {
			++j;
		}
	}
	res[num]++;
}

typedef map<ullong, unsigned int> factmap;


// Функция вычисления биномиального коэффициента.

ullong C(int n, int k) {

	ullong result = 1uLL;

	// Факториальная формула биномиального коэффициента приводится к виду
	// C(n, k) = П(n, k+1) / П(1, n-k)

	int lim_numer_min = k + 1;  // минимальное число в произведении для числителя
	int lim_numer_max = n;      // максимальное число в произведении для числителя
	int lim_denom_min = 2;      // минимальное число в произведении для знаменателя
	int lim_denom_max = n - k;  // максимальное число в произведении для знаменателя

	// контейнеры для разложения на простые множители числителя и знаменателя
	factmap numerator, denominator;

	// раскладываем все сомножители числителя на простые множители
	for (int i = lim_numer_min; i <= lim_numer_max; ++i) {
		factorize(ullong(i), numerator);
	}
	// раскладываем все сомножители знаменателя на простые множители
	for (int i = lim_denom_min; i <= lim_denom_max; ++i) {
		factorize(ullong(i), denominator);
	}

	// производим сокращение числителя на знаменатель
	for (auto i = denominator.begin(); i != denominator.end(); i++) {
		numerator[i->first] -= i->second;
	}

	// вычисляем значение биномиального коэффициента по
	// несокращённым простым множителям числителя
	for (auto i = numerator.begin(); i != numerator.end(); i++) {
		// возведение в степень в целых числах умножением
		for (ullong p = 0; p < i->second; ++p) {
			result *= i->first;
		}
	}

	return result;
}

Тестирование

Головная программа была сделана наподобие программы из первой статьи:

int main() {
	int n;
	for (n = 34; n < 68; ++n) {
		cout << "C(" << n << ", " << n/2 << ") = " << C(n, n/2) << endl;
	}
	return 0;
}