/** This demonstrates using a binary splitting algorithm to calculate very
large factorials. This can be faster than using the naive calculation.
Frink caches factorials up to a certain size (currently the largest is
10000!) and returns them directly from a cache. When building this cache,
all factorials smaller than that factorial are calculated and stored.
http://www.oocities.org/hjsmithh/Numbers/BinSplit.html
This file is no longer necessary as Frink now uses a binary splitting
algorithm to calculate giant factorials that are too big to fit into its
cache. Its algorithm, though, follows this very closely.
This function actually calculates a! / b! so it can be used to calculate
ratios of factorials, or the factorial n! by calling p[n, 0]
*/
p[a,b] :=
{
diff = a-b
// println["$a $b $diff"]
if diff > 3
{
m = (a+b) div 2
return p[a, m] * p[m,b]
}
if diff == 0
return 1
if diff == 1
return a
// The following cases are not strictly necessary but they reduce the
// amount of recursion necessary.
if diff == 2
return a(a-1)
if diff == 3
return a(a-1)(a-2)
return 0
}
// Calculate the factorial n! using a binary splitting algorithm. p[a,b]
// is equal to a! / b! so we can use it to calculate a! by setting b to 0.
factorialSplit[n] := p[n,0]
// Attempt to calculate binomial coefficents, that is
// m! / ((m-n)! n!)
// Which we can maybe do faster using binary spliting because p[a,b] calculates
// a! / b! so we've done a lot of dividing out already (by never dividing at
// all; there are just multiplications in calculating a! / b! )
//
// This implementation is no longer necessary. As of the 2017-04-04 release
// of Frink, binomial coefficients are calculated using binary splitting.
// However, Frink's implementation follows this closely.
binomialSplit[m,n] :=
{
if m < n
return 0
// Since binomial[m, n] == binomial[m, m-n], use the one that will give
// the larger divisor b! and thus the smaller initial value that needs to
// be divided later.
if n < (m-n)
n = m - n
inter = p[m,n]
// println["Inter is $inter"]
divisor = (m-n)!
// println["Dividing by " + (m-n) + "! = $divisor"]
return inter / divisor
}
/* Timing and tests */
/*
bss = 0 s
bs = 0 s
for i = 1 to 20000
{
m = random[0,15000]
n = random[0,m]
s = now[]
a = binomialSplit[m,n]
e = now[]
bss = bss + (e-s)
s = now[]
b = binomial[m,n]
e = now[]
bs = bs + (e-s)
if b != m!/((m-n)! n!)
println["$m $n"]
}
println["Time in binomialSplit: " + bss]
println["Time in binomial: " + bs]
*/