Расчет биномиальных коэффициентов на Си (С++)
Расчет биномиальных коэффициентов с использованием Фурье-преобразований
Вычисление биномиальных коэффициентов… вручную
В последней статье автор продемонстрировал чисто математическую оптимизацию вычисления биномиальных коэффициентов. Оказалось, что для их вычисления и компьютер-то не особенно нужен. Хватит бумаги, ручки, калькулятора
Алгоритм
В комментариях к статье автору попеняли за отсутствие формализованного алгоритма. После недолгих размышлений нарисовалось следующее.
При вычислении биномиального коэффициента после приведения факториальной формулы к виду:
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;
}
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
Комментарии (10)
PavelSandovin
12.01.2016 18:35-1Дело облегчается тем, что в данном случае, и числитель и знаменатель являются произведениями небольших чисел
Хорошо, давайте посчитаем этим методом, например,
С(394959569696694800489675432,2340540340005932332490003325025398054778666543).
Просто любопытно, сколько времени займет факторизация всех входящих в числитель и знаменатель множителей?
a_batyr
12.01.2016 18:51+1С(394959569696694800489675432,2340540340005932332490003325025398054778666543).
У вас k > n, как вы собираетесь считать факториал отрицательного числа?PavelSandovin
12.01.2016 18:59+1В данном случае — никак. По определению, С(n, k) = 0 при k > n.
И что-то я такую проверку в коде в статье не увидел :)a_batyr
12.01.2016 19:23+1Да, действительно, не проверяет :)
Но и длинной арифметики он не реализовывал, тогда достаточно запустить C(2147483647, 1073741824)
Меня вот больше беспокоит, что в функции main() переменная k не используется. А ещё при n==k неправильно тоже.
cranium256
12.01.2016 21:01Да, корректность аргументов функция не проверяет. И длинную арифметику не реализовывал. И для факторизации был взят самый примитивный алгоритм. Цель другая была: формализовать алгоритм, описанный в статье V2008n.
Тем не менее, для целой 64-разрядной арифметики всё работает быстро, точно и без переполнения (естественно, в пределах разумного).
Кстати, при n == k работает правильно. Циклы факторизации числителя и знаменателя не выполняются, поскольку условия продолжения циклов ложны. В итоге функция возвращает 1, как и ожидается.
Переменная k не используется — да, грешен. Проглядел. Осталась от предыдущей версии функции main(). Вся надежда на оптимизирующий компилятор ;)
Если же рассматривать вопрос вычисления биномиальных коэффициентов с применением длинной арифметики, то проблема переполнения, с которой, собственно, и началась эта серия статей, отпадает сразу. Останутся традиционные вопросы: скорость работы и занимаемая память. Но это как бы уже далеко за рамками данной статьи.
PS. Переменную k из кода убрал.a_batyr
12.01.2016 21:19На самом деле вот этот комментарий уже раскрывал вопрос разложения на множители.
Ваш код тоже можно было б оставить в виде комментария под спойлером, не будь вы read-only.
Добро пожаловать.
PavelSandovin
13.01.2016 11:08Попробуйте реализовать еще одну версию, вычисляющую биномиальный коээффициент «в лоб», без факторизации и сравнить производительность обоих версий.
Делаю ставку на то, что даже при вычислениях в пределах диапазона типа INT без факторизации будет быстрее.
А за его пределами, как в моем примере — только если поменять местами аргументы — факторизацией пользоваться вообще нереально.
dcc0
12.01.2016 22:04Что означает «самый примитивный»? Т.е. какой именно алгоритм?
Метод Ферма тоже довольно прозрачен, и операций в нем немного, я даже на php с ним игрался.
Код C# понимаю очень слабо, поэтому спросил.
Дело в том, что у Вас реализация алгоритма, формализация на языке, конкретном, а хотелось бы
на естественном… или по шагам.cranium256
13.01.2016 03:27Наверное этот алгоритм называется «перебор делителей». Хотя я над этим не задумывался.
«Самый примитивный» — потому что вообще без оптимизаций.
Вообще-то код у меня на кондовом C++ ;)
Примерно так будет выглядеть эта функция на PHP (сохранение простых делителей заменено просто на вывод):
function factorize($num) { $j = 2; while ($num / 2 >= $j) { if ($num % $j == 0) { echo $j; $num /= $j; $j = 2; } else { ++$j; } } echo $num; }
Переменная j — это проверяемый делитель. В самом плохом случае, если num постое, в цикле while она пробежит все значения от 2 до num / 2 (условие продолжения цикла). Если же num — не простое, то при каком-то значении j остаток от деления станет равным 0, следовательно текущее значение является делителем num. В это случае выводим найденное значение, num делим на это значение и j сбрасываем на начальное значение (2). Таким образом находим все делители, кроме последнего, в порядке возрастания. Последний делитель — значение num после завершения цикла.
a_batyr
Астрологи объявили неделю биномиальных коэффициентов. Количество математиков выросло.
Вангую вычисление БК на ардуино за 5 долларов.