Register forum user name Search FAQ

Gammon Forum

Notice: Any messages purporting to come from this site telling you that your password has expired, or that you need to verify your details, confirm your email, resolve issues, making threats, or asking for money, are spam. We do not email users with any such messages. If you have lost your password you can obtain a new one by using the password reset link.

Due to spam on this forum, all posts now need moderator approval.

Entire forum ➜ MUSHclient ➜ Lua ➜ Memoized Factorial function

Memoized Factorial function

It is now over 60 days since the last post. This thread is closed. Refresh page

Pages: 1 2

Posted by WillFa USA (525 posts) Bio

Date

Fri 03 Sep 2010 06:20 PM (UTC)

Amended on Fri 03 Sep 2010 06:26 PM (UTC) by WillFa

Message

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.

Posted by Twisol USA (2,257 posts) Bio

Date

Reply #1 on Fri 03 Sep 2010 07:30 PM (UTC)

Amended on Fri 03 Sep 2010 08:13 PM (UTC) by Twisol

Message

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 ;) )

'Soludra' on Achaea

Blog: http://jonathan.com/
GitHub: http://github.com/Twisol

Posted by WillFa USA (525 posts) Bio

Date Reply #2 on Fri 03 Sep 2010 08:01 PM (UTC)

Message

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...

Posted by

David Haley USA (3,881 posts) Bio

Date

Reply #3 on Fri 03 Sep 2010 08:09 PM (UTC)

Message

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.

David Haley aka Ksilyan
Head Programmer,
Legends of the Darkstone

http://david.the-haleys.org

Posted by Twisol USA (2,257 posts) Bio

Date Reply #4 on Fri 03 Sep 2010 08:10 PM (UTC)

Message

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

'Soludra' on Achaea

Blog: http://jonathan.com/
GitHub: http://github.com/Twisol

Posted by

Twisol USA (2,257 posts) Bio

Date

Reply #5 on Fri 03 Sep 2010 08:11 PM (UTC)

Message

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. =/

'Soludra' on Achaea

Blog: http://jonathan.com/
GitHub: http://github.com/Twisol

Posted by Twisol USA (2,257 posts) Bio

Date

Reply #6 on Fri 03 Sep 2010 08:27 PM (UTC)

Amended on Fri 03 Sep 2010 08:56 PM (UTC) by Twisol

Message

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

'Soludra' on Achaea

Blog: http://jonathan.com/
GitHub: http://github.com/Twisol

Posted by

WillFa USA (525 posts) Bio

Date

Reply #7 on Fri 03 Sep 2010 08:28 PM (UTC)

Message

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).

Posted by

Twisol USA (2,257 posts) Bio

Date

Reply #8 on Fri 03 Sep 2010 08:37 PM (UTC)

Message

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...

'Soludra' on Achaea

Blog: http://jonathan.com/
GitHub: http://github.com/Twisol

Posted by

David Haley USA (3,881 posts) Bio

Date

Reply #9 on Fri 03 Sep 2010 08:48 PM (UTC)

Message

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.

David Haley aka Ksilyan
Head Programmer,
Legends of the Darkstone

http://david.the-haleys.org

Posted by WillFa USA (525 posts) Bio

Date

Reply #10 on Fri 03 Sep 2010 08:48 PM (UTC)

Amended on Fri 03 Sep 2010 08:49 PM (UTC) by WillFa

Message

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. ;)

Posted by

Twisol USA (2,257 posts) Bio

Date

Reply #11 on Fri 03 Sep 2010 09:15 PM (UTC)

Message

Looks good! :)

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

'Soludra' on Achaea

Blog: http://jonathan.com/
GitHub: http://github.com/Twisol

Posted by David Haley USA (3,881 posts) Bio

Date

Reply #12 on Fri 03 Sep 2010 10:25 PM (UTC)

Amended on Fri 03 Sep 2010 10:26 PM (UTC) by David Haley

Message

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.

David Haley aka Ksilyan
Head Programmer,
Legends of the Darkstone

http://david.the-haleys.org

Posted by

Twisol USA (2,257 posts) Bio

Date

Reply #13 on Fri 03 Sep 2010 10:44 PM (UTC)

Amended on Fri 03 Sep 2010 10:45 PM (UTC) by Twisol

Message

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.

'Soludra' on Achaea

Blog: http://jonathan.com/
GitHub: http://github.com/Twisol

Posted by

David Haley USA (3,881 posts) Bio

Date

Reply #14 on Fri 03 Sep 2010 11:10 PM (UTC)

Message

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.

David Haley aka Ksilyan
Head Programmer,
Legends of the Darkstone

http://david.the-haleys.org

The dates and times for posts above are shown in Universal Co-ordinated Time (UTC).

To show them in your local time you can join the forum, and then set the 'time correction' field in your profile to the number of hours difference between your location and UTC time.

82,870 views.

This is page 1, subject is 2 pages long: 1 2

It is now over 60 days since the last post. This thread is closed. Refresh page

Go to topic: Search the forum

Information and images on this site are licensed under the Creative Commons Attribution 3.0 Australia License unless stated otherwise.