Rebol [
title: "Rosetta code: Combinations and permutations"
file: %Combinations_and_permutations.r3
url: https://rosettacode.org/wiki/Combinations_and_permutations
]
perm: function [
"Falling factorial: x * (x-1) * ... * (x-y+1), i.e. x! / (x-y)!"
x [integer!] y [integer!]
][
;; Computed in floating point so large intermediate products
;; lose precision rather than overflowing an integer!
z: 1.0
for i (x - y + 1) x 1 [z: z * i]
try/with [to integer! z][z]
]
fact: function ["Factorial of x, in floating point" x [integer!]] [
z: 1.0
for i 2 x 1 [z: z * i]
]
comb: function ["Binomial coefficient C(x,y)" x [integer!] y [integer!]] [
if (x - y) < y [y: x - y]
c: (perm x y) / (fact y)
try/with [to integer! c][c]
]
lstirling: function/with [
"Natural log of n! via Stirling's series; recurses up to n >= 10 first for better accuracy there"
n [integer!]
][
if n < 10 [return (lstirling (n + 1)) - (log-e (n + 1))]
(0.5 * log-e (2 * pi * n)) + (n * log-e ((n / e) + (1 / (12 * e * n))))
][ e: exp 1 ]
tolog: function [
"Convert a natural-log value v into 'mantissa e exponent' scientific notation, e.g. '1.234e567'"
v [number!]
][
h: to integer! (v / (log-e 10))
mantissa: exp (v - (h * (log-e 10)))
ajoin [mantissa "e" h]
]
permf: function [
"Scientific-notation approximation of P(n,k) for large n, via Stirling's series"
n [integer!] k [integer!]
][
tolog ((lstirling n) - (lstirling (n - k)))
]
combf: function [
"Scientific-notation approximation of C(n,k) for large n, via Stirling's series"
n [integer!] k [integer!]
][
tolog ((lstirling n) - (lstirling (n - k)) - (lstirling k))
]
print as-yellow "=> Exact results:"
for n 1 12 1 [
p: n // 3
print ajoin ["P(" n "," p ") = " as-green perm n p]
]
;; double has 53 bits for integer
for n 10 50 10 [
p: n // 3
print ajoin ["C(" n "," p ") = " as-green comb n p]
]
print ""
print as-yellow "=> Floating point approximations:"
foreach n [5 50 500 1000 5000 15000] [
p: n // 3
print ajoin ["P(" n "," p ") = " as-green permf n p]
]
for n 100 1000 100 [
p: n // 3
print ajoin ["C(" n "," p ") = " as-green combf n p]
]
=> 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
=> Floating point approximations:
P(5,1) = 5.0e0
P(50,16) = 1.03017326223237e26
P(500,166) = 3.53487492178935e434
P(1000,333) = 5.96932628851383e971
P(5000,1666) = 6.85674575725791e6025
P(15000,5000) = 9.64985398884416e20469
C(100,33) = 2.94692433197593e26
C(200,66) = 7.26975256420019e53
C(300,100) = 4.15825146640008e81
C(400,133) = 1.25794868458534e109
C(500,166) = 3.92602838684241e136
C(600,200) = 2.50601778345756e164
C(700,233) = 8.1032035638234e191
C(800,266) = 2.64562336278924e219
C(900,300) = 1.74335637334568e247
C(1000,333) = 5.77613455326658e274