О вычислении биномиальных коэффициентов в целочисленной арифметике

Вычислить их сколько-то заранее, и хранить в виде таблицы. Самый быстрый вариант :-)

Стоит ещё сделать проверку на корректность аргументов, иначе при к>н оно зависает.

Ну, достаточно заменить (= k n) на (>= k n). По крайней мере, уже на зависнет.

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

#include <stdio.h>

typedef unsigned long long int  bc_t;
typedef unsigned char           index_t;

static bc_t
binomial_iter(bc_t n, bc_t k, bc_t i, bc_t prev) {
  return i == k ? prev : binomial_iter(n, k, i + 1, (n - i) * prev / (i + 1));
}

bc_t
binomial(bc_t n, bc_t k) {
  if ((k == 0) || (n == k)) return 1;
  return k < n-k ? binomial_iter(n, k, 0, 1) : binomial_iter(n, n - k, 0, 1);
}

int main() {
  index_t n, k;
  unsigned int i = 1;
  for (n = 1; n < 100; n++)
    for (k = 1; k < 100; k++)
      if (n >= k) {
        printf("%d: %d %d %llu\n", i, n, k, binomial(n, k));
        i++;
      }
  return 0;
}

Даст меньше 2000 ответов, а мемоизированный вот:

#include <stdio.h>
#include <assert.h>

typedef unsigned long long int  bc_t;
typedef unsigned char           index_t;

#define	BC_TABLE_SIZE  68

/* ~35 KB of memory */
bc_t bc_table[BC_TABLE_SIZE][BC_TABLE_SIZE];

void
fill_bc_table(void) {
  index_t i, j;
  for (i = 0; i < BC_TABLE_SIZE; i++) {
    bc_table[i][0] = 1;
    bc_table[i][i] = 1;
  }
  for (i = 1; i < BC_TABLE_SIZE; i++)
    for (j = 1; j < i; j++)
      bc_table[i][j] = bc_table[i-1][j-1] + bc_table[i-1][j];
}

bc_t
binomial(bc_t n, bc_t k) {
  assert((n >= k) && (n < BC_TABLE_SIZE));
  return bc_table[n][k];
}

int main() {

  fill_bc_table();

  index_t n, k;
  unsigned int i = 1;
  for (n = 1; n < 68; n++)
    for (k = 1; k < 68; k++)
      if (n >= k) {
        printf("%d: %d %d %llu\n", i, n, k, binomial(n, k));
        i++;
      }
  return 0;
}

даёт больше, причём на заполнение таблицы тратится время порядка миллисекунд, а все дальнейшие операции binomial - O(1).

$ gcc -O2 binomial_1.c -o binomial_1
$ ./binomial_1 > binomial_1.log

$ gcc -O2 binomial_2.c -o binomial_2
$ ./binomial_2 > binomial_2.log

В принципе, 35 KB на статический массив это нормально, тем более что больше тут никак не может быть, так как unsigned long long int начинает переполняется. Если бы имело смысл делать больше и экономить - можно было бы сделать linked-list содержащие realloc-двумерные-массивы, и время от времени туда добавлять элементов (в зависимости от гранулировки шага мемоизации).

А вот если брать bignum алгоритмы, то (опять же - насколько я знаю) лучше и эффективнее реализуется:

     N!            1 * ... * N          (K+1) * ... * N
  --------  =  --------------------  =  ---------------,  N>1
  K!(N-K)!     1 * ... * K * (N-K)!         (N-K)!

Т.е. произведение ряда + целочисленное деление + факториал. К сожалению при фиксированных числах тут очень быстро получится переполнение разрядной сетки.

use ММИ, Luke!

>>use ММИ, Luke!

Математик, наверное? Потому что как из анекдота: твой совет абсолютно точен и абсолютно бесполезен (c).

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

Например, в случае мультипликативной реализации, представленной ниже, при n от 1 до 200, только в 4 случаях возникли проблемы: при (n,k) равных (150,136) (150,137) (151,138) (181,176).

Здесь, функция binomial-overflows? проверяет, произошло ли переполнение при вычислении C(n,k). Переполнением считается тот случай, когда конечный результат ещё укладывается в интервал [0,2^64), но один из промежуточных результатов лежит вне его.

Ключевой момент этой функции в том, что при вычислении C(n,k) мы пользуемся той из двух формул:

C(n,k)=n/(n-k)*C(n-1,k)

C(n,k)=(n-k+1)/k*C(n,k-1)

которая дает промежуточный результат (рациональное число) с меньшим «весом». Эта стратегия достаточно прямолинейна, и может быть улучшена.

Тестовая функция enumerate-overflows в качестве ответа возвращает набор пар, вида (n (k[1] k[2] ... k[s])) таких, что при вычислении C(n,k[i]) произошло переполнение в указанном выше смысле.

Соответствующая функция binomial имеет вид:

(define (binomial n k)
  (define (iter result n k)
    (cond ((or (= k 0) (= k n)) result)
          (else
           (let ((q1 (make-rat n (- n k)))
                 (q2 (make-rat (+ 1 (- n k)) k)))
             (let ((r1 (mul result q1))
                   (r2 (mul result q2)))
               (if (< (weight r1) (weight r2))
                   (iter r1 (- n 1) k)
                   (iter r2 n (- k 1))))))))

  (numer (iter (make-rat 1 1) n k)))

И как мы видим, соответствующая часть, реализующая упрощённый интерфейс для рациональных чисел, занимает всего 12 строк.

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

Будем рассматривать упорядоченные наборы натуральных чисел, длины K.

Увеличение каждого элемента набора X на s обозначим так: X+s Циклический сдвиг набора X влево на s обозначим так: X<=s

Будем говорить, что набор X соответствует набору Y, если Y[i] делится на X[i] для всех i.

Рассмотрим наборы натуральных чисел Q=[1 2 ... K] и R=[1 2 ... K].

Теорема. Q<=N соответствует набору R+N для любого N.

Доказательство:

Q<=0 соответствует R+0, по условию.

Предположим Q<=(N-1) соответствует R+(N-1). Обозначим Q<=(N-1) как Y, R+(N-1) как X.

Покажем, что Q<=N соответствует набору R+N. Обозначим Q<=N как S, R+N как T.

По предположению индукции для 1<=i<=K выполнено X[i] делится на Y[i]. По определению S совпадает с Y<=1 и T совпадает с X+1. Значит для 1<=i<=(K-1) выполнено T[i] делится на S[i]. Осталось показать, что T[K] делится на S[K].

T[K] = X[K]+1 = (X[1]+K-1)+1 = X[1]+K, S[K] = Y[1] по определению.

Но X[1] по предположению индукции делится на Y[1]. Теорема доказана.

Поскольку произвольный набор K подряд идущих чисел [N+1, ... N+K] равен Q+N, то набор R<=N соответствует ему. Но из этого следует, что произведение (N+1)*...(N+K) делится на произведение R[1]*...*R[K]=K!.

Бред! Доказательство неправильное!

Почему бы не использовать библиотеки? Неужели для scheme нет рациональных чисел, filter и map из какой-нибудь коробки?

Ещё в scheme есть большие числа, поэтому эмулировать на ней проверку переполнения не совсем правильно - лучше сразу использовать заведомо лучшие алгоритмы для больших чисел.

В мультипликативных алгоритмах написанных на си с использованием целочисленной арифметики (сабж) переполнение возникает раньше чем при использовании той аддитивной рекурсивной формулы и ручной мемоизации.

Хотя насчёт решениях в рациональных дробях, представимых в виде, например, массива uint64 rat[2] - возможно действительно можно получить ещё больше ответов. Есть смысл попробовать (ну и Кнута посмотреть).

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

Вы понимаете, почему именно на Scheme написан этот пример? Не для того, чтоб расчитывать C(n,k), используя длинные числа. В этом примере существенно не используются возможности Scheme для работы с длинными числами.

Вы можете переписать его на C с незначительными изменениями при проверке переполнения (надеюсь знаете, как по двум числам заранее сказать - будет ли переполнение при перемножении, не перемножая их?)

Внимание!

Этот пример показывает, что при правильно реализованном мультипликативном алгоритме очень редки случаи, когда C(n,k) находится в интервале [0,2^64), но один из промежуточных результатов вне его.

Тот факт, что пример реализован на Scheme, здесь ни при чем.

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

Короче говоря, в этом примере длинные числа Scheme используются только для удобной проверки переполнения результата.

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

Вы видите, что при n от 1 до 100 в моей реализации мультипликативного алгоритма никогда не бывает так, чтоб конечный результат был в [0,2^64) а промежуточные выходили из него?

А в вашей реализации такие случаи есть. Вот и весь сказ.

Вы можете переписать его на C с незначительными изменениями

Фух, тогда ладно, а то я испугася уже.

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

Это не мой алгоритм, это мультипликативный алгоритм в целочисленной арифметике (то что было нужно TC), который был приведён немного ранее. «Мой» (опять же, никакой он не мой, если подумать) алгоритм аддитивный, но он даёт правильные ответы что-то вроде до C(67,_), если в рациональных числах вида (int64, uint64) можно добраться до C(100, _) тогда, по любому, TC нужно брать «ваш» алгоритм.

Вы видите, что при n от 1 до 100 в моей реализации мультипликативного алгоритма никогда не бывает так, чтоб конечный результат был в [0,2^64) а промежуточные выходили из него?

Хотя нет, вы меня запутали. У меня-то алгоритм аддитивный - если что-то начало выходить за границы разрядной сетки, значит сами C(n,k) начали выходить за её границы (а не промежуточные результаты, как при мультипликативном алгоритме).

Например C(68, 34) = 28453041475240576740 что уже не влезает в 64 (unsigned прошу заметить) бита. Поэтому «мой» алгоритм считает до C(67, _). А как «ваш» получает ответы «от 1 до 100», когда уже начиная с 68-ого цикла сами коэффициенты не влезают в 64 бита?

Мой алгоритм отслеживает переполнения только при вычислениях тех C(n,k) которые сами укладываются в машинное слово. Как мы видим, при n>68 такие C(n,k) существуют. Например, C(150,136)=17910906207136060650, что меньше чем 2^64=18446744073709551616.

Так вот, оказалось, что при вычислениях тех C(n,k) которые сами укладываются в машинное слово, хорошо реализованный мультипликативный алгоритм дает переполнения в очень редких случаях.

Я надеюсь ещё немного подумать, и придумать мультипликативный алгоритм, возможно с откатами, ограниченными по длине, который будет так же корректен, как и аддитивный алгоритм, т.е. при вычислении тех C(n,k) которые сами укладываются в машинное слово, он не будет давать переполнений.

Я не отслеживаю возникновение переполнений при вычислении C(n,k) которые сами не укладываются в машинное слово. Для таких C(n,k) не имеет значения, какой алгоритм использовать.

В то же время, существуют C(n,k) которые при n>67 укладываются в машинное слово. Например, C(150,136)=17910906207136060650 что меньше, чем 2^64=18446744073709551616.

Таким образом мой алгоритм «получает» не все ответы при n от 1 до 100, а только те, которые сами укладываются в машинное слово. При этом, переполнений промежуточных результатов в нем не происходит.

Мой алгоритм отслеживает переполнения только при вычислениях тех C(n,k) которые сами укладываются в машинное слово. Как мы видим, при n>68 такие C(n,k) существуют.

Ну да - числа C(n, k) при данном n в интервале [1, 2, ..., k <= n] растут от минимальных значений (от единицы) до некоторых максимальных и уменьшаются обратно к минимальным (и к единице). Т.е. находить их можно при сколь угодно больших n, но при n > 67 «центральные» k начинают обрезаться, а k близкие к нулю и к n - влезают в границы. С аддитивным алгоритмом тоже самое. Вот например:

#include <stdio.h>
#include <stdint.h>
#include <assert.h>

/*
 * Binomial coefficients
 */

/* find from C(1, _) to C(99, _) */
#define BC_TABLE_SIZE  100
#define UNDEFINED      0

typedef uint64_t       bc_t;
typedef uint16_t       index_t;

bc_t bc_table[BC_TABLE_SIZE][BC_TABLE_SIZE];

void
fill_bc_table(void) {
  index_t i, j;
  for (i = 0; i < BC_TABLE_SIZE; i++) {
    bc_table[i][0] = 1;
    bc_table[i][i] = 1;
  }
  for (i = 1; i < BC_TABLE_SIZE; i++)
    for (j = 1; j < i; j++) {
      bc_t sum = bc_table[i-1][j-1] + bc_table[i-1][j];
      /* check overflow */
      if ((sum < bc_table[i-1][j-1])        || (sum < bc_table[i-1][j]) ||
          (bc_table[i-1][j-1] == UNDEFINED) || (bc_table[i-1][j] == UNDEFINED))
        bc_table[i][j] = UNDEFINED;
      else
        bc_table[i][j] = sum;
    }
}

bc_t
binomial(bc_t n, bc_t k) {
  assert((n >= k) && (n < BC_TABLE_SIZE));
  return bc_table[n][k];
}

int main() {

  /* some milli-seconds for table, some 10^1 kb's */
  fill_bc_table();

  index_t n, k, i = 1;
  for (n = 1; n < BC_TABLE_SIZE; n++)
    for (k = 1; k < BC_TABLE_SIZE; k++)
      if (n >= k) {
        printf("%d: %d %d %llu\n", i, n, k, binomial(n, k));
        i++;
      }

  return 0;
}

пишет нули вместо значений больше 2^64, получается такая картинка - http://pastebin.com/q8ue1gGW (~100kb).

Можно в принципе искать дальше C(100, _), но там полезных значений меньше, а памяти аддитивный алгоритм просит - по 80kb на каждые 100 «циклов» коэффициентов.

При этом, переполнений промежуточных результатов в нем не происходит.

Где доказательство отсутствия? :)

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

Но мультипликативный алгоритм при 1<=n<= 500, 0<=k<=n дает около 50 промежуточных переполнений, а при n<=100 не дает промежуточных переполнений даже в простейшем виде, и не требует 80 килобайт на каждые 100 значений параметра n.

Где доказательства отсутствия переполнений?

Т.е. необходимо предоставить формальное доказательство того, что функция (binomial-overflows? n k) возвращает true тогда и только тогда, когда:

1) C(n,k) принадлежит интервалу [0,2^64), и

2) во время её работы происходит переполнение промежуточного результата, т.е. один из параметров вызываемых в теле функций или одна из временных переменных в окружении let лежит вне интервала [0,2^64)

Правильно я понимаю? Подождите, подготовлю - это сделать не трудно.

Функция (enumerate-overflows from to) реализует вложенные циклы:

for (n = from; n <= to; n++) {
    for (k = 0; k <= n; k++)
        // если при вычислении С(n,k) произошло переполнение
        // промежуточного результата, то
        // добавить пару (n,k) в возвращаемый список
        // (на самом деле - результат формируется немного хитрее,
        // но суть именно в этом)
}

Так вот, в после вызова (enumerate-overflows 1 100) результирующий список пуст. Думаю, в корректности двух циклов вы не сомневаетесь, поэтому осталось доказать корректность функции binomial-overflows?

Это я сделаю чуть позже, как написал в предыдущем посте.