Как укротить биномиальный коэффициент?

Изменено: . Разместил: FadeDemon | Меню
Количество просмотров2`695 / FD1159-6-0

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


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


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


\[ \binom{n}{k} = {n! \over {k! \cdot (n - k)!}}, \, \, \, k \leqslant n \]
Определение биномиального коэффициента

Как нетрудно подметить, в формуле фигурирует факториал n, а факториал, как мы помним, функция неприлично быстро растущая, и при попытке посчитать тот самый факториал 13-ти мы заметим, что полученное число  6`227`020`800 уже не умещается в 4-х байтовое число типа long, превосходя его максимальное значение примерно на 4 миллиарда. Для unsigned long будет впрочем уже ≈2 миллиарда, но я думаю легче нам от этого не станет, поскольку факториал 14-ти представляет из себя уже примерно 87 миллиардов, факториал 15 ≈ 1307 млрд., и так далее, что на первый взгляд, неиллюзорно заставляет нас уронить челюсть, и впасть в глубочайший траур… Стоп. У нас же есть еще 64-битные числа, носящие имена Int64, или long long. Конечно, небольшую отдушинку они дают. До тех самых пор, пока факториал n не станет больше чем 263, то есть, приблизительно до n = 21. И всё, дальше снова переполнение. Признаться, когда я понял это впервые, данный факт несколько опечалил меня. При таком прямолинейном подходе получалось, что нельзя посчитать биномиальный коэффициент для множества, содержащего более 20 элементов. Но печаль быстро рассеялась, и в очередном творческом порыве, решение сией проблемы быстро было найдено.


Давайте мы с вами вспомним, как собственно говоря определяется функция факториала:


\[ n! = \prod_{i=1}^{n} i \]
Факториал

Тобишь, это просто напросто произведение всех целых чисел от 1 до n, и ничего более. Впрочем все мы это знаем, и ничего необычного в этом нет. Другое дело, что не сразу можно подметить одну особенность, что в факториалах любых двух чисел всегда (ну, или почти всегда) есть одинаковые сомножители, что вытекает из собственно говоря самого определения факториала. Давайте приведем простой, канонiческий я бы сказал пример, вычисления биномиального коэффициента для n = 6 и k = 3:


\[ {n! \over {k! \cdot (n - k)!}} = {1 \cdot 2 \cdot 3 \cdot 4 \cdot 5 \cdot 6 \over 1 \cdot 2 \cdot 3 \cdot (6 - 3)!} = {1 \cdot 2 \cdot 3 \cdot 4 \cdot 5 \cdot 6 \over \left (1 \cdot 2 \cdot 3\right ) \cdot \left (1 \cdot 2 \cdot 3\right )} = {720 \over 6 \cdot 6} = 20 \]

Конечно, считая в лоб мы получаем 20, впрочем ничего другого мы не могли получить. Но, давайте взглянем на 3-ю дробь, ведь совсем не трудно заметить, что там у нас минимум 3 раза мелькает последовательность 1 ⋅ 2 ⋅ 3, которую ведь можно просто-напросто взять, и сократить:


\[ {n! \over {k! \cdot (n - k)!}} = {\cancel 1 \cdot \cancel 2 \cdot \cancel 3 \cdot 4 \cdot 5 \cdot 6 \over \left (\cancel 1 \cdot \cancel 2 \cdot \cancel 3\right ) \cdot \overstrike{\left (1 \cdot 2 \cdot 3\right )}} = {4 \cdot 5 \cdot 6 \over 1 \cdot 2 \cdot 3} = {120 \over 6} = 20 \]


И снова получилось 20 (хотя вряд-ли стоило ждать чего-то другого), но на этот раз мы оперировали числами, в 6 раз меньшими, чем в предыдущем примере. Занимательно, не правда-ли? Хотя на самом деле всё это очевидно, и мы, кстати говоря, уже почти решили проблему переполнения из-за факториалов, по крайней мере для k значительно меньших чем n (пример, n = 100, k = 5), поскольку сокращается в этом случае то, что стоит в скобках, и сокращается ровно на 95 множителей. Такой подход, собственно говоря, и позволил мне решить ту самую задачу. Но, однако, поразмыслив ещё, мы можем прийти к ещё более малым числам, практически полностью решив проблему переполнения для любых n и k (конечно же, при условии, что само значение биномиального коэффициента не больше чем 231 или 263):


\[ {n! \over {k! \cdot (n - k)!}} = {4 \cdot 5 \cdot 6 \over 1 \cdot 2 \cdot 3} = {4 \over 1} \cdot {5 \over 2} \cdot {6 \over 3} = 4 \cdot 2.5 \cdot 2 = 20 \]


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


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


\[ \binom{n}{k} = {\prod\limits_{i=1}^{n}i \over {\prod\limits_{i=1}^{k}i \cdot \prod\limits_{i=1}^{n - k}i}} \]
Если взглянуть глубже

И теперь мы разберем простую ситуацию, когда k < (n − k). Очевидно в этом случае, что нижняя правая Π в формуле сокращается, точно также, как мы это делали на простых примерах, и в следствии этого, у верхнего произведения меняется стартовое значение i, причем так, что оно становится равным (n − k) + 1, в чем нетрудно убедиться на простых примерах выше, поскольку вверху там у нас была четверочка, а k равнялся троечке. Именно нижнее правое произведение мы сокращаем по очевидной причине, что раз k < (n − k), значит (n − k) > k, а значит и его факториал будет больше, и соответственно сократить мы сможем больше множителей, чем если бы сокращали левое произведение с пределом в k. Запишем теперь, что у нас получается:


\[ {\prod\limits_{i=1}^{n}i \over {\prod\limits_{i=1}^{k}i \cdot \prod\limits_{i=1}^{n - k}i}} = {\prod\limits_{i = n - k + 1}^{n}i \over {\prod\limits_{i=1}^{k}i}} \]
Преобразование первое

Теперь наша запись представляет собой более-менее строгое описания ситуации в простых примерах, когда в числителе у нас получилось 120, оптимальной для k < (n − k). Но, мы веть смогли придумать и более оптимальный алгоритм, верно? Так вот, давайте мы и его теперь попробуем строго записать, отталкиваясь от того, что получилось у нас выше. Для этого, сначала, приведем оба произведения в числителе и знаменателе к одинаковым пределам, для того, чтобы потом можно было перейти к одному произведению с дробью. Предел нижнего Π у нас неизменен, и все еще равен 1 … k, вот и попробуем в верхнем Π нарисовать такой же предел, вместо той страхоёбины того выражения, которое у нас получилось, и я уже заранее думаю, что это у нас прекрасно получиться:


\[ \begin{aligned}& \prod\limits_{i = n - k + 1}^{n}i = \prod\limits_{i = \cancel n - k + 1 - \cancel n}^{n - n}\left (i + n\right ) = \prod\limits_{i = - k + 1}^{0}\left (i+n\right ) = \prod\limits_{i = -\cancel k + 1 + \cancel k}^{0 + k}\left (i+n - k\right ) = \\ & = \prod\limits_{i = 1}^{k}\left (i+n - k\right ) \end{aligned} \]


Полагаю, что сии действия всё-таки нуждаются в некотором разъяснении. Наверно это сразу не очевидно, но если от начального и конечного пределов отнять некоторое число, и одновременно ко всем i, которые встречаются под знаком произведения (или суммы) это-же число прибавлять (или наоборот), то смысл формулы от этого совершенно не измениться, она будет считать абсолютно то-же самое, что и до этой операции, поскольку значения i в каждой итерации будут верно корректироваться, не меняясь в итоге. Вернее было бы наоборот сказать, что на самом деле это совершенно очевидно, и, я надеюсь, никто не будет требовать строгого доказательства этого очень простого факта, на котором мы и основывали наше преобразование пределов, и, если всмотреться, то сделали мы две операции: сначала от верхнего и нижнего пределов отняли n, прибавив одновременно его к i, а затем прибавили к пределам k, одновременно отняв его от i под знаком произведения.


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


\[ {\prod\limits_{i = n - k + 1}^{n}i \over {\prod\limits_{i=1}^{k}i}} = {\prod\limits_{i = 1}^{k}\left (i+n - k\right ) \over {\prod\limits_{i=1}^{k}i}} = \prod\limits_{i = 1}^{k}\left ({i+n - k \over i}\right ) = \prod\limits_{i = 1}^{k}\left ({1+{n - k \over i}}\right )\]
Преобразование второе

Тот факт, что знаки произведения с одинаковыми пределами можно вынести за дробь, я думаю также не должен вызывать затруднений, мы этим по сути дела просто меняем порядок операций одного и того-же приоритета. Последнее преобразование направлено на то, что единичку, как мне кажется, все-таки проще прибавлять, чем i, поскольку на ассемблере это в теории может скомпилировать в одну инструкцию inc, хотя в данном случае теория довольно далека от практики. Но это впрочем не важно. Важно тут то, что рассмотренный нами случай оптимален только для ситуации, когда, как мы помним, k меньше, чем (n − k), в противном случае множители факториалов очевидным образом сокращаются не оптимально. Как же быть в этом случае? В общем-то, по аналогии. Нужно просто вывести формулу из исходного представления биномиального коэффициента, только сокращая не правое, а левое произведение, и получиться формула, очень похожая на выведенную выше, но оптимальная для случая k > (n − k). А ещё можно её не выводить, а внимательно всмотреться в существующую, заодно обозначив выражение (n − k) буквой m. А что же мы там можем увидеть? А то, что получается так, что k и m – те самые значения, которые стояли под факториалом в исходном выражении, и я думаю, что мы не поделим случайно на ноль, не создадим черной дыры, не сожмем пространство в сингулярность, и что самое главное, не нарушим равенства, если попытаемся k и m поменять местами:


\[ \binom{n}{k} = \prod\limits_{i = 1}^{k}\left ({1+{m \over i}}\right ) = \prod\limits_{i = 1}^{m}\left ({1+{k \over i}}\right )\]
Главное равенство формул

Итак. Мы, изрядно попотев, подошли к финальной стадии построения математической модели алгоритма. Мы помним, что m = (n − k), что первая формула оптимальна для случая k < m, а вторая формула оптимальна когда наоборот, k > m. Собрав все эти факты воедино, мы, наконец, можем выразиться окончательной формулой, по которой можно оптимально посчитать наш коэффициент в любых ситуациях, и при любой погоде:


\[ \binom{n}{k} = \left\lbrace \begin{aligned} & \prod\limits_{i = 1}^{k}\left ({1+{m \over i}}\right ), & k \leqslant m \\ & \prod\limits_{i = 1}^{m}\left ({1+{k \over i}}\right ), & k > m\end{aligned} \right . \]
Оптимальная формула биномиального коэффициента

Мда… Количество наукообразия в полученной формуле, конечно слегка зашкаливает, или вернее, близится к летальному. Ну да ладно, нам уже не привыкать. Если коротко пояснить, то формула означет, что если k ≤ m, то используем первую формулу, иначе, если k > m – вторую. Нестрогое неравенство в первом случае деликатно обыгрывает ситуацию, когда k = m. Можем себя еще раз похвалить.


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


Код: C++
#include <cstdlib>
#include <iostream>
#include <conio.h>
#include <math.h>

using namespace std;

double binomial(long n, long k){
    //если k больше n, очевидно что выбрать у нас ничего не получиться
    if(n < k || 0 == n){
        return 0;
    }
    
    //если они равны, то способ только один, как не поверни
    if(n == k){
        return 1;
    }    
    
    //считаем m, он будет постоянен
    long m = n - k;
    long i;
    
    //тут будет храниться результат
    double b = 1;
    
    //эти переменные вводятся, чтобы исключить преобразование типов в каждой итерации
    double k1 = k;
    double m1 = m;
    
    if(k <= m){
        //если первый случай, k меньше чем m
        for(i = 1; i <= k; i++){
            //перемножаем наш ряд
            b *= 1 + m1 / i;
        }
    }else{
        //второй случай, противоположный
        for(i = 1; i <= m; i++){
            //перемножаем другой ряд
            b *= 1 + k1 / i;
        }
    }
    
    //возвращаем округлённое значение
    return round(b);
}

int main(int argc, char *argv[])
{
    long n = 1, k = 1;
    double b;
    
    //зацикливаем все это дело
    while(n && k){
    
        cout << "Please input n and k: ";
        scanf("%d %d", &n, &k);
    
        //считаем коэффициент тут же его выводя, без знаков после запятой
        printf("Binominal coefficient for (%d, %d) is %0.0f\n", n, k, binomial(n, k));
    }
    
    return EXIT_SUCCESS;
}

Данный перл, принимает на вход два значения, собственно n и k, и на выходе печатает рассчитанное значение коэффициента для них. Надо отметить, что точность всё-таки ограничена не 64-битными Int64, а точностью double, то есть немного меньше заявленной выше. Но и этого на самом деле нам уже вполне достаточно для повседневного использования. Попробуем теперь это запустить:


Проверяем правильность подсчёта
Рис. №439. Проверяем правильность подсчёта

Мда, кааак тут много цыфер… Как сразу можем заметить, канонiчный пример для (6, 3) неизменно равен 20, что уже сигнализирует о том, что программа вроде-бы работает правильно. Чтобы проверить остальные значения, мы привлечём небезызвестный WolframAlpha, положившись на его вычислительные способности, как вердикт в последней инстанции. Вы, если хотите, можете сделать это самостоятельно, а мы приведем результат проверки для последнего запроса (45, 26):


Вольфрам в своём репертуаре
Рис. №440. Вольфрам в своём репертуаре

Как видим, все в точности совпадает, хотя чего ещё мы могли ожидать?… Для всех остальных приведенных значений n и k результаты программы и вольфрама также сошлись, за исключением одного случая (1310, 7) — там как раз произошел по всей видимости выход за границу точности чисел double. Программа насчитала 1`292`691`869`666`776`800, а вольфрамчик сказал, что будет 1`292`691`869`666`776`920, и я думаю у нас есть больше оснований верить вольфраму... Кстати последний ещё сказал, что это число звучит как 1 квинтиллион, 292 квадриллиона, 691 триллион, 869 миллиардов (биллионов), 666 миллионов, 776 тысяч 920, и имеет (как может отметить Капитан Очевидность) 18 нулей, точнее представляется как ≈1.292 × 1018. Жутковато как-то. Кстати в этом случае n = 1310, и как снова мы можем заметить, если бы мы считали коэффициент в лоб… Да, конец был бы немного предсказуем. В общем мы можем с уверенностью сказать, что получили программу, с помощью которой можно считать биномиальные коэффициенты для любых n и k (удовлетворяющих условию самого определения коэффициента, разумеется), не больше чем ≈1015, невзирая на то самое факториальное переполнение, которое нас так жестко обламывало в самом начале…


Вот она где, математика, на службе программирования, если так можно выразиться… Ещё один эксперимент, ещё одно рассуждение, еще один результат, ещё один опыт. С Новым 2013 Годом всех читающих! :)

Ключевые слова Программизм, математика, эксперимент, биномиальный коэффициент, факториал, переполнение, произведение, предел, WolframAlpha, оптимизация
QR-код
Интересное
The Gauss cannon. Начало истории... Когда-то давно, примерно года 3-4 назад, когда мысли мои еще были захвачены конструированием самодельных воздушек, стреляющих пластмассовыми пульками, а программирование я только только начинал осваивать, довелось мне узнать про существование электромагнитного оружия, а именно – пушки Гаусса.

Читать »»
Случайные фото