Memoized Factorial function

Posted by WillFa on Fri 03 Sep 2010 06:20 PM — 28 posts, 105,278 views.

WillFa USA Fri 03 Sep 2010 06:20 PM #0

Me and the Guys (that I went to highschool with) have gotten back into Magic the Gathering on tuesday/card nights. So I've been trying to figure out some statistics for any given deck that I've made. I've been scripting out the math for the different combinations of opening hands for any given deck, and thought this function might be useful as a demonstration of the memoizing technique. I think Nick started a thread on it years ago, but hey, hopefully this might start another calculator/geeky fun with code discussion...

(And yes... Mudding, Scripting, Magic... "My Nerdery knows no bounds" </Val Kilmer as Doc Holliday voice>)


function factorial_enclosure()
    local memofac ={[0]=1}
    return function (x)
        if memofac[x] then
            return memofac[x]
        else
            memofac[x] = x * factorial(x-1)
            return memofac[x]
        end
    end
end
factorial = factorial_enclosure()

Why?:
Well, when figuring out combinations of cards, recalculating 60 factorial (60!) (60*59*58...*2*1) over and over again is just a waste of processor cycles. Figure it out once, and store it in a table for later.

For now, let's just assume that every card is unique and there's no duplicates in there. Your opening hand of 7 cards+1 drawn will be 60!/(53!*8!) giving you ~48.3 million possible combinations. Of course, most decks have multiple basic lands and up to 4 of the same card, so the real number of combinations are less, but you see why recalculating that value over and over again can get expensive...

How?:
The factorial_enclosure function is mainly there to have the upvalue table for storing computed results. The recursion of the factorial function also means that factors (smaller factorial values) also get stored. So once you call factorial(60), the calls to factorial (53) and factorial(8) return the memoized value immediately without recursing again.

Twisol USA Fri 03 Sep 2010 07:30 PM #1

WillFa said:
The factorial_enclosure function is mainly there to have the upvalue table for storing computed results.

You can create an artificial scope without resorting to a function:

do
    local memofac ={[0]=1}
    factorial = function (x)
        if memofac[x] then
            return memofac[x]
        else
            memofac[x] = x * factorial(x-1)
            return memofac[x]
        end
    end
end

You can also reduce the number of times you access into the table by using an intermediate local:

do
    local memofac ={[0]=1}
    factorial = function (x)
        local num = memofac[x]
        if not num then
            num = x * factorial(x-1)
            memofac[x] = num
        end
        return memofac[x]
    end
end

You can also work tail recursion into it, except you lose the "memoize every single step aspect of your calculation". :(

--[[ Rules:
  factorial(1, acc) = acc+1
  factorial(x, acc) = factorial(x-1, acc+x)
  factorial(x) = factorial(x-1, x)
--]]

do
  local inner_factorial
  inner_factorial = function(x, result)
    if x == 0 then
      return result
    else
      return inner_factorial(x - 1, result * x)
    end
  end
  
  local memofac = {}
  factorial = function(x)
    if x < 0 then error("X must be greater than 0.") end
    local result = memofac[x]
    if not result then
      result = inner_factorial(x, 1)
      memofac[x] = result
    end
    return result
  end
end

If I could just figure out how to mix tail recursion and memoization here... :S I think this is good enough for most purposes though.

Mudding + scripting = epic win. (I've never played Magic so I can't comment there ;) )

WillFa USA Fri 03 Sep 2010 08:01 PM #2

How about:


do
    local memofac ={[0]=1}
    function factorial (x)
        assert(type(x) == 'number' and x>=0, "X must be a non-negative integer.")
        local result = memofac[x] or x * factorial(x-1)
        memofac[x] = result
        return result
    end
end

I definitely wanted the factors memoized, because I'm to getting into the combinations and hypergeometric distributions... fun fun...

David Haley USA Fri 03 Sep 2010 08:09 PM #3

You can't mix memoization and tail recursion here, because by construction you require the result of the intermediate call in order to store it.

I don't think tail recursion really matters, though, unless you're computing very large factorials without having first computed smaller ones.

If you really needed to work around large stack depth recursion, you could convert it to an iteration and still memoize each step of the way.

Twisol USA Fri 03 Sep 2010 08:10 PM #4

I guess it's a bit of a trade-off then! I like tail recursion, but if you absolutely need memoization, there you go. :)

I prefer to put the assertion outside the main recursive block, because after the first check we only pass valid input to the next level. And you're setting the result to memofac[x] even if you don't need to (i.e. if result = memofac[x] you turn right around and do memofac[x] = result again).

do
    local memofac ={[0]=1}
    local function inner_factorial (x)
        local result = memofac[x]
        if not result then
          result = x * inner_factorial(x-1)
          memofac[x] = result
        end
        return result
    end
    function factorial (x)
      assert(type(x) == 'number' and x>=0, "X must be a non-negative integer.")
      return inner_factorial(x)
    end
end

Twisol USA Fri 03 Sep 2010 08:11 PM #5

David Haley said:
You can't mix memoization and tail recursion here, because by construction you require the result of the intermediate call in order to store it.

Yeah, that's what I thought. =/

Twisol USA Fri 03 Sep 2010 08:27 PM #6

Hey, I cracked it. I just reversed the order of multiplication: instead of going down from the top, I start at 1 and go up. That way the accumulator always contains a valid factorial computation which can be memoized.

do
  local memofac = {[0] = 1}
  
  local inner_factorial
  inner_factorial = function(x, acc, last)
    acc = acc * x
    memofac[x] = acc
    
    if x == last then
      return acc
    else
      return inner_factorial(x+1, acc, last)
    end
  end
  
  factorial = function(x)
    assert(x >= 0, "X must be greater than or equal to 0.")
    if x <= #memofac then
      return memofac[x]
    end
    return inner_factorial(#memofac+1, memofac[#memofac], x)
  end
end

Example:

tprint(memofac)
--[[
  0=1
--]]

Note(factorial(5)) -- 120

tprint(memofac)
--[[
  1=1
  2=2
  3=6
  4=24
  5=120
  0=1
--]]

Note(factorial(6)) -- 720

tprint(memofac)
--[[
  1=1
  2=2
  3=6
  4=24
  5=120
  6=720
  0=1
--]]

WillFa USA Fri 03 Sep 2010 08:28 PM #7

Twisol said:

...
I prefer to put the assertion outside the main recursive block, because after the first check we only pass valid input to the next level. And you're setting the result to memofac[x] even if you don't need to (i.e. if result = memofac[x] you turn right around and do memofac[x] = result again).
...

Good point about the assertion.

I didn't really mind setting the value to the same thing. My thinking was that an assignment to an existing table key is basically one op. 'if not' check being 2 ops. Probably a case of over-optimization before profiling. I'm not too worried about stack depth because a deck is typically 60 cards, maybe 100 for an EDH deck (which I don't play).

Twisol USA Fri 03 Sep 2010 08:37 PM #8

WillFa said:
I'm not too worried about stack depth because a deck is typically 60 cards, maybe 100 for an EDH deck (which I don't play).

Fair enough. It seems to me that tail recursion would be faster, no? It only has to return from effectively one frame, rather than many. Extremely minor, perhaps, but since it's there...

David Haley USA Fri 03 Sep 2010 08:48 PM #9

You know, if you were really worried about speed and wanting to avoid the function calls, you wouldn't be using recursion in the first place here, with or without tail calls.

WillFa USA Fri 03 Sep 2010 08:48 PM #10

Using an iteration and your previous code loop...


do
  local memofac = {[0] = 1}
  
  local function build_memoized(x)
    while #memofac < x do
      local newkey = #memofac+1
      memofac[newkey] = memofac[#memofac] * (newkey)
    end
    return memofac[x]
  end
  
  function factorial(x)
    assert(type(x) == 'number' and x>=0, "X must be a non-negative integer.")
    return memofac[x] or build_memoized(x)
  end
end

I want things memoized for when I start doing this stuff:

http://stattrek.com/Lesson2/Hypergeometric.aspx?Tutorial=Stat (bottom of the page, "Cumulative Hypergeometric Probability")

I came around to David's way of thinking. ;)

Twisol USA Fri 03 Sep 2010 09:15 PM #11

Looks good! :)

(R.I.P. tail-recursive memoized factorial function)

David Haley USA Fri 03 Sep 2010 10:25 PM #12

Incidentally, Twisol, your tail-recursive function is actually slower than the non-tail-recursive function...

$ cat test.lua

do
    local memofac ={[0]=1}
    local function inner_factorial (x)
        local result = memofac[x]
        if not result then
          result = x * inner_factorial(x-1)
          memofac[x] = result
        end
        return result
    end
    function factorial (x)
      assert(type(x) == 'number' and x>=0, "X must be a non-negative integer.")
      return inner_factorial(x)
    end
end

print(os.time())

for i = 1, 10000000 do
    local x = factorial(i)
end

print(os.time())

do
  local memofac = {[0] = 1}
  
  local inner_factorial
  inner_factorial = function(x, acc, last)
    acc = acc * x
    memofac[x] = acc
    
    if x == last then
      return acc
    else
      return inner_factorial(x+1, acc, last)
    end
  end
  
  factorial = function(x)
    assert(x >= 0, "X must be greater than or equal to 0.")
    if x <= #memofac then
      return memofac[x]
    end
    return inner_factorial(#memofac+1, memofac[#memofac], x)
  end
end

for i = 1, 10000000 do
    local x = factorial(i)
end

print(os.time())

$ lua test.lua
1283549016
1283549025
1283549037

so not slower by much, but basically this shows that seeking tail recursion is not always a goal in and of itself. Sometimes the extra overhead you introduce by tracking stuff in other ways outweighs the problem of tail recursion.

Now it turns out that the iterative code WillFa wrote is slower yet:

$ cat test2.lua

do
  local memofac = {[0] = 1}
  
  local function build_memoized(x)
    while #memofac < x do
      local newkey = #memofac+1
      memofac[newkey] = memofac[#memofac] * (newkey)
    end
    return memofac[x]
  end
  
  function factorial(x)
    assert(type(x) == 'number' and x>=0, "X must be a non-negative integer.")
    return memofac[x] or build_memoized(x)
  end
end

for i = 1, 10000000 do
    x = factorial(i)
end

$ time lua test2.lua 
lua test2.lua  17.46s user 0.26s system 98% cpu 17.915 total
$

However, this test happens to be pathological for the iterative code as it has to set up the loop structure each time for just a single increment.

Consider this version, which does things in the other order, and adds a few optimizations:


$ cat test3.lua 

require 'os'

do
  local memofac = {[0] = 1}
  
  function factorial(x)
      assert(type(x) == 'number' and x>=0, "X must be a non-negative integer.")

      if memofac[x] then
          return memofac[x]
      end
      
      for y = #memofac, x-1 do
          memofac[y+1] = memofac[y] * (y+1)
      end
      return memofac[x]
  end
end

print (os.time())

for i = 10000000, 1, -1 do
    x = factorial(i)
end

print (os.time())

$ time lua test3.lua
1283552603
1283552611
lua test3.lua  8.02s user 0.23s system 98% cpu 8.383 total
$

I suspect that were you to run this with JIT, you'd get even better relative performance for the iterative version.

Incidentally, if I take test3.lua and make it compute going forward, we get:

$ time lua test3.lua
1283552677
1283552688
lua test3.lua  11.10s user 0.21s system 98% cpu 11.478 total
$

which is 6 seconds faster than WillFa's.

Twisol USA Fri 03 Sep 2010 10:44 PM #13

David Haley said:
Incidentally, Twisol, your tail-recursive function is actually slower than the non-tail-recursive function...

[...]

so not slower by much, but basically this shows that seeking tail recursion is not always a goal in and of itself. Sometimes the extra overhead you introduce by tracking stuff in other ways outweighs the problem of tail recursion.

They run at the same speed on my computer, actually (10s). Same thing if I run it a second time (10s).

Using GetInfo(232), the high-resolution timer, I get:

53166.273594832
53176.736115121 (10.462s)
53188.700880857 (11.964s)

But if I run it again, I get:

53247.548003798
53257.361734383 (9.813s)
53266.815252998 (9.453s)

There may be some bias towards the non-tail-recursive version, but the difference is probably so minor that there's no reason to prefer one over the other just based on speed.

David Haley USA Fri 03 Sep 2010 11:10 PM #14

Quote:
There may be some bias towards the non-tail-recursive version, but the difference is probably so minor that there's no reason to prefer one over the other just based on speed.

Well... yes. So prefer the simpler one, namely the one without tail recursion. :-) There's less going on, so it's easier to understand. It's a very straightforward implementation of factorial, as you see it in textbooks etc.

Twisol USA Fri 03 Sep 2010 11:29 PM #15

*shrug* It really just counts up instead of down, which is necessary to permit memoization down the line. A non-memoizing tail-recursive function is dead easy: my final code snippet in my first post here does this, except it does memoize the final result.

I'd personally prefer the tail-recursive one, because I know I can pass just about any number and not risk a stack overflow. (Whether Lua can represent the result accurately is a different matter.) Since the speed difference is minute, it's a very self-contained black box, and it's not at risk of overflows, it's the one I'd prefer.

Then again, I'd probably use the iterative one over any of the others.

Nick Gammon Australia Forum Administrator Sat 04 Sep 2010 12:02 AM #16

You guys are only scratching the surface here. See:

http://en.wikipedia.org/wiki/Factorial

:P

David Haley USA Sat 04 Sep 2010 12:09 AM #17

I knew about the optimizations for very large factorials by finding the prime factors, however I don't think the computation here would go that far for it to be useful.

In fact, I think that the memoization described here is not optimal; you are much more likely to be calling the same factorial many times than a whole bunch of them. So it's ok to memoize the "standard" way, namely, "have I seen this question before, if yes, return the answer, if not, compute it from scratch" -- that would considerably improve performance in the case of repeated identical questions, but obviously would perform worse if the common pattern is to do things like I did in my test case.

WillFa USA Sat 04 Sep 2010 12:59 AM #18

Good morning , Nick.

And again, I'm agreeing with David regarding standard memoization versus building a lookup table.

I've started building the combination and hypergeometric distribution functions, and I'm seeing trends on common values. To list a few major ones: 60 card deck. 1,2,3, or 4 instances of most cards. 17-25 lands... It's easier to just memoize those values.

As an aside, since David and Nick appear to know more about this Math type stuff... ;)

I'm stumped about combining the probabilities though. i.e. I want to track the probability that I'll have an opening hand of 2 lands and a "1 drop" and a "2 drop" (1 mana and 2 mana casting cost) creatures in my first 9 cards (two turns).

I have the algorithm for finding the probability of having 2 lands in the first 2 turns. h(21,60,9,2) which is (c(21, 2) * c(60-21, 9-2))/ c(60,9) . c computes combinations.
Likewise either of the critters can be found. h(9,60,9,1)

But how would I combine those 3 probabilities to find the intersection of them?

Nick Gammon Australia Forum Administrator Sat 04 Sep 2010 04:46 AM #19

I am no great expert in probability, but I think you multiply them. For example, if there probability of being male is 0.5 and the probability of having red hair is 0.3 the the probability of a random person being a male red-head is 0.5 * 0.3, namely 0.15.

David Haley USA Sat 04 Sep 2010 06:49 AM #20

Unfortunately it is not so easy. Multiplication only works when the events are independent. If you assume that the events "child is male" and "child is red-haired" have any correlation, then a straight multiplication will misstate the probabilities.

There are a few ways to approach this; one is called the inclusion/exclusion principle:
http://en.wikipedia.org/wiki/Inclusion-exclusion_principle

and in particular:
http://en.wikipedia.org/wiki/Inclusion-exclusion_principle#In_probability

Another is to compute conditional probabilities; they are defined as:

P(A|B) = P(A and B) / P(B)

Or in other words, the probability of event A occurring given that B has occurred is equal to the probability that A and B occur at the same time divided by the probability that B occurs on its own.

You can rearrange this to say that the probability of A and B is equal to P(A|B) * P(B). Then you need only compute P(A|B), which in and of itself might be non-obvious.

By the way, let's see how this applies in the red-haired male child example, assuming independence of masculinity and red-hair.

P(Red and Male) = P(Red given Male) / P(Male)
= 0.3 / 0.5
= 0.15

which is exactly the result of multiplying the probabilities. Note that because we assumed independence, we were able to write that P(Red given Male) = P(Red). Had we not assumed independence, we would not have been able to write that.

Of course in the card-drawing case independence won't apply. Clearly, if you drew 7 lands, you wouldn't be able to draw 5 other cards in your first hand. So the fact that one of them held true will affect the other.

I'd have to sit down and think about this for a while longer and get all the parameters before I could give an actual answer to this, and besides combinatorics were never my strong suit, but it'd be fun to give it a shot. You might even be able to derive an answer empirically if you are willing to make enough simplifying assumptions. (The more simplifications you make, the fewer combinations there are to look at.)

Nick Gammon Australia Forum Administrator Sat 04 Sep 2010 07:16 AM #21

WillFa said:

For now, let's just assume that every card is unique and there's no duplicates in there. ... Of course, most decks have multiple basic lands and up to 4 of the same card ...

If you don't mind my pointing it out, these two statements completely contradict each other.

I'm not totally sure how probability is going to help you here. If there are quite a few lands in the pack, then drawing one doesn't really lower the probability of getting another (by much). However if there is a unique card, then yes, drawing it lowers the probability of another to zero.

Besides the rules can change. To quote from them: "Whenever a card's text directly contradicts the rules, the card takes precedence."

So woe betide the player who draws a card that states that the laws of probability have become invalid.

WillFa USA Sat 04 Sep 2010 02:12 PM #22

would

Quote:
For illustration, let's just assume that every card is unique and there's no duplicates in there. ... In reality, most decks have multiple basic lands and up to 4 of the same card ...

Have been better?

David Haley USA Sat 04 Sep 2010 05:36 PM #23

But even then, you're looking for the probability of drawing two lands -- so even if the lands are different, there is at least some notion of categorization going on. So I think you can't assume total uniqueness.

Quote:
So woe betide the player who draws a card that states that the laws of probability have become invalid.

Made me smile. :-)

Nick Gammon Australia Forum Administrator Sat 04 Sep 2010 08:49 PM #24

WillFa said:

would [ .. a rewording of the problem ..] Have been better?

I'm just trying to gently suggest that if you look up probability ideas for (say) a deck of 52 playing cards (which are indeed unique) and try to extrapolate to the Magic cards, where you acknowledge duplicates, the results may not hold.

WillFa USA Sat 04 Sep 2010 11:34 PM #25

But looking up probabilities for decks of cards leads to urn problems which give examples of multivariate hypergeometric distributions... which is what I needed. :)

(I had the hypergeometric functions... just needed to find the multivariate ones to make headway. :) )

WillFa USA Mon 06 Sep 2010 08:14 AM #26

If anyone's curious... (which I don't think anyone is...)



do
  local memofac = {[0] = 1}

  local function memo_fac(x)
    while #memofac < x do
      local newkey = #memofac+1
      memofac[newkey] = memofac[#memofac] * (newkey)
    end
    return memofac[x]
  end

  function factorial(x)
    assert(type(x) == 'number' and x>=0, "X must be a non-negative integer.", 2)
    return memofac[x] or memo_fac(x)
  end
end


function comb (n,k)
    assert (type(n) == "number" and type(k) == "number" and n >= k, 
      "Expected numbers for Population and Samples. Samples must be smaller than Population")
    return factorial(n)/(factorial(n-k)*factorial(k))
end

function draw (deckSize, handsize, ...)
    local numInDeck = 1
    local numDrawn = 2
    local desired = {...}
    local padding = {deckSize, 0}
    local padDraws = 0
    local prob = 0

    for _, v in ipairs(desired) do
        padding[numInDeck] = padding[numInDeck] - v[numInDeck]      --number of cards not in the requested sets
        padding[numDrawn] = padding[numDrawn] + v[numDrawn]         --total number of cards requested
    end

    for _,i in ipairs(desired) do
        local workingProb = 0
        for j = 0, math.min( i[numInDeck]-i[numDrawn], handsize - padding[numDrawn] ) do
            workingProb = comb(padding[numInDeck] , handsize - padding[numDrawn] - j)
            for _,k in ipairs(desired) do
                if k == i then
                    workingProb = workingProb * comb(k[numInDeck], k[numDrawn] + j)
                else
                    workingProb = workingProb * comb(k[numInDeck], k[numDrawn])
                end
            end
            workingProb = workingProb / comb(deckSize, handsize)
            prob = prob + workingProb
        end
    end

    return prob
end

This gives a "draw" function that will take in the size of the deck, how many cards you're drawing, and then any number of tables with the number of cards in the deck and how many you want...


print(draw(60, 7, {17,2}, {8,2}, {4,1}, {8,1}))

0.038349597671632

The formula btw is:


 comb(17,2) * comb(8,2) * comb(4,1) * comb(8,1) * comb(60 -17-8-4-8 , 5 -2-2-1-1) / comb(60, 5)

That figures out the probability of getting EXACTLY 2 of a, 2 of b, 1 of c, and 1 of d, and 1 not in the previous sets. The loop sums the probabilities so you're getting AT LEAST 2,2,1,1 of the requested cards... (i.e. 3,2,1,1,0 ; 2,3,1,1,0 ; 2,2,2,1,0; 2,2,1,2,0...)

Twisol USA Mon 06 Sep 2010 08:31 AM #27

Very nice! :)

Incidentally...

local factorial = setmetatable({[0] = 1}, {
  __index = function(t, x)
    assert(type(x) == 'number' and x>=0, "X must be a non-negative integer.", 2)
    while #t < x do
      local newkey = #t+1
      t[newkey] = rawget(t, #t) * newkey
    end
    return rawget(t, x)
  end,
  
  __call = function(t, x)
    return t[x]
  end,
})

Note(factorial(50))

This form greatly appeals to me. ^_^