BinarySplittingFactorial.frink

Download or view BinarySplittingFactorial.frink in plain text format


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


Download or view BinarySplittingFactorial.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 20136 days, 3 hours, 49 minutes ago.