How I could use tail call optimisation on this combination function?

Question

My exercise, that you can see here, says that I need to implement a recursive version of C(n, k).

That's my solution:

module LE1.Recursao.Combinacao where

combina :: Integral a => a -> a -> a
combina n k | k == 0         = 1
            | n == k         = 1
combina n k | k > 0 && n > k = (combina (n - 1) k) + (combina (n - 1) (k - 1))
combina _ _                  = 0

Now I want to create a tail-recursive version of this function, so I don't get stack overflows for large numbers and also calculate the combination quicker!

I'm new to tail call optimisation, but I did that in Elixir for the Fibonacci series:

defmodule Fibonacci do
  def run(n) when n < 0, do: :error
  def run(n), do: run n, 1, 0
  def run(0, _, result), do: result
  def run(n, next, result), do: run n - 1, next + result, next
end

I understand this Fibonacci code and I think that the combination algorithm isn't too different, but I don't know how to start.

Answer 1

Daniel Wagner writes

combina n k = product [n-k+1 .. n] `div` product [1 .. k]

This can get rather inefficient for large n and medium-size k; the multiplications get huge. How might we keep them smaller? We can write a simple recurrence:

combina :: Integer -> Integer -> Integer
combina n k
  | k > n || k < 0 = 0
  | otherwise = combina' n k'
  where
    -- C(n,k) and C(n,n-k) are the same,
    -- so we choose the cheaper one.
    k' = min k (n-k)

-- Assumes 0 <= k <= n
combina' _n 0 = 1
combina' n k
  = -- as above
  product [n-k+1 .. n] `div` product [1 .. k]
  = -- expanding a bit
  (n-k+1) * product [n-k+2 .. n] `div` (product [1 .. k-1] * k)
  = -- Rearranging, and temporarily using real division
  ((n-k+1)/k) * (product [n-(k-1)+1 .. n] / product [1 .. k-1]
  = -- Daniel's definition
  ((n-k+1)/k) * combina' n (k-1)
  = -- Rearranging again to go back to integer division
  ((n-k+1) * combina' n (k-1)) `quot` k

Putting that all together,

combina' _n 0 = 1
combina' n k = ((n-k+1) * combina' n (k-1)) `quot` k

Just one problem remains: this definition is not tail recursive. We can fix that by counting up instead of down:

combina' n k0 = go 1 1
  where
    go acc k
      | k > k0 = n `seq` acc
      | otherwise = go ((n-k+1)*acc `quot` k) (k + 1)

Don't worry about the n `seq` ; it's of very little consequence.

Note that this implementation uses O(min(k,n-k)) arithmetic operations. So if k and n-k are both very large, it will take a long time. I don't know if there's any efficient way to get exact results in that situation; I believe that binomial coefficients of that sort are usually estimated rather than calculated precisely.