Combinations and permutations

This page uses content from Wikipedia. The original article was at Combination. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

This page uses content from Wikipedia. The original article was at Permutation. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

Task

Implement the combination (ⁿC_k) and permutation (ⁿP_k) operators in the target language:

$^{n}\operatorname {C} _{k}={\binom {n}{k}}={\frac {n(n-1)\ldots (n-k+1)}{k(k-1)\dots 1}}$
$^{n}\operatorname {P} _{k}=n\cdot (n-1)\cdot (n-2)\cdots (n-k+1)$

See the wikipedia articles for a more detailed description.

To test, generate and print examples of:

A sample of permutations from 1 to 12 and Combinations from 10 to 60 using exact Integer arithmetic.
A sample of permutations from 5 to 15000 and Combinations from 100 to 1000 using approximate Floating point arithmetic.
This 'floating point' code could be implemented using an approximation, e.g., by calling the Gamma function.

See Also:

Evaluate binomial coefficients

The number of samples of size k from n objects.

With combinations and permutations generation tasks.
	Order Unimportant	Order Important
Without replacement	${\binom {n}{k}}=^{n}\operatorname {C} _{k}={\frac {n(n-1)\ldots (n-k+1)}{k(k-1)\dots 1}}$	$^{n}\operatorname {P} _{k}=n\cdot (n-1)\cdot (n-2)\cdots (n-k+1)$
Without replacement	Task: Combinations	Task: Permutations
With replacement	${\binom {n+k-1}{k}}=^{n+k-1}\operatorname {C} _{k}={(n+k-1)! \over (n-1)!k!}$	$n^{k}$
With replacement	Task: Combinations with repetitions	Task: Permutations with repetitions

ALGOL 68

Works with: ALGOL 68 version Revision 1 - one minor extension to language used - PRAGMA READ, similar to C's #include directive.

Works with: ALGOL 68G version Any - tested with release algol68g-2.6.

File: prelude_combinations_and_permutations.a68<lang algol68># -*- coding: utf-8 -*- #

COMMENT REQUIRED by "prelude_combinations_and_permutations.a68" CO

 MODE CPINT = #LONG# ~;
 MODE CPOUT = #LONG# ~; # the answer, can be REAL #
 MODE CPREAL = ~; # the answer, can be REAL #
 PROC cp fix value error = (#REF# CPARGS args)BOOL: ~;

PROVIDES:#
OP C = (CP~,CP~)CP~: ~ #
OP P = (CP~,CP~)CP~: ~ #

END COMMENT

MODE CPARGS = STRUCT(CHAR name, #REF# CPINT n,k);

PRIO C = 8, P = 8; # should be 7.5, a priority between *,/ and **,SHL,SHR etc #

I suspect there is a more reliable way of doing this using the Gamma Function approx #

OP P = (CPINT n, r)CPOUT: (

 IF n < r ORF r < 0 THEN IF NOT cp fix value error(CPARGS("P",n,r)) THEN stop FI FI;
 CPOUT out := 1;

basically nPk = (n-r+1)(n-r+2)...(n-2)(n-1)n = n!/(n-r)! #

 FOR i FROM n-r+1 TO n DO out *:= i OD;
 out

);

OP P = (CPREAL n, r)CPREAL: # 'ln gamma' requires GSL library #

 exp(ln gamma(n+1)-ln gamma(n-r+1));

basically nPk = (n-r+1)(n-r+2)...(n-2)(n-1)n = n!/(n-r)! #

COMMENT # alternate slower version # OP P = (CPREAL n, r)CPREAL: ( # alternate slower version #

 IF n < r ORF r < 0 THEN IF NOT cp fix value error(CPARGS("P",ENTIER n,ENTIER r)) THEN stop FI FI;
 CPREAL out := 1;

basically nPk = (n-r+1)(n-r+2)...(n-2)(n-1)n = n!/(n-r)! #

 CPREAL i := n-r+1;
 WHILE i <= n DO
   out*:= i;

a crude check for underflow #

   IF i = i + 1 THEN IF NOT cp fix value error(CPARGS("P",ENTIER n,ENTIER r)) THEN stop FI FI;
   i+:=1
 OD;
 out

); END COMMENT

basically C(n,r) = nCk = nPk/r! = n!/(n-r)!/r! #

OP C = (CPINT n, r)CPOUT: (

 IF n < r ORF r < 0 THEN IF NOT cp fix value error(("C",n,r)) THEN stop FI FI;
 CPINT largest = ( r > n - r | r | n - r );
 CPINT smallest = n - largest;
 CPOUT out := 1;
 INT smaller fact := 2;
 FOR larger fact FROM largest+1 TO n DO

try and prevent overflow, p.s. there must be a smarter way to do this #
Problems: loop stalls when 'smaller fact' is a largeish co prime #

   out *:= larger fact;
   WHILE smaller fact <= smallest ANDF out MOD smaller fact = 0 DO 
     out OVERAB smaller fact;
     smaller fact +:= 1
   OD
 OD;
 out # EXIT with: n P r OVER r P r #

);

OP C = (CPREAL n, CPREAL r)CPREAL: # 'ln gamma' requires GSL library #

 exp(ln gamma(n+1)-ln gamma(n-r+1)-ln gamma(r+1));

basically C(n,r) = nCk = nPk/r! = n!/(n-r)!/r! #

COMMENT # alternate slower version # OP C = (CPREAL n, REAL r)CPREAL: (

 IF n < r ORF r < 0 THEN IF NOT cp fix value error(("C",ENTIER n,ENTIER r)) THEN stop FI FI;
 CPREAL largest = ( r > n - r | r | n - r );
 CPREAL smallest = n - largest;
 CPREAL out := 1;
 REAL smaller fact := 2;
 REAL larger fact := largest+1;
 WHILE larger fact <= n DO # todo: check underflow here #

try and prevent overflow, p.s. there must be a smarter way to do this #

   out *:= larger fact;
   WHILE smaller fact <= smallest ANDF out > smaller fact DO 
     out /:= smaller fact;
     smaller fact +:= 1
   OD;
   larger fact +:= 1
 OD;
 out # EXIT with: n P r OVER r P r #

); END COMMENT

SKIP</lang>File: test_combinations_and_permutations.a68<lang algol68>#!/usr/bin/a68g --script #

-*- coding: utf-8 -*- #

CO REQUIRED by "prelude_combinations_and_permutations.a68" CO

 MODE CPINT = #LONG# INT;
 MODE CPOUT = #LONG# INT; # the answer, can be REAL #
 MODE CPREAL = REAL; # the answer, can be REAL #
 PROC cp fix value error = (#REF# CPARGS args)BOOL: (
   putf(stand error, ($"Value error: "g(0)gg(0)"arg out of range"l$, 
                        n OF args, name OF args, k OF args)); 
   FALSE # unfixable #
 );

PROVIDES:#
OP C = (CP~,CP~)CP~: ~ #
OP P = (CP~,CP~)CP~: ~ #

PR READ "prelude_combinations_and_permutations.a68" PR;

printf($"A sample of Permutations from 1 to 12:"l$); FOR i FROM 4 BY 1 TO 12 DO

 INT first = i - 2,
     second = i - ENTIER sqrt(i);
 printf(($g(0)" P "g(0)" = "g(0)$, i, first, i P first, $", "$));
 printf(($g(0)" P "g(0)" = "g(0)$, i, second, i P second, $l$))

OD;

printf($l"A sample of Combinations from 10 to 60:"l$); FOR i FROM 10 BY 10 TO 60 DO

 INT first = i - 2,
     second = i - ENTIER sqrt(i);
 printf(($"("g(0)" C "g(0)") = "g(0)$, i, first, i C first, $", "$));
 printf(($"("g(0)" C "g(0)") = "g(0)$, i, second, i C second, $l$))

OD;

printf($l"A sample of Permutations from 5 to 15000:"l$); FOR i FROM 5 BY 10 TO 150 DO

 REAL r = i,
      first = r - 2,
      second = r - ENTIER sqrt(r);
 printf(($g(0)" P "g(0)" = "g(-real width,real width-5,-1)$, r, first, r P first, $", "$));
 printf(($g(0)" P "g(0)" = "g(-real width,real width-5,-1)$, r, second, r P second, $l$))

OD;

printf($l"A sample of Combinations from 10 to 190:"l$); FOR i FROM 100 BY 100 TO 1000 DO

 REAL r = i,
      first = r - 2,
      second = r - ENTIER sqrt(r);
 printf(($"("g(0)" C "g(0)") = "g(0,1)$, r, first, r C first, $", "$));
 printf(($"("g(0)" C "g(0)") = "g(0,1)$, r, second, r C second, $l$))

OD</lang>Output:

A sample of Permutations from 1 to 12:
4 P 2 = 12, 4 P 2 = 12
5 P 3 = 60, 5 P 3 = 60
6 P 4 = 360, 6 P 4 = 360
7 P 5 = 2520, 7 P 5 = 2520
8 P 6 = 20160, 8 P 6 = 20160
9 P 7 = 181440, 9 P 6 = 60480
10 P 8 = 1814400, 10 P 7 = 604800
11 P 9 = 19958400, 11 P 8 = 6652800
12 P 10 = 239500800, 12 P 9 = 79833600

A sample of Combinations from 10 to 60:
(10 C 8) = 45, (10 C 7) = 120
(20 C 18) = 190, (20 C 16) = 4845
(30 C 28) = 435, (30 C 25) = 142506
(40 C 38) = 780, (40 C 34) = 3838380
(50 C 48) = 1225, (50 C 43) = 99884400
(60 C 58) = 1770, (60 C 53) = 386206920

A sample of Permutations from 5 to 15000:
5 P 3 =  6.0000000000e1, 5 P 3 =  6.0000000000e1
15 P 13 =  6.538371840e11, 15 P 12 =  2.179457280e11
25 P 23 =  7.755605022e24, 25 P 20 =  1.292600837e23
35 P 33 =  5.166573983e39, 35 P 30 =  8.610956639e37
45 P 43 =  5.981111043e55, 45 P 39 =  1.661419734e53
55 P 53 =  6.348201677e72, 55 P 48 =  2.519127650e69
65 P 63 =  4.123825296e90, 65 P 57 =  2.045548262e86
75 P 73 =  1.24045704e109, 75 P 67 =  6.15306072e104
85 P 83 =  1.40855206e128, 85 P 76 =  7.76318374e122
95 P 93 =  5.16498924e147, 95 P 86 =  2.84666515e142
105 P 103 =  5.40698379e167, 105 P 95 =  2.98003957e161
115 P 113 =  1.46254685e188, 115 P 105 =  8.06077407e181
125 P 123 =  9.41338588e208, 125 P 114 =  4.71650327e201
135 P 133 =  1.34523635e230, 135 P 124 =  6.74020139e222
145 P 143 =  4.02396303e251, 145 P 133 =  1.68014597e243

A sample of Combinations from 10 to 190:
(100 C 98) = 4950.0, (100 C 90) = 17310309456438.8
(200 C 198) = 19900.0, (200 C 186) = 1179791641436960000000.0
(300 C 298) = 44850.0, (300 C 283) = 2287708142022840000000000000.0
(400 C 398) = 79800.0, (400 C 380) = 2788360983670300000000000000000000.0
(500 C 498) = 124750.0, (500 C 478) = 132736424690773000000000000000000000000.0
(600 C 598) = 179700.0, (600 C 576) = 4791686682467800000000000000000000000000000.0
(700 C 698) = 244650.0, (700 C 674) = 145478651313640000000000000000000000000000000000.0
(800 C 798) = 319600.0, (800 C 772) = 3933526871034430000000000000000000000000000000000000.0
(900 C 898) = 404550.0, (900 C 870) = 98033481673646900000000000000000000000000000000000000000.0
(1000 C 998) = 499500.0, (1000 C 969) = 76023224077705100000000000000000000000000000000000000000000.0

C

Using big integers. GMP in fact has a factorial function which is quite possibly more efficient, though using it would make code longer. <lang c>#include <gmp.h>

void perm(mpz_t out, int n, int k) { mpz_set_ui(out, 1); k = n - k; while (n > k) mpz_mul_ui(out, out, n--); }

void comb(mpz_t out, int n, int k) { perm(out, n, k); while (k) mpz_divexact_ui(out, out, k--); }

int main(void) { mpz_t x; mpz_init(x);

perm(x, 1000, 969); gmp_printf("P(1000,969) = %Zd\n", x);

comb(x, 1000, 969); gmp_printf("C(1000,969) = %Zd\n", x); return 0; }</lang>

D

Translation of: Ruby

<lang d>import std.stdio, std.mathspecial, std.range, std.algorithm,

      std.bigint, std.conv;

enum permutation = (in uint n, in uint k) =>

   reduce!q{a * b}(1.BigInt, iota(n - k + 1, n + 1));

enum combination = (in uint n, in uint k) =>

   n.permutation(k) / reduce!q{a * b}(1.BigInt, iota(1, k + 1));

enum bigPermutation = (in uint n, in uint k) =>

   exp(logGamma(n + 1) - logGamma(n - k + 1));

enum bigCombination = (in uint n, in uint k) =>

   exp(logGamma(n + 1) - logGamma(n - k + 1) - logGamma(k + 1));

void main() {

   12.permutation(9).writeln;
   12.bigPermutation(9).writeln;
   60.combination(53).writeln;
   145.bigPermutation(133).writeln;
   900.bigCombination(450).writeln;
   1_000.bigCombination(969).writeln;
   15_000.bigPermutation(73).writeln;
   15_000.bigPermutation(1185).writeln;
   writefln("%(%s\\\n%)", 15_000.permutation(74).text.chunks(50));

}</lang>

Output:

79833600
7.98336e+07
386206920
1.68015e+243
2.24747e+269
7.60232e+58
6.00414e+304
6.31335e+4927
89623761385296782623991723856543314935307441602519\
77843015933352436993580407381279508723841971598849\
05490054194835376498534786047382445592358843238688\
90331846707057518455295399761517897302775271453951\
38931598154729489879215876713997904109589031888166\
84444202526779550201576117111844818124800000000000\
000000000

J

It looks like this task wants a count of the available combinations or permutations (given a set of 3 things, there are three distinct combinations of 2 of them) rather than a representation of the available combinations or permutations (given a set of three things, the distinct combinations of 2 of them may be identified by <0,1>, <0,2> and <1,2>)).

Also, this task allows a language to show off its abilities to support floating point numbers outside the usual range of 64 bit IEEE floating point numbers. We'll neglect that part.

Implementation:

C is a primitive, but we will give it a name for this task. P could be defined using %&!~ but that approach would overflow for large arguments even when the result would have been representable.

Example use:

<lang J> P table 1+i.12 ┌──┬─────────────────────────────────────────────────────────────┐ │P │1 2 3 4 5 6 7 8 9 10 11 12│ ├──┼─────────────────────────────────────────────────────────────┤ │ 1│1 2 6 24 120 720 5040 40320 362880 3628800 39916800 479001600│ │ 2│0 1 3 12 60 360 2520 20160 181440 1814400 19958400 239500800│ │ 3│0 0 1 4 20 120 840 6720 60480 604800 6652800 79833600│ │ 4│0 0 0 1 5 30 210 1680 15120 151200 1663200 19958400│ │ 5│0 0 0 0 1 6 42 336 3024 30240 332640 3991680│ │ 6│0 0 0 0 0 1 7 56 504 5040 55440 665280│ │ 7│0 0 0 0 0 0 1 8 72 720 7920 95040│ │ 8│0 0 0 0 0 0 0 1 9 90 990 11880│ │ 9│0 0 0 0 0 0 0 0 1 10 110 1320│ │10│0 0 0 0 0 0 0 0 0 1 11 132│ │11│0 0 0 0 0 0 0 0 0 0 1 12│ │12│0 0 0 0 0 0 0 0 0 0 0 1│ └──┴─────────────────────────────────────────────────────────────┘

  C table 10+10*i.6x

┌──┬─────────────────────────────────────────────────────────────────┐ │C │10 20 30 40 50 60│ ├──┼─────────────────────────────────────────────────────────────────┤ │10│ 1 184756 30045015 847660528 10272278170 75394027566│ │20│ 0 1 30045015 137846528820 47129212243960 4191844505805495│ │30│ 0 0 1 847660528 47129212243960 118264581564861424│ │40│ 0 0 0 1 10272278170 4191844505805495│ │50│ 0 0 0 0 1 75394027566│ │60│ 0 0 0 0 0 1│ └──┴─────────────────────────────────────────────────────────────────┘

  5 P 100

7.77718e155

  100 P 200

8.45055e216

  300 P 400

2.09224e254

  700 P 800

3.18349e287

  5 C 100

75287520

  100 C 200

9.05485e58

  300 C 400

2.24185e96

  700 C 800

3.41114e129</lang>

МК-61/52

Input: x ^ n ^ k В/О С/П, where x = 8 for permutations; 20 for permutations with repetitions; 26 for combinations; 44 for combinations with repetitions.

Printing of test cases is performed incrementally, which is associated with the characteristics of the device output.

Perl

Although perl can handle arbitrarily large numbers using Math::BigInt and Math::BigFloat, it's native integers and floats are limited to what the computer's native types can handle.

As with the perl6 code, some special handling was done for those values which would have overflowed the native floating point type.

<lang perl>use strict; use warnings;

showoff( "Permutations", \&P, "P", 1 .. 12 ); showoff( "Combinations", \&C, "C", map $_*10, 1..6 ); showoff( "Permutations", \&P_big, "P", 5, 50, 500, 1000, 5000, 15000 ); showoff( "Combinations", \&C_big, "C", map $_*100, 1..10 );

sub showoff { my ($text, $code, $fname, @n) = @_; print "\nA sample of $text from $n[0] to $n[-1]\n"; for my $n ( @n ) { my $k = int( $n / 3 ); print $n, " $fname $k = ", $code->($n, $k), "\n"; } }

sub P { my ($n, $k) = @_; my $x = 1; $x *= $_ for $n - $k + 1 .. $n ; $x; }

sub P_big { my ($n, $k) = @_; my $x = 0; $x += log($_) for $n - $k + 1 .. $n ; eshow($x); }

sub C { my ($n, $k) = @_; my $x = 1; $x *= ($n - $_ + 1) / $_ for 1 .. $k; $x; }

sub C_big { my ($n, $k) = @_; my $x = 0; $x += log($n - $_ + 1) - log($_) for 1 .. $k; exp($x); }

sub eshow { my ($x) = @_; my $e = int( $x / log(10) ); sprintf "%.8Fe%+d", exp($x - $e * log(10)), $e; } </lang>

Since the output is almost the same as perl6's, and this is only a Draft RosettaCode task, I'm not going to bother including the output of the program.

Perl 6

Perl 6 can't compute arbitrary large floating point values, thus we will use logarithms, as is often needed when dealing with combinations. We'll also use a Stirling method to approximate $\ln(n!)$ :

$\ln n!\approx {\frac {1}{2}}\ln(2\pi n)+n\ln \left({\frac {n}{e}}+{\frac {1}{12en}}\right)$

Notice that Perl6 can process arbitrary long integers, though. So it's not clear whether using floating points is useful in this case.

<lang perl6>multi P($n, $k) { [*] $n - $k + 1 .. $n } multi C($n, $k) { P($n, $k) / [*] 1 .. $k }

sub lstirling(\n) {

   n < 10 ?? lstirling(n+1) - log(n+1) !!
   .5*log(2*pi*n)+ n*log(n/e+1/(12*e*n))

}

role Logarithm {

   method gist {

my $e = (self/10.log).Int; sprintf "%.8fE%+d", exp(self - $e*10.log), $e;

} multi P($n, $k, :$float!) {

   (lstirling($n) - lstirling($n -$k))
   but Logarithm

} multi C($n, $k, :$float!) {

   (lstirling($n) - lstirling($n -$k) - lstirling($k))
   but Logarithm

}

say "Exact results:"; for 1..12 -> $n {

   my $p = $n div 3;
   say "P($n, $p) = ", P($n, $p);

}

for 10, 20 ... 60 -> $n {

   my $p = $n div 3;
   say "C($n, $p) = ", C($n, $p);

}

say; say "Floating point approximations:"; for 5, 50, 500, 1000, 5000, 15000 -> $n {

   my $p = $n div 3;
   say "P($n, $p) = ", P($n, $p, :float);

}

for 100, 200 ... 1000 -> $n {

   my $p = $n div 3;
   say "C($n, $p) = ", C($n, $p, :float);

}</lang>

Output:

Exact results:
P(1, 0) = 1
P(2, 0) = 1
P(3, 1) = 3
P(4, 1) = 4
P(5, 1) = 5
P(6, 2) = 30
P(7, 2) = 42
P(8, 2) = 56
P(9, 3) = 504
P(10, 3) = 720
P(11, 3) = 990
P(12, 4) = 11880
C(10, 3) = 120
C(20, 6) = 38760
C(30, 10) = 30045015
C(40, 13) = 12033222880
C(50, 16) = 4923689695575
C(60, 20) = 4191844505805495

Floating point approximations:
P(5, 1) = 5.00000000E+0
P(50, 16) = 1.03017326E+26
P(500, 166) = 3.53487492E+434
P(1000, 333) = 5.96932629E+971
P(5000, 1666) = 6.85674576E+6025
P(15000, 5000) = 9.64985399E+20469
C(100, 33) = 2.94692433E+26
C(200, 66) = 7.26975256E+53
C(300, 100) = 4.15825147E+81
C(400, 133) = 1.25794868E+109
C(500, 166) = 3.92602839E+136
C(600, 200) = 2.50601778E+164
C(700, 233) = 8.10320356E+191
C(800, 266) = 2.64562336E+219
C(900, 300) = 1.74335637E+247
C(1000, 333) = 5.77613455E+274

Python

Using library scipy <lang python>from __future__ import print_function

from scipy.misc import factorial as fact from scipy.misc import comb

def perm(N, k, exact=0):

   return comb(N, k, exact) * fact(k, exact)

exact=True print('Sample Perms 1..12') for N in range(1, 13):

   k = max(N-2, 1)
   print('%iP%i =' % (N, k), perm(N, k, exact), end=', ' if N % 5 else '\n')

print('\n\nSample Combs 10..60') for N in range(10, 61, 10):

   k = N-2
   print('%iC%i =' % (N, k), comb(N, k, exact), end=', ' if N % 50 else '\n')

exact=False print('\n\nSample Perms 5..1500 Using FP approximations') for N in [5, 15, 150, 1500, 15000]:

   k = N-2
   print('%iP%i =' % (N, k), perm(N, k, exact))

print('\nSample Combs 100..1000 Using FP approximations') for N in range(100, 1001, 100):

   k = N-2
   print('%iC%i =' % (N, k), comb(N, k, exact))</lang>

Output:

Sample Perms 1..12
1P1 = 1, 2P1 = 2, 3P1 = 3, 4P2 = 12, 5P3 = 60
6P4 = 360, 7P5 = 2520, 8P6 = 20160, 9P7 = 181440, 10P8 = 1814400
11P9 = 19958400, 12P10 = 239500800, 

Sample Combs 10..60
10C8 = 45, 20C18 = 190, 30C28 = 435, 40C38 = 780, 50C48 = 1225
60C58 = 1770, 

Sample Perms 5..1500 Using FP approximations
5P3 = 60.0
15P13 = 653837184000.0
150P148 = 2.85669197822e+262
1500P1498 = inf
15000P14998 = inf

Sample Combs 100..1000 Using FP approximations
100C98 = 4950.0
200C198 = 19900.0
300C298 = 44850.0
400C398 = 79800.0
500C498 = 124750.0
600C598 = 179700.0
700C698 = 244650.0
800C798 = 319600.0
900C898 = 404550.0
1000C998 = 499500.0

Racket

Racket's "math" library has two functions that compute nCk and nPk. They work only on integers, but since Racket supports unlimited integers there is no need for a floating point estimate:

lang racket

(require math) (define C binomial) (define P permutations)

(C 1000 10) ; -> 263409560461970212832400 (P 1000 10) ; -> 955860613004397508326213120000 </lang>

(I'll spare this page from yet another big listing of samples...)

REXX

The hard part of this REXX program was coding the DO loops for the various ranges. <lang rexx>/*REXX program to compute a sampling of combinations and permutations. */ numeric digits 100 /*use hundred digits of precision*/

   do j=1  to  12                     /*show permutations from 1──► 12 */
   _=;  do k=1  to  j                 /*step through all J permutations*/
        _=_ 'P('j","k')='perm(j,k)" " /*add an extra blank between #s. */
        end   /*k*/
   say strip(_)                       /*show a horizontal line of PERMs*/
   end           /*j*/

say

   do j=10  to  60  by 10             /*show some combinations 10──► 60*/
   _=;  do k=1  to  j  by j%5         /*step through some combinations.*/
        _=_ 'C('j","k')='comb(j,k)" " /*add an extra blank between #s. */
        end   /*k*/
   say strip(_)                       /*show a horizontal line of COMBs*/
   end           /*j*/

say numeric digits 20 /*force floating point for big #s*/

   do j=5  to 15000  by 1000          /*show a few permutations, big #s*/
   _=;  do k=1  to  j  by j%10  for 5 /*go through some J permutations.*/
        _=_ 'P('j","k')='perm(j,k)" " /*add an extra blank between #s. */
        end   /*k*/
   say strip(_)                       /*show a horizontal line of PERMs*/
   end           /*j*/

say

   do j=100  to 1000  by 100          /*show a few combinations, big #s*/
   _=;  do k=1  to  j  by j%5         /*step through some combinations.*/
        _=_ 'C('j","k')='comb(j,k)" " /*add an extra blank between #s. */
        end   /*k*/
   say strip(_)                       /*show a horizontal line of COMBs*/
   end           /*j*/

exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────COMB subroutine─────────────────────*/ comb: procedure; parse arg x,y /*args: X things, Y at-a-time.*/ if y>x then return 0 /*oops-say, to big a chunk. */ if x=y then return 1 /* X things same as chunk size. */ if x-y<y then y=x-y /*switch things around for speed.*/ call .cmbPrm /*call sub to do heavy lifting. */ return _/!(y) /*perform one last division. */ /*──────────────────────────────────PERM subroutine─────────────────────*/ perm: procedure; parse arg x,y; call .cmbPrm; return _ /*──────────────────────────────────.CMBPRM sugroutine──────────────────*/ .cmbPrm: _=1; do j=x-y+1 to x; _=_*j; end; return _ /*──────────────────────────────────! subroutine────────────────────────*/ !: procedure; parse arg x; !=1; do j=2 to x; !=!*j; end; return !</lang> output

P(1,1)=1
P(2,1)=2  P(2,2)=2
P(3,1)=3  P(3,2)=6  P(3,3)=6
P(4,1)=4  P(4,2)=12  P(4,3)=24  P(4,4)=24
P(5,1)=5  P(5,2)=20  P(5,3)=60  P(5,4)=120  P(5,5)=120
P(6,1)=6  P(6,2)=30  P(6,3)=120  P(6,4)=360  P(6,5)=720  P(6,6)=720
P(7,1)=7  P(7,2)=42  P(7,3)=210  P(7,4)=840  P(7,5)=2520  P(7,6)=5040  P(7,7)=5040
P(8,1)=8  P(8,2)=56  P(8,3)=336  P(8,4)=1680  P(8,5)=6720  P(8,6)=20160  P(8,7)=40320  P(8,8)=40320
P(9,1)=9  P(9,2)=72  P(9,3)=504  P(9,4)=3024  P(9,5)=15120  P(9,6)=60480  P(9,7)=181440  P(9,8)=362880  P(9,9)=362880
P(10,1)=10  P(10,2)=90  P(10,3)=720  P(10,4)=5040  P(10,5)=30240  P(10,6)=151200  P(10,7)=604800  P(10,8)=1814400  P(10,9)=3628800  P(10,10)=3628800
P(11,1)=11  P(11,2)=110  P(11,3)=990  P(11,4)=7920  P(11,5)=55440  P(11,6)=332640  P(11,7)=1663200  P(11,8)=6652800  P(11,9)=19958400  P(11,10)=39916800  P(11,11)=39916800
P(12,1)=12  P(12,2)=132  P(12,3)=1320  P(12,4)=11880  P(12,5)=95040  P(12,6)=665280  P(12,7)=3991680  P(12,8)=19958400  P(12,9)=79833600  P(12,10)=239500800  P(12,11)=479001600  P(12,12)=479001600

C(10,1)=10  C(10,3)=120  C(10,5)=252  C(10,7)=120  C(10,9)=10
C(20,1)=20  C(20,5)=15504  C(20,9)=167960  C(20,13)=77520  C(20,17)=1140
C(30,1)=30  C(30,7)=2035800  C(30,13)=119759850  C(30,19)=54627300  C(30,25)=142506
C(40,1)=40  C(40,9)=273438880  C(40,17)=88732378800  C(40,25)=40225345056  C(40,33)=18643560
C(50,1)=50  C(50,11)=37353738800  C(50,21)=67327446062800  C(50,31)=30405943383200  C(50,41)=2505433700
C(60,1)=60  C(60,13)=5166863427600  C(60,25)=51915437974328292  C(60,37)=23385332420868600  C(60,49)=342700125300

P(5,1)=5  P(5,1)=5  P(5,1)=5  P(5,1)=5  P(5,1)=5
P(1005,1)=1005  P(1005,101)=9.1176524923776877363E+300  P(1005,201)=1.2772738260896333926E+594  P(1005,301)=6.9244021230613662196E+881  P(1005,401)=2.4492580742838357278E+1163
P(2005,1)=2005  P(2005,201)=1.6533543480914610058E+659  P(2005,401)=3.0753126526205309249E+1305  P(2005,601)=7.9852540678709597130E+1940  P(2005,801)=8.0516979630356802995E+2563
P(3005,1)=3005  P(3005,301)=1.1935689764209015622E+1040  P(3005,601)=1.5619600532077469150E+2062  P(3005,901)=1.0291767405881430479E+3068  P(3005,1201)=1.5669988662999668720E+4055
P(4005,1)=4005  P(4005,401)=4.3808609526948101266E+1435  P(4005,801)=2.3060742016678396933E+2848  P(4005,1201)=2.2044072986703009755E+4239  P(4005,1601)=2.8973897505902543204E+5605
P(5005,1)=5005  P(5005,501)=1.4180262672357809801E+1842  P(5005,1001)=2.8239356430641573722E+3656  P(5005,1501)=3.6832518654277810594E+5443  P(5005,2001)=3.9303728189857603162E+7199
P(6005,1)=6005  P(6005,601)=2.4482219222979658097E+2257  P(6005,1201)=1.0247320583108167487E+4482  P(6005,1801)=1.0131515875595211375E+6674  P(6005,2401)=4.8762853043004294329E+8828
P(7005,1)=7005  P(7005,701)=7.6900396347210828241E+2679  P(7005,1401)=1.2659048848290269952E+5322  P(7005,2101)=1.7753130713788487191E+7926  P(7005,2801)=7.2114365718218704695E+10486
P(8005,1)=8005  P(8005,801)=3.9769062582658855959E+3108  P(8005,1601)=4.3272401446508603181E+6174  P(8005,2401)=1.4466844005282015778E+9197  P(8005,3201)=8.3354759867215982278E+12169
P(9005,1)=9005  P(9005,901)=5.6135384805755901099E+3542  P(9005,1801)=1.1194389175552115248E+7038  P(9005,2701)=2.4737530806300682806E+10484  P(9005,3601)=5.6056491332873455398E+13874
P(10005,1)=10005  P(10005,1001)=5.3580936683833889197E+3981  P(10005,2001)=1.3407350644082770778E+7911  P(10005,3001)=1.3407953461588097193E+11786  P(10005,4001)=8.1811569565040010437E+15598
P(11005,1)=11005  P(11005,1101)=1.1340564277915775963E+4425  P(11005,2201)=7.9753039151558717610E+8792  P(11005,3301)=8.0842724022079710248E+13100  P(11005,4401)=2.9749926937463675736E+17340
P(12005,1)=12005  P(12005,1201)=2.1391703159775094656E+4872  P(12005,2401)=3.7994859265471124812E+9682  P(12005,3601)=3.5081603331307865953E+14427  P(12005,4801)=6.9968644993020337359E+19097
P(13005,1)=13005  P(13005,1301)=1.6832659142209713962E+5323  P(13005,2601)=3.1719005022749408296E+10579  P(13005,3901)=1.1206120643200361213E+15765  P(13005,5201)=5.0883658993886790178E+20869
P(14005,1)=14005  P(14005,1401)=2.9074578200382556975E+5777  P(14005,2801)=1.2835011192416517281E+11483  P(14005,4201)=3.8312600202917546343E+17112  P(14005,5601)=8.7456467698123057261E+22654

C(100,1)=100  C(100,21)=2.0418414110621321255E+21  C(100,41)=2.0116440213369968048E+28  C(100,61)=9.0139240300346304925E+27  C(100,81)=1.3234157293921226741E+20
C(200,1)=200  C(200,41)=8.0006165666286406037E+42  C(200,81)=2.4404128184470558197E+57  C(200,121)=1.0891098528606695394E+57  C(200,161)=5.0935602365182339252E+41
C(300,1)=300  C(300,61)=3.5574671252567510894E+64  C(300,121)=3.3878557197772169409E+86  C(300,181)=1.5098730832156289128E+86  C(300,241)=2.2510943427454545270E+63
C(400,1)=400  C(400,81)=1.6703771503944415835E+86  C(400,161)=4.9770797199515347150E+115  C(400,241)=2.2166247162163128334E+115  C(400,321)=1.0537425948749981954E+85
C(500,1)=500  C(500,101)=8.0859177660929770887E+107  C(500,201)=7.5447012685604958486E+144  C(500,301)=3.3587706644089915087E+144  C(500,401)=5.0915068227892187414E+106
C(600,1)=600  C(600,121)=3.9913554739811382543E+129  C(600,241)=1.1667430218545073615E+174  C(600,361)=5.1927067085306791157E+173  C(600,481)=2.5101559893540422495E+128
C(700,1)=700  C(700,141)=1.9971304197729039382E+151  C(700,281)=1.8294137562513979560E+203  C(700,421)=8.1403842518866639077E+202  C(700,561)=1.2548814134936695871E+150
C(800,1)=800  C(800,161)=1.0093242166331589874E+173  C(800,321)=2.8977104736455704539E+232  C(800,481)=1.2892100651978213666E+232  C(800,641)=6.3378002682503352943E+171
C(900,1)=900  C(900,181)=5.1402207737939392620E+194  C(900,361)=4.6256239559539226532E+261  C(900,541)=2.0577328996911473555E+261  C(900,721)=3.2260054093505652100E+193
C(1000,1)=1000  C(1000,201)=2.6336937554862900107E+216  C(1000,401)=7.4293352412781479131E+290  C(1000,601)=3.3046738011675400053E+290  C(1000,801)=1.6522236106515115238E+215

Ruby

Float calculation as Tcl. <lang ruby>include Math

class Integer

 def permutation(k)
   (self-k+1 .. self).inject( :*)
 end

 def combination(k)
   self.permutation(k) / (1 .. k).inject( :*)
 end

 def big_permutation(k)
   exp( lgamma_plus(self) - lgamma_plus(self -k)) 
 end

 def big_combination(k)
   exp( lgamma_plus(self) - lgamma_plus(self - k) - lgamma_plus(k))
 end

 private
 def lgamma_plus(n)
   lgamma(n+1)[0]  #lgamma is the natural log of gamma
 end

end

p 12.permutation(9) #=> 79833600 p 12.big_permutation(9) #=> 79833600.00000021 p 60.combination(53) #=> 386206920 p 145.big_permutation(133) #=> 1.6801459655817956e+243 p 900.big_combination(450) #=> 2.247471882064647e+269 p 1000.big_combination(969) #=> 7.602322407770517e+58 p 15000.big_permutation(73) #=> 6.004137561717704e+304

That's about the maximum of Float:

p 15000.big_permutation(74) #=> Infinity

Integer has no maximum:

p 15000.permutation(74) #=> 896237613852967826239917238565433149353074416025197784301593335243699358040738127950872384197159884905490054194835376498534786047382445592358843238688903318467070575184552953997615178973027752714539513893159815472948987921587671399790410958903188816684444202526779550201576117111844818124800000000000000000000 </lang>

Tcl

Tcl doesn't allow the definition of new infix operators, so we define $P$ and $C$ as ordinary functions. There are no problems with loss of significance though: Tcl has supported arbitrary precision integer arithmetic since 8.5.

Library: Tcllib (Package: math)

<lang tcl># Exact integer versions proc tcl::mathfunc::P {n k} {

   set t 1
   for {set i $n} {$i > $n-$k} {incr i -1} {

set t [expr {$t * $i}]

   }
   return $t

} proc tcl::mathfunc::C {n k} {

   set t [P $n $k]
   for {set i $k} {$i > 1} {incr i -1} {

set t [expr {$t / $i}]

   }
   return $t

}

Floating point versions using the Gamma function

package require math proc tcl::mathfunc::lnGamma n {math::ln_Gamma $n} proc tcl::mathfunc::fP {n k} {

   expr {exp(lnGamma($n+1) - lnGamma($n-$k+1))}

} proc tcl::mathfunc::fC {n k} {

   expr {exp(lnGamma($n+1) - lnGamma($n-$k+1) - lnGamma($k+1))}

}</lang> Demonstrating: <lang tcl># Using the exact integer versions puts "A sample of Permutations from 1 to 12:" for {set i 4} {$i <= 12} {incr i} {

   set ii [expr {$i - 2}]
   set iii [expr {$i - int(sqrt($i))}]
   puts "$i P $ii = [expr {P($i,$ii)}], $i P $iii = [expr {P($i,$iii)}]"

} puts "A sample of Combinations from 10 to 60:" for {set i 10} {$i <= 60} {incr i 10} {

   set ii [expr {$i - 2}]
   set iii [expr {$i - int(sqrt($i))}]
   puts "$i C $ii = [expr {C($i,$ii)}], $i C $iii = [expr {C($i,$iii)}]"

}

Using the approximate floating point versions

puts "A sample of Permutations from 5 to 15000:" for {set i 5} {$i <= 150} {incr i 10} {

   set ii [expr {$i - 2}]
   set iii [expr {$i - int(sqrt($i))}]
   puts "$i P $ii = [expr {fP($i,$ii)}], $i P $iii = [expr {fP($i,$iii)}]"

} puts "A sample of Combinations from 100 to 1000:" for {set i 100} {$i <= 1000} {incr i 100} {

   set ii [expr {$i - 2}]
   set iii [expr {$i - int(sqrt($i))}]
   puts "$i C $ii = [expr {fC($i,$ii)}], $i C $iii = [expr {fC($i,$iii)}]"

}</lang>

Output:

A sample of Permutations from 1 to 12:
4 P 2 = 12, 4 P 2 = 12
5 P 3 = 60, 5 P 3 = 60
6 P 4 = 360, 6 P 4 = 360
7 P 5 = 2520, 7 P 5 = 2520
8 P 6 = 20160, 8 P 6 = 20160
9 P 7 = 181440, 9 P 6 = 60480
10 P 8 = 1814400, 10 P 7 = 604800
11 P 9 = 19958400, 11 P 8 = 6652800
12 P 10 = 239500800, 12 P 9 = 79833600
A sample of Combinations from 10 to 60:
10 C 8 = 45, 10 C 7 = 120
20 C 18 = 190, 20 C 16 = 4845
30 C 28 = 435, 30 C 25 = 142506
40 C 38 = 780, 40 C 34 = 3838380
50 C 48 = 1225, 50 C 43 = 99884400
60 C 58 = 1770, 60 C 53 = 386206920
A sample of Permutations from 5 to 15000:
5 P 3 = 59.9999999964319, 5 P 3 = 59.9999999964319
15 P 13 = 653837183936.7548, 15 P 12 = 217945727984.54794
25 P 23 = 7.755605021026223e+24, 25 P 20 = 1.2926008369145724e+23
35 P 33 = 5.166573982873315e+39, 35 P 30 = 8.610956638634269e+37
45 P 43 = 5.981111043018166e+55, 45 P 39 = 1.6614197342883882e+53
55 P 53 = 6.348201676661335e+72, 55 P 48 = 2.5191276496660396e+69
65 P 63 = 4.123825295988996e+90, 65 P 57 = 2.0455482620718488e+86
75 P 73 = 1.2404570405684596e+109, 75 P 67 = 6.153060717624475e+104
85 P 83 = 1.4085520572027225e+128, 85 P 76 = 7.763183737477006e+122
95 P 93 = 5.164989244208789e+147, 95 P 86 = 2.846665148075141e+142
105 P 103 = 5.406983791334563e+167, 105 P 95 = 2.980039567808848e+161
115 P 113 = 1.462546846791721e+188, 115 P 105 = 8.060774068156828e+181
125 P 123 = 9.413385884788385e+208, 125 P 114 = 4.716503269639238e+201
135 P 133 = 1.345236353714729e+230, 135 P 124 = 6.74020138809567e+222
145 P 143 = 4.0239630289197437e+251, 145 P 133 = 1.6801459658196038e+243
A sample of Combinations from 100 to 1000:
100 C 98 = 4950.000000564707, 100 C 90 = 17310309460118.861
200 C 198 = 19900.000002250566, 200 C 186 = 1.1797916416885855e+21
300 C 298 = 44850.00000506082, 300 C 283 = 2.287708142503998e+27
400 C 398 = 79800.00000901309, 400 C 380 = 2.788360984244711e+33
500 C 498 = 124750.00001405331, 500 C 478 = 1.327364247175741e+38
600 C 598 = 179700.00002031153, 600 C 576 = 4.7916866834178515e+42
700 C 698 = 244650.00002750417, 700 C 674 = 1.454786513417567e+47
800 C 798 = 319600.0000360682, 800 C 772 = 3.933526871778561e+51
900 C 898 = 404550.0000452471, 900 C 870 = 9.803348169192494e+55
1000 C 998 = 499500.0000564987, 1000 C 969 = 7.602322409167201e+58

It should be noted that for large values, it can be much faster to use the floating point version (at a cost of losing significance). In particular expr C(1000,500) takes approximately 1000 times longer to compute than expr fC(1000,500)