Accurate Binomial Coefficients

Binomial coefficients / Pascal Triangle elements for / N choose K look-up table (LUT) with no numerical errors for 0 <= N <=1029.

Syntax

b = binomial(n)

b = binomial(n, k)

b = binomial(n, [])

Description

b = binomial(n) generates a table with the binomial coefficients up to order 'n', with 0 <= n <= 1029.

b = binomial(n, k) returns the 'N' over 'K' coefficients in n(:) and k(:) or in n(:, 1), n(:, 2).

b = binomial(n, []) returns the 'N' + 1 -th row of the Pascal Triangle, or the rows in n(:) + 1.

The values are pre-computed using arbitrary numerical precision and therefore are exempt of overflows. Due to the numerical precision of the 'double' class type, the values have been saved as a cell array with strings. Conversion to double can be easily achieved through the str2double command, at the cost of loss of numerical precision.

Examples

% Pascal Triangle for n <= 4
b = binomial(4);

b =

5×5 cell array

{'1'} {'0'} {'0'} {'0'} {'0'}

{'1'} {'1'} {'0'} {'0'} {'0'}

{'1'} {'2'} {'1'} {'0'} {'0'}

{'1'} {'3'} {'3'} {'1'} {'0'}

{'1'} {'4'} {'6'} {'4'} {'1'}

% Binomial expansion terms for n = 4
b = binomial(4, []);

b =

1×5 cell array

{'1'} {'4'} {'6'} {'4'} {'1'}

% Accurate N choose K (1000, 500)
b = round(vpa(binomial(1000, 500)));

b =

270288240945436569515614693625980000000000000092092817939260465262396897314854756510301009414481288655741593673629116029259908377031297694048974886875144010566190142039764976194007777780262204642682840529729586602474213985385651649048525534270439972471247027710107116837053571858737469793700061118464

% Table of binomial coefficients
x = [3:3:12; 0,3,6,8];
b = binomial(x);     % Equal to: binomial(x'), or binomial(x(:, 1), x(:, 2))

b =

4×1 cell array

{'1' }

{'20' }

{'84' }

{'495'}

% Some binomial elements for a given 'n'
b = binomial(10, 0 : 2 : 10);

b =

1×6 cell array

{'1'} {'45'} {'210'} {'210'} {'45'} {'1'}

% Overflow and loss of numerical precision in 'double' types
b = binomial(60, 30);
err = str2double(b) - nchoosek(60, 30);
b = binomial(65, 35);
res = sym(b) - str2double(b);

err =

16

res =

-32