BinarySplittingFactorial.frink

View or download 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]
*/


View or download 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 17480 days, 13 hours, 38 minutes ago.