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

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

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

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

Алгоритм

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

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

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, k; 	for (n = 34; n < 68; ++n) { 		cout << "C(" << n << ", " << n/2 << ") = " << C(n, n/2) << endl; 	} 	return 0; } 

ссылка на оригинал статьи http://habrahabr.ru/post/274911/