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

в 9:06, , рубрики: c++, Алгоритмы, Биномиальные коэффициенты, математика, метки:

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

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

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

Алгоритм

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

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

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;
}

Результат выполнения:
C(34, 17) = 2333606220
C(35, 17) = 4537567650
C(36, 18) = 9075135300
C(37, 18) = 17672631900
C(38, 19) = 35345263800
C(39, 19) = 68923264410
C(40, 20) = 137846528820
C(41, 20) = 269128937220
C(42, 21) = 538257874440
C(43, 21) = 1052049481860
C(44, 22) = 2104098963720
C(45, 22) = 4116715363800
C(46, 23) = 8233430727600
C(47, 23) = 16123801841550
C(48, 24) = 32247603683100
C(49, 24) = 63205303218876
C(50, 25) = 126410606437752
C(51, 25) = 247959266474052
C(52, 26) = 495918532948104
C(53, 26) = 973469712824056
C(54, 27) = 1946939425648112
C(55, 27) = 3824345300380220
C(56, 28) = 7648690600760440
C(57, 28) = 15033633249770520
C(58, 29) = 30067266499541040
C(59, 29) = 59132290782430712
C(60, 30) = 118264581564861424
C(61, 30) = 232714176627630544
C(62, 31) = 465428353255261088
C(63, 31) = 916312070471295267
C(64, 32) = 1832624140942590534
C(65, 32) = 3609714217008132870
C(66, 33) = 7219428434016265740
C(67, 33) = 14226520737620288370

Process returned 0 (0x0) execution time: 0.051 s

Автор: cranium256

Источник

* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js