Расчет биномиальных коэффициентов на Си (С++) и Python

[stack_trace]

Извините, но C++ тут совсем непричём.

[Teivaz]

Здесь этот тэг можно применить по аналогии со Stack Overflow, которые пишут:

Use this tag for questions about code (to be) compiled with a C++ compiler.

[stack_trace]

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

[zelserg]

Я все примеры компилировал на С++. Заголовком я хотел сказать «можно брать функции расчета, вставлять их в программы на языках С, С++ и все будет работать». На Бейсике или Паскале это не получится. Что Вас не устраивает?

[stack_trace]

Я же написал, что ничего против не имею, почему вы у меня спрашиваете, что меня не устраивает?

[zelserg]

А хочу понять логику высказывания «С++ непричем». Первый отзыв и сразу «непричем»… Если мне нужно найти некий алгоримХ на языке С++, я пишу в поисковике «алгоритмХ С++» и если поисковик возвращает мне алгоритмХ на PHP, то вот этот РНР — «непричем». А если он возвращает алгоритм на С++ или С, то это то, что надо. Я его беру и использую.

[stack_trace]

Да потому что широко распространено заблуждение, что если вы можете что-то написать на C то код на C++ будет выглядеть также. Это заблуждение и это чаще всего неверно. То, что код для С компилируется плюсовым компилятором почти всегда и так всем известно, это не делает его кодом на С++. Этот же код, например, отлично скомпилировался бы и Objective C компилятором, может и этот хаб/тег добавите?

[zelserg]

«То, что код для С компилируется плюсовым компилятором почти всегда и так всем известно». В том то и дело, что раз «почти», то есть исключения. Вот я и написал в заголовке С++ в скобках, показывающее, что примеры в статье к исключениям не относятся.
В заголовке указано: на Си (С++), Objective C не подходит ни под одну из категорий?
Если вы трактуете простой и понятный текст своебразным образом, то причем тут я, заголовок, С и С++?

[stack_trace]

Послушайте, я никак не трактую ваш текст, я просто сделал замечание, с которым, я думаю, согласится любой, для кого C++ — основной язык разработки. Код на С != код на С++, равно как С != С++. Заходя в хаб С++ я ожидаю увидеть код на С++, а заходя в хаб С — код на С. При этом я прекрасно отдаю себе отчёт в том, что почти любой код из хаба С скомпилируется как С++.

Обычно, указывают что код на C++ когда он явно на C++ и ничем больше не скомпилируется. Так можно ещё поставить теги C++11, C++14, C++17. А что? Компиляторами, поддерживающими эти стандарты он тоже скомпилируется.

Вообще, не стоит воспринимать моё замечание так в штыки. Вы просили пояснений — я их дал, я не собирался с вами спорить.

В заголовке указано: на Си (С++), Objective C не подходит ни под одну из категорий?

Конечно же, нет.

[zelserg]

Насчет заблуждений. Я занимаюсь программированием не один десяток лет и давно избавился от заблуждения, которое в грубой формулировке звучит примерно так: «если при написании кода можно выпендриться, то нужно выпендриться».

[stack_trace]

Да причём тут «выпендриться»? Вы о чём вообще?

[zelserg]

Похоже, у нас разный менталитет и мы не понимаем друг друга. Вы в тегах С++ ищете то, что будет работать только на С++, а я в тегах С++ ищу то, что будет работать на С++.

[stack_trace]

Ну, хабр это же не pastbin, тут люди обычно не код ищут а статьи на конкретную тему, узнать что-то новое про конкретную область. Так вот, к теме С++ данная статья имеет крайне слабое отношение, и ничего нового про С++ из этой статьи узнать не получиться, равно как и про Objective C. Про С — возможно, про алгоритмы — наверняка. ИМХО, правильными хабами для этой статьи были бы C и Алгоритмы, но никак не C++.

Только не обижайтесь, ради бога. Никак не хотел вас задеть своим замечанием.

[zelserg]

А разве не вы мне тут минусов наставили? Или какой-то идиот c чувством неполноценности влез в чужое обсуждение?

[zelserg]

Мне даже жалко этого идиота. То, что ему вернется, явно будет несопоставимо с какими-то «минусами». Хотя он, скорее всего, не свяжет свои мелкие гадости с этим возвратом. Что тут скажешь — сам делал, никто не принуждал…

[stack_trace]

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

[semenyakinVS]

Вот-вот. Вот если бы "<название статьи> на этапе компиляции" — тогда другое дело

[zelserg]

Т.е «Расчет биномиальных коэффициентов на Си (С++) на этапе компиляции»? Мне это кажется странным, мягко говоря.

[Eivind]

Нет ничего проще:

#include <type_traits>
#include <limits>

template <class K, class N, bool, bool>
struct bci_helper {
    static constexpr int value = 1;
};

template <class T, T K, T N>
using bci = bci_helper<std::integral_constant<T,K>, std::integral_constant<T,N>, K!=N, K!=0>;

template <class T, T K, T N>
struct bci_helper<std::integral_constant<T,K>, std::integral_constant<T,N>, true, true> {
    static_assert(bci<T, K, N-1>::value <= std::numeric_limits<T>::max() - bci<T, K-1, N-1>::value,"Integer overflow");
    static constexpr T value = bci<T, K, N-1>::value + bci<T,K-1,N-1>::value;
};

[semenyakinVS]

… и вот как-то так были устранены претензии к названию статьи

[Eivind]

Немного хардкора с длинной арифметикой

#include <iostream>
#include <type_traits>

template <char ... Chars>
struct number {};

std::ostream& operator<< (std::ostream& os, number<>){
    return os;
}

template <char C, char ... Chars>
std::ostream& operator<< (std::ostream& os, number<C, Chars...>){
    os << C << number<Chars...>{};
    return os;
}

template <class A, class B>
struct append;

template <char ... As, char ... Bs>
struct append<number<As...>, number<Bs...>> {
    using type = number<As..., Bs...>;
};

template <class Source, class Result = number<>>
struct reverse;

template <char A, char ... As, char ... Bs>
struct reverse<number<A,As...>, number<Bs...>> {
    using type = typename reverse<number<As...>,number<A,Bs...>>::type;
};

template <char ... Bs>
struct reverse<number<>,number<Bs...>> {
    using type = number<Bs...>;
};

template <class T>
struct remove_zeroes {
    using type = T;
};

template <>
struct remove_zeroes<number<'0'>> {
    using type = number<'0'>;
};

template <char ... As>
struct remove_zeroes<number<'0',As...>> {
    using type = typename remove_zeroes<number<As...>>::type;
};

template <char A, char B, char C>
struct basic_sum {
    static constexpr char H = '0' + ((A-'0') + (B-'0') + (C -'0'))/10;
    static constexpr char L = '0' + ((A-'0') + (B-'0') + (C -'0'))%10;
};

template <class A, class B, char C>
struct partial_sum;

template <char A, char ... As, char B, char ... Bs, char C>
struct partial_sum<number<A,As...>, number<B, Bs...>, C> {
    using L = number<basic_sum<A,B,C>::L>;
    using H = typename partial_sum<number<As...>,number<Bs...>,basic_sum<A,B,C>::H>::type;
    using type = typename append<L,H>::type;
};
template <char A, char ... As, char C>
struct partial_sum<number<A,As...>,number<>,C> {
    using L = number<basic_sum<A,'0',C>::L>;
    using H = typename partial_sum<number<As...>,number<>,basic_sum<A,'0',C>::H>::type;
    using type = typename append<L,H>::type;
};

template <char B, char ... Bs, char C>
struct partial_sum<number<>,number<B,Bs...>,C> {
    using L = number<basic_sum<'0',B,C>::L>;
    using H = typename partial_sum<number<>,number<Bs...>,basic_sum<B,'0',C>::H>::type;
    using type = typename append<L,H>::type;
};

template <char C>
struct partial_sum<number<>,number<>,C> {
    using L = number<C>;
    using H = number<>;
    using type = typename append<L,H>::type;
};

template <class A, class B>
using sum = typename remove_zeroes<typename reverse<typename partial_sum<typename reverse<A>::type,typename reverse<B>::type,'0'>::type>::type>::type;

template <class>
struct decrement_helper;

template <char A, char... As>
struct decrement_helper<number<A,As...>> {
    using L = number<(A=='0')?'9':(A-1)>;
    using H = typename std::conditional<A=='0',typename decrement_helper<number<As...>>::type,number<As...>>::type;
    using type = typename append<L,H>::type;
};

template <>
struct decrement_helper<number<>> {
    using type = number<>;
};

template <class T>
using decrement = typename remove_zeroes<typename reverse<typename decrement_helper<typename reverse<T>::type>::type>::type>::type;

template <class K, class N>
struct bcl_helper {
    using type = sum<typename bcl_helper<K,decrement<N>>::type,typename bcl_helper<decrement<K>,decrement<N>>::type>;
};

template <class N>
struct bcl_helper<N,N> {
    using type = number<'1'>;
};

template <class N>
struct bcl_helper<number<'0'>,N> {
    using type = number<'1'>;
};

template <typename CharT, CharT ...String> auto operator "" _n() {
    return number<String...>();
}

template <class K, class N>
auto bcl(K,N) {
    return typename bcl_helper<K,N>::type();
}

int main(int /*argc*/, char** /*argv*/)
{
    std::cout << bcl("50"_n,"100"_n) << std::endl;
}

[encyclopedist]

Вот такая функция считает в 32 битных числах правильно до 33 включительно (проверяем, делится ли, и если да то сначала делим а потом умножаем):

unsigned bci2(int n, int k) 
{
    if (k>n/2) k=n-k; // возьмем минимальное из k, n-k.. В силу симметричность C(n,k)=C(n,n-k)
    if (k==1)  return n;
    if (k==0)  return 1;
    unsigned r=1;
    for (int i=1; i<=k;++i) {
        if(r % i == 0)
        {
    	    r/=i;
    	    r*=n-i+1;
        }
        else
        {
            r*=n-i+1;
            r/=i;
        }
    }
    return r;
}

[zelserg]

Интересное замечание, но хочу заметить:
1. Возможно, при каких-то k и n<=33 все же возникнет переполнение
2. Граница отодвинулась, но для n=34 все же осталось переполнение.
3. Переход на 64 бита кардинально отодвигает эту границу
Можно, конечно, в числителе и знаменмтеле поискать общие множители и посокращать их, но я не стал двигаиться в этом напралении.

[encyclopedist]

Можно, конечно, в числителе и знаменмтеле поискать общие множители и посокращать их, но я не стал двигаиться в этом напралении.

— См. обновление habrahabr.ru/post/274689/#comment_8729927

[encyclopedist]

И ещё лучший выриант (Live):

template<typename T>
T gcd(T a, T b)
{
    while (b > 0)
    {
        a = a % b;
        std::swap(a, b);
    }
    return a;
}

template <typename T>
T bci2(T n, T k) 
{
    if (k > n/2)
        k = n - k;
    if (k == 1)
        return n;
    if (k == 0)
        return 1;
    T r = 1;
    for (T i = 1; i <= k; ++i) {
        T f = gcd(r, i);
        r /= f;
        r *= n - i + 1;
        r /= (i/f);
    }
    return r;
}

С unsigned int работает до 33, с unsigned long работает до 66.

[encyclopedist]

Обновление. Теперь точно все множители сокращает:

template <typename T>
T bci2(T n, T k) 
{
    if (k > n/2)
        k = n - k;
    if (k == 1)
        return n;
    if (k == 0)
        return 1;
    T r = 1;
    for (T i = 1; i <= k; ++i) {
        T a = n - i + 1;
        T b = i;
        // compute r = r * a / b;
        T f = gcd(r, i);
        r /= f;
        b /= f;
        T f1 = gcd(a, b);
        r *= a/f1;
        r /= b/f1;
    }
    return r;
}

Работает до 34 и 67 соответственно.

[zelserg]

«Простая задача часто имеет простое и понятное решение. Увы, неправильное...» Как увидеть эту неправильность и довести до ума — вот в чем вопрос. Способы могут быть разными. У меня работает и при n=67 и не тратит время на вычисление НОДов.

[encyclopedist]

Все зависит от требований. У меня без доп памяти. С доп памятью надо ещё посмотреть что лучше (на вскидку и моё решение и динамическое программирование работают за O(n^2)).

[encyclopedist]

Моё даже лучше: gcd считается за логарифмическое время, поэтому O(n log n).

[encyclopedist]

Вот без НОДов (трюк r*a/b = (r/b)*a + (r%b)*a/b из этого комментария Mrrl ). Считает без переполнения пока это теоретически возможно. (на Coliru)

Заголовок спойлера

template <typename T>
T bci2(T n, T k) 
{
    if (k > n/2)
        k = n - k;
    T r = 1;
    for (T i = 1; i <= k; ++i)
    {
        // r * a / b = (r/b)*a + (r%b)*a/b
        T a = n - i + 1;
        r = (r/i)*a + ((r%i)*a)/i;
    }
    return r;
}

[zelserg]

Спасибо за трюк и за ссылку. Надо же — перед написанием статьи искал на хабре «биномиальные», а надо было «мультиномиальные»

[a_batyr]

Естественно, пришлось использовать «длинную арифметику» (свою). Для этой задачи я написал «рекурсию с памятью»: вычисленные коэффициенты запоминались в массиве и экспоненциального роста времени расчета не было.

Во-первых стоило сразу привести в статье длинную арифметику. Это грустно, когда оперируют с double при расчёте целочисленных значений

return ceil(r-0.2); // округлим до ближайшего целого, отбросив дробную часть

А во-вторых, можно привести вычисление коэффициентов без факториалов, а динамически оперируя только суммой.

[zelserg]

Напишите, как. Кроме приведенной в статье рекурсии, естественно.

[a_batyr]

Расчёт БК без умножения, динамически.
Храним массив БК для элементов C(i,0), C(i,1), ..., C(i,i). Инициализируем для i=0 как {1}. Каждую следующую итерацию вычисляем через предыдущие значения C(i,k) = C(i-1,k)+C(i-1,k-1) вплоть до заданного C(n,k).
Вы написали что-то похожее, только рекурсивно, что слишком прожорливо. Вам не нужно повторно вычислять промежуточные значения, а в случае рекурсии вы это делаете. Храня уже известные значения вы на порядок улучшите производительность.

[zelserg]

Я сделал добавление в статье — рекурсивная функция bcr с запоминанием.
Она, конечно, самая быстрая, но значимого выигрыша в скорости для int и long long не дает. Значимого — т.е. чтобы bci/bcl, были узким местом в какой-то задаче, а новая bcr расшивала это узкое место.

[chersanya]

Странное у вас понимание «значимости» ускорения :) Если функция работает например в несколько раз быстрее — это уже значимое ускорение, хотя при правильном решении на олимпиаде это обычно не повлияет, т.к. таймлимит обычно ставится значительно больше необходимого. Измерили бы скорости всех вариантов да указали в статье, было бы интересно сравнить.

[zelserg]

«Во-первых стоило сразу привести в статье длинную арифметику. Это грустно, когда оперируют с double при расчёте целочисленных значений»
Сначала о грусти. Для практических расчетов биномиальных коэффициентов (БК) double — самое простое и эффективное решение. Если вы моделируете случайные процессы, то точности в 16 знаков «за глаза» хватает и выборки из 1000 элементов, скорее всего, тоже хватит. Кроме того, БК в реальных расчетах живут не сами по себе, а умножаются на вероятности, которые сами по себе double. Грустно городить длинную арифметику для расчета БК, а потом умножать их на double.
Длинная арифметика нужна при расчете БК только в одном случае — при решении олимпиадных задач и это тема для отдельной статьи (не моей). Я считаю, что каждый потенциальный олимпиец должен написать свою длинную арифметику, не такая уж она и сложная. А коль не может — значит олимпиады не для него.

[a_batyr]

Дальнейшее повышение точности и расчет при n>67

… используем double ...

Расхождение расчета bcd с точным значением начинается с n=57.

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

Для экстремалов и «олимпийцев»

и совершенно не раскрываете тему статьи для олимпийцев. И уберите кавычки с этого слова ради бога, олимпийцы это живые люди!
Что ж, некоторые вообще считают нормальным снимать вертикальное видео и я на это никак не могу повлиять :-(

[zelserg]

«И уберите кавычки с этого слова ради бога». Облимпийцы без кавычек — это Олимпийские боги (уж они то не одобрят убирание кавычек) в греческой мифологии и спортсмены-участники Олимпийских игр. К участникам прочих олимпиад слово олимпийцы не относится.

В моей статье я сделал 2 небольших замечания для желающих заниматься олимпиадным программированием.
1. Привел улучшенный вариант функции bcl, который правильно считает все БК, попадающие в unsigned long long
2. Для тех, кому потребуется точное вычисления очень больших биномиальных коэффициентов («длинная арифметика»), посоветовал использовать рекурсивный расчет с запоминанием. Это все, я хотел сказать и сказал. Кто не понял — я не виноват.

[zelserg]

12 человек добавили в избранное, но только 2 плюсанули. Итог — минус 2. Пожалуй, не стоит писать подобные статьи (это был пробный камень. так сказать).

[samo-delkin]

Ну, вот ты написал

		r=bcl(n-1,k);
		r+=+bcl(n-1,k-1);

это же можно в одну строку записать.
То есть не знаешь ни C, ни C++. Алгоритм тоже неоптимальный.
На олимпиады ссылаешься, а зачем они нужны? Пользы от них никакой нет.

[a_batyr]

Это тот случай когда в комментариях вопрос раскрыт намного лучше, чем в самой статье.

[zelserg]

Каждый видит и понимает «по способностям».

[zelserg]

В каждой статье есть шероховатости. Раздельный расчет нужен был для анализа, где переходить к рекурсивному расчету.

Критиковать гораздо легче, чем написать свою статью. На любую, самую замечательныую статью найдутся критики и недоброжелатели. Не навится — не читайте и не комментируйте. Значит она не для вас. Кто-то комментирует, чтобы статья стала лучше (и я вношу исправления по ходу обсуждения), а кто-то с целью нагадить и выпятить свое Я.

В данному моменту 26 человек поставили ее в избранное, значит посчитали ее интересной и/или полезной. Вот для таких она и написана. А любителям выискивать соринки в чужих глазах посоветую заняться своими.

«На олимпиады ссылаешься, а зачем они нужны? Пользы от них никакой нет.» Вопрос спорный. Не смог решить задачу, тебя не взяли в команду и обида гложет?

[samo-delkin]

«На олимпиады ссылаешься, а зачем они нужны? Пользы от них никакой нет.» Вопрос спорный. Не смог решить задачу, тебя не взяли в команду и обида гложет?

Не, я просто решал такие задачи и не получил от этого ничего. Зря время только потратил.
Можно, конечно, сидеть и ждать, когда тебя в команду возьмут, а можно просто сидеть и делать программы, которые потом работают вместо тебя месяцами. (Хотя я пользуюсь уже несколько лет.)

[KvanTTT]

А еще есть: Методы вычисления мультиномиальных коэффициентов.

[zelserg]

Хорошая статья, добавил ее в избранное. Метод рекурсивного расчета с запоминанием я когда-то делал и для int и для long long. Он быстрый, ест не так уж много памяти. Но в результате остановился на приведенной к конце статьи bcl. Запоминание осталось только для длинной арифметики, когда надо было посчитать очень много очень больших коэффицентов (да и операция деления в длинной арифмиетике — то еще «удовольствие»).

Треугольник Паскаля интересная штука — очень быстрая, если запоминать промежуточные результаты и очень медленная, если не запоминать.

[middle]

А почему вы не считаете треугольник Паскаля? Простое быстрое итеративное решение.

[zelserg]

Что вы имеете в виду под «почему вы не считаете треугольник Паскаля»?

Я написал в статье функцию bcr, поторая считает по треугольнику паскаля. И написал, что элемент из 34 строки считался около минуты. Также в конце статьи написал, что при расчете коэффициентов с длинной арифметикой использовал рекурсию (т.е. треугольник паскаля) с запоминанием посчитанных коэффициентов.

[middle]

Я имел в виду, почему вы не используете итеративный алгоритм построения треугольника Паскаля, храня строку в массиве? 70 строк у меня на ноуте считаются за долю секунды.

Upd: и это при том, что я считаю каждый элемент отдельно — чтобы жизнь мёдом не казалась. А если печатать всю посчитанную строку целиком — будет вообще миллисекунды.

[zelserg]

Я писал процедуры расчета по треугольнику с запоминанием, но не стал приводить их в статье.
1. Для int и long long заметной глазу разницы в скорости нет (ну не замечает мой глаз миллисекунды).
2. Для олимпиадных задач, в которых надо было считать коэффициенты разницы тоже не было.
3. Код с запоминанием не такой простой и чуть более объемный (зачем искать/приводить путь в два шага, если есть путь в один шаг).

[chersanya]

1. Для int и long long заметной глазу разницы в скорости нет (ну не замечает мой глаз миллисекунды).

2. Для олимпиадных задач, в которых надо было считать коэффициенты разницы тоже не было.

Странно видеть рядом с упоминанием олимпиадных задач оценку скорости на глаз :)

3. Код с запоминанием не такой простой и чуть более объемный (зачем искать/приводить путь в два шага, если есть путь в один шаг)

А по-моему, он почти не отличается от кода без запоминания, если будем хранить весь треугольник. Заодно автоматически получается кеширование — если несколько раз попросим одно и то же значение, оно не будет пересчитываться.

Для сравнения — ваш код из статьи:

unsigned bcr(int n,int k) 
{
	if (k>n/2) k=n-k;
	if (k==1)  return n;
	if (k==0)  return 1;
	return bcr(n-1,k)+bcr(n-1,k-1);
}

И простейший вариант с запоминанием (не компилировал, но должно работать):

unsigned bcr_cache[N_MAX][K_MAX] = {0};

unsigned bcr(int n,int k) 
{
	if (k>n/2) k=n-k;
	if (k==1)  return n;
	if (k==0)  return 1;
        if (bcr_cache[n][k] == 0)
            bcr_cache[n][k] = bcr(n-1,k)+bcr(n-1,k-1);
        return bcr_cache[n][k];
}

Кстати, в вашей строчке

return ceil(r-0.2); // округлим до ближайшего целого, отбросив дробную часть

неправильный либо код, либо комментарий.

[zelserg]

1. Я не сравниваю «на глаз» с олимпиадным программированианием. «На глаз» относится к замечанию «на моем ноутбуке», про миллисекунды — типа юмор. Про олимпиадное — то что тесты проходят и оптимизации скорости расчета БК не требуется.
2. При запоминании фиксированные размеры мне как-то не очень… Но вариант вполне простой и рабочий.
return ceil(r-0.2); // округлим до ближайшего целого, отбросив дробную часть
неправильный либо код, либо комментарий.

Комментарий мог быть таким: я сравнивал double с точными расчетами и «экспериментально» нашел, что это значение подходит (чтобы расхождение наступило как можно позже). Наверно, можно заменить на 0.4 или другую константу, но лучше не будет.

Кстати, в статье есть еще одно «экспериментально найденная» константа: if(n+k)>90 для перехода к рекурсивной ветке.

[chersanya]

Ну фиксированные размеры я поставил только для простоты, и чтобы код оставался кодом на Си (у вас как я понял старательно не используется ничего из С++). Можно сделать вектором и добавлять туда по надобности, хотя раз N всё равно < 70, то можно независимо от задачи делать размер 70х70.

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

Комментарий в таком случае вводит в заблуждение — видя округление к ближайшему целому читатель ожидает ceil(r — 0.5) или floor(r + 0.5), или лучше всего round. И кстати не понятно, зачем это вообще — если N такие, что переполнения нет, то можно считать с int64. Если же N большие и мы используем double, то значит точность до единицы не нужна — она всё равно не будет получена при больших N.

Кстати, если так убрать ваш ceil и возвращать просто r, то функцию можно сразу сделать шаблонной — чтобы не было такого, что «заменим unsigned на double», а одна функция работала и для int32, и для int64, и для double и для длинной арифметики.

[zelserg]

Если просто возращать r, то у биномиальных коэффициентов появляются дробные хвосты и мне это не понравилось (все же БК — по определению целые).

[zelserg]

Изменил комментарий в статье.

[chersanya]

И ещё, для как вы называете «олимпийцев» вычисление биномиальных коэффициентов неплохо написано на e-maxx.ru/algo/binomial_coeff, причём в разных вариантах.

[zelserg]

Интересная стстья. При «длинных» расчетах запоминать последние две строки приходилсь, а вот запоминать факториалы — никогда. Как-то выкручивался или мне просто не попалась задача, где это необходимо

[zelserg]

Интересно, кто и за что заминусовал этот комментарий. Похоже, кто-то просто наставил минусов из желания нагадить. Ладно, они ему так или иначе вернутся.

[middle]

Ну хоть упомяньте их в статье, а то складывается впечатление, что либо быстро, но с переполнением, либо медленно, но без него. Разбирать задачу и даже не упомянуть характеристики столь важного варианта — как-то удивительно.

Про int и long long я ничего не говорил :)

[middle]

Собственно, о сложности столь простого куска кода и рассуждать-то странно:

  // typedef uint32_t bino_int;
  // typedef uint64_t bino_int;

  buffer[0] = 1;

  for (size = 1; size <= n; ++size) {
    bino_int prev = 1; // buffer[0]
    for (size_t off = 1; off < size; ++off) {
      bino_int tmp = buffer[off];
      buffer[off] += prev;
      prev = tmp;
    }
    buffer[size] = 1;
  }

[zelserg]

Написал добавление в статье

[cybrid]

Достаточно странная статья. Рассмотрены только «кривые» решения, а рекурсивное решение с запоминанием и прямое решение с подсчетом простых множителей не приведены.

[zelserg]

Привел рекурсивное с запоминанием, а подсчет множителей я не считаю «прямым».

[dprotopopov]

Мой ответ на Вашу публикацию habrahabr.ru/post/274729 — Расчет биномиальных коэффициентов с использованием Фурье-преобразований

[zelserg]

Давайте сравним скорость расчета C(67,33). Можете кинуть проект, который считает этот коэффицент и который можно открыть в VS2010?

[zelserg]

Или сами сравните. Полагаю, перевести bcr с запоминанием на c# несложно.

[dprotopopov]

взятое возвращаю кармой.
Можете сравнить и описать на Хабре эксперимент для разных типов данных, для C++ и C# — многих заинтересует.

[Ryder95]

мне кажется, первый код работать не будет, так как переменная r в цикле инициализируется, то после цикла её нет, поэтому вернуть r невозможно

[zelserg]

Спасибо, поправил. Намудрил, когда переводил переделал факториал из рекурсии в цикл.

[ShashkovS]

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

Если известны ограничения на n, то простые числа можно захардкодить, иначе добыварь решетом Эратосфена.
После этого останется только вычесть показатели степеней и перемножить. При этом гарантируется, что самое большое число, которое встретится в вычислениях — это ответ.

Заодно такая техника позволит вычислять очень большие C_n^k

Вот код на Python:

from timeit import timeit

def pow_binomial(n, k):
    def eratosthenes(N):  # Возвращает список простых чисел от 2 до N
        simp = [2]  # Начинаем с 2
        nonsimp = set()  # Составные числа. Пока пусто
        for i in range(3, N + 1, 2):  # Проходим все нечётные числа
            if i not in nonsimp:  # Если оно не составное, значит, новое простое
                nonsimp |= {j for j in range(i * i, N + 1, 2 * i)}  # Добавляем все кратные
                simp.append(i)  # А само число в список простых
        return simp
    def calc_pow_in_factorial(a, p):  # Вычисляет кратность вхождения простого числа p в число a!
        res = 0
        while a:
            a //= p
            res += a
        return res
    # Не паримся и используем формулу n!/(k!(n-k)!)
    simple = eratosthenes(n+1)  # Берём список простых
    n_fact_pows = {p: calc_pow_in_factorial(n, p) for p in simple}
    k_fact_pows = {p: calc_pow_in_factorial(k, p) for p in simple}
    nmk_fact_pows = {p: calc_pow_in_factorial(n - k, p) for p in simple}
    ans = 1
    for p in simple:
        ans *= p ** (n_fact_pows[p] - k_fact_pows[p] - nmk_fact_pows[p])
    return ans


def naive_factorial(n):
    res = 1
    for i in range(2, n + 1):
        res *= i
    return res


def naive_binomial(n, k):
    return naive_factorial(n) // naive_factorial(k) // naive_factorial(n - k)


REPEATS = 10
N = 50000
K = 10000
setup = 'from __main__ import N, K, naive_binomial, pow_binomial'
print(timeit(stmt='pow_binomial(N, K)', setup=setup, number=REPEATS)/REPEATS)
print(timeit(stmt='naive_binomial(N, K)', setup=setup, number=REPEATS)/REPEATS)

>>> 
0.02103983737743497
1.1220947057368764

То есть C₅₀₀₀₀¹⁰⁰⁰⁰ вычисляется за 0.02 секунды, а там ещё есть, что ускорять.

При n <= 71 если убрать комментарии и сократить сопли, получается 11 строк:

def pow_binomial(n, k):
    def calc_pow(a, p):
        res = 0
        while a:
            a //= p
            res += a
        return res
    ans = 1
    for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71]:
        ans *= p ** (calc_pow(n, p) - calc_pow(k, p) - calc_pow(n - k, p))
    return ans

[sci_nov]

Дальнейшее повышение точности и расчет при n>67

Для экстремалов и «олимпийцев»

www.artlebedev.ru/kovodstvo/sections/136
Абзацы 07 и 08.

[zelserg]

Добавил в статью свой расчет на Python, чтобы сравнить с Си и расчетом в: комментарии