/** 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] */