Saturday, December 31, 2011

Picomath

The Picomath project holds some reusable math functions inspired by John D. Cook's Stand-alone code for numerical computing, including:

  • Error function
  • Phi (standard normal CDF)
  • Phi inverse
  • Gamma
  • Log Gamma
  • exp(x) - 1 (for small x)
  • log(n!)

These functions are implemented in an impressive list of languages: Ada, C++, C#, D, Erlang, Go, Haskell, Java, Javascript, Lua, Pascal, Perl, PHP, Python (2.x and 3.x), Ruby, Scheme, and Tcl.

And now Factor!

You can find the code (and a bunch of tests) for this on my Github.

Friday, December 30, 2011

Slot Machines

Playing slot machines can be pretty fun, but don't be fooled by claims that the casino has the "loosest slots", odds are probably still against you. I thought it would be fun (and cheaper!) to build a slot machine simulator using Factor.

Our slot machine is going to be a console application that will look something like this:

Spinning

Even though our slot machine is text-only, we can still make use of some nice unicode characters to be our symbols:

CONSTANT: SYMBOLS "☀☁☂☃"

Each spin chooses a different symbol at random (each being equally likely):

: spin ( value -- value' )
    SYMBOLS remove random ;

To reproduce the feel of spinning a slot machine, we will introduce a slight delay so that the wheel spins fast at the beginning and then slower and slower until it stops on a symbol:

: spin-delay ( n -- )
    15 * 25 + milliseconds sleep ;

Spinning the slot machine takes a spin number, delays for a bit, then rotates each wheel (we stop spinning the first column after 10 spins, the second after 15, and the last after 20):

: spin-slots ( a b c n -- a b c )
    {
        [ spin-delay ]
        [ 10 < [ [ spin ] 2dip ] when ]
        [ 15 < [ [ spin ]  dip ] when ]
        [ drop spin ]
    } cleave ;

Display

Each "spin" of the slot machine will be printed out. Using ANSI escape sequences, we move the cursor to the top left ("0,0") of the screen and then issue a clear screen instruction. Then we print out the current display and flush the output to screen:

: print-spin ( a b c -- a b c )
    "\e[0;0H\e[2J" write
    "Welcome to the Factor slot machine!" print nl
    "  +--------+" print
    "  | CASINO |" print
    "  |--------| *" print
    3dup "  |%c |%c |%c | |\n" printf
    "  |--------|/" print
    "  |    [_] |" print
    "  +--------+" print flush ;

Playing

The player will have won if, after all the spins, the "pay line" shows three of the same characters:

: winner? ( a b c -- )
    3array all-equal? nl "You WIN!" "You LOSE!" ? print nl ;

Playing the slot machine consists of spinning the wheels 20 times, displaying each spin to the user, then checking if the user has won the game.

: play-slots ( -- )
    f f f 20 iota [ spin-slots print-spin ] each winner? ;

Since our casino wants the user to keep playing, we make it really easy to just hit ENTER to continue:

: continue? ( -- ? )
    "Press ENTER to play again." write flush readln ;

And, to finish it off, we define a "MAIN" entry point that will be run when the script is executed:

: main-slots ( -- )
    [ play-slots continue? ] loop ;

MAIN: main-slots

The code for this is on my Github.

Wednesday, December 28, 2011

Elementology

Question: What do the words bamboo, crunchy, finance, genius, and tenacious have in common? I'll give you a hint: its the same thing they have in common with the words who, what, when, where, and how?

Stumped? Well, it's not that these are all English words.

Answer: All of these words can be spelled using elements from the periodic table!

I was recently inspired by the Periodic GeNiUS T-shirt from ThinkGeek and a website that can "make any words out of elements in the periodic table". I thought it would be fun to use Factor to see how many other words can be spelled using the symbols for chemical elements.

First, we need a list of elements:

: elements ( -- assoc )
    H{
        { "H" "Hydrogen" }
        { "He" "Helium" }
        { "Li" "Lithium" }
        { "Be" "Beryllium" }
        { "B" "Boron" }
        { "C" "Carbon" }
        ...
        { "Uut" "Ununtrium" }
        { "Uuq" "Ununquadium" }
        { "Uup" "Ununpentium" }
        { "Uuh" "Ununhexium" }
        { "Uus" "Ununseptium" }
        { "Uuo" "Ununoctium" }
    } [ [ >lower ] dip ] assoc-map ;

Next, a word that checks if a particular substring is the symbol of an element:

: element? ( from to word -- ? )
    2dup length > [ 3drop f ] [ subseq elements key? ] if ;

We know that symbols are only ever one, two, or three characters. A word is considered "periodic" if it can be composed of any number of (possibly repeating) element symbols. We build a recursive solution that starts with the first character and continues as long as element symbols are a match or until the end of the word is reached:

: (periodic?) ( word from -- ? )
    {
        [ swap length = ]
        [
            { 1 2 3 } [
                dupd + [ pick element? ] keep
                '[ dup _ (periodic?) ] [ f ] if
            ] with any? nip
        ]
    } 2|| ;

: periodic? ( word -- ? )
    >lower 0 (periodic?) ;

It's easy to get a list of dictionary words from most Unix systems:

: dict-words ( -- words )
    "/usr/share/dict/words" ascii file-lines ;

And then a list of all "periodic words":

: periodic-words ( -- words )
    dict-words [ periodic? ] filter ;

So, how many words are "periodic words"? About 13.7% of them.

IN: scratchpad dict-words length .
235886

IN: scratchpad periodic-words length .
32407

The code for this is on my Github.

Thursday, November 10, 2011

wc -l

The wc program is a utility that can count the number of lines in a file. It has a number of options that are described on its man page that can change its function to count characters or word or produce the length of the longest line.

For the last month, Joe Groff has been improving the performance of Factor's I/O libraries. Yesterday, we were investigating slow performance when doing lots of small reads from a file. Joe was able to make a number of nice speedups, which will be in the next release of Factor.

After suggesting that we use wc -l as a benchmark to aspire to, we came up with several approaches with varying performance. I want to demonstrate these, using timing information from running it on my computer. Although the Factor image file is binary data, we are going to count the number of newline characters in it.

Note: some of these require the latest development version of Factor to run.

Our "gold standard" will be wc -l, which takes just over 0.1 seconds:

$ time wc -l factor.image
 6149212 factor.image

real 0m0.111s
user 0m0.090s
sys  0m0.020s

Our first attempt in Factor is the shortest amount of code but takes 5.8 seconds (you can time this by running "USE: system [ image wc-file-lines ] time"):

: wc-file-lines ( path -- n )
    binary file-lines length ;

We don't really need the lines, just their count, so perhaps just increment each-line in a loop. This is an improvement at just over 3 seconds:

: wc-each-line ( path -- n )
    binary [ 0 [ drop 1 + ] each-line ] with-file-reader ;

Trying to use read-until to look for the next newline, is a bit slower at 3.4 seconds:

: wc-read-until ( path -- n )
    binary [
        0 [ "\n" read-until [ drop 1 + ] dip ] loop
    ] with-file-reader ;

Instead of reading each line at a time, we can just read 65,536 byte blocks and count the number of newlines. This takes about 1.5 seconds:

: wc-each-block ( path -- n )
    binary [
        0 [ [ CHAR: \n = ] count + ] each-block
    ] with-file-reader ;

Since we are only counting characters in each block, we don't need to allocate and copy the bytes out of the I/O buffer. Instead, we can look at a slice. This takes about 1 second:

: wc-each-block-slice ( path -- n )
    binary [
        0 [ [ CHAR: \n = ] count + ] each-block-slice
    ] with-file-reader ;

The stream functions (such as read) operate on an input-stream dynamic variable, which introduces some overhead. If we remove that, the compiler can eliminate some of the dynamic dispatches. taking 0.320 seconds:

: wc-fast ( path -- n )
    binary <file-reader> [
        0 swap [
            [ CHAR: \n = ] count +
        ] each-stream-block-slice
    ] with-disposal ;

If we make an assumption that the number of lines in a file will fit into a fixnum, then we can get a bit faster at 0.240 seconds:

: wc-faster ( path -- n )
    binary <file-reader> [
        0 swap [
            [ CHAR: \n = ] count + >fixnum
        ] each-stream-block-slice
    ] with-disposal ;

And, if we cheat and just run the wc -l process directly, we can get to 0.210 seconds:

: wc-system ( path -- n )
    "wc -l " prepend utf8 [
        readln " " split harvest first string>number
    ] with-process-reader ;

Overall, not too bad! Factor is getting within shouting distance of programs written in "faster" languages. Probably with a few more hours, we could close the gap, but thats enough for today.

Update: using SIMD, Joe was able to get the time down to 0.120 seconds!
: aligned-slices ( seq -- head tail )
   dup length 15 bitnot bitand cut-slice ; inline

: wc-simd ( path -- n )
   [
       0 swap binary <file-reader> &dispose [
           aligned-slices [
               uchar-16 cast-array swap
               [ 10 uchar-16-with v= vcount + >fixnum ] reduce
           ] [ [ 10 = ] count + >fixnum ] bi*
       ] each-stream-block-slice
   ] with-destructors ;

Wednesday, October 26, 2011

Disassemble

A neat trick that Factor provides is the ability to disassemble functions into the machine code that is generated by the compiler. In 2008, Slava Pestov created a disassembler, and has improved it a bit since then (switching to udis86 for its implementation).

Constant Folding

The compiler performs constant folding, using the compiler.tree.debugger vocabulary, you can output the optimized form of a quotation:

( scratchpad ) [ 2 2 + ] optimized.
[ 4 ]

Using the disassembler, you can see the machine code this generates:

( scratchpad ) [ 2 2 + ] disassemble
011c1a5530: 4983c608        add r14, 0x8
011c1a5534: 49c70640000000  mov qword [r14], 0x40
011c1a553b: c3              ret 
011c1a553c: 0000            add [rax], al
011c1a553e: 0000            add [rax], al

Local Variables

One of the questions that comes up sometimes is whether local variables affect performance. We can examine two words that add numbers together, one using locals and one just using the stack:

( scratchpad ) : foo ( x y -- z ) + ;

( scratchpad ) :: bar ( x y -- z ) x y + ;

The "optimized output" looks a little different:

( scratchpad ) \ foo optimized.
[ + ]

( scratchpad ) \ bar optimized.
[ "COMPLEX SHUFFLE" "COMPLEX SHUFFLE" R> + ]

But, the machine code that is generated is identical:

( scratchpad ) \ foo disassemble
01115de7b0: 488d1d05000000  lea rbx, [rip+0x5]
01115de7b7: e9e49439ff      jmp 0x110977ca0 (+)
01115de7bc: 0000            add [rax], al
01115de7be: 0000            add [rax], al

( scratchpad ) \ bar disassemble
01115ef620: 488d1d05000000  lea rbx, [rip+0x5]
01115ef627: e9748638ff      jmp 0x110977ca0 (+)
01115ef62c: 0000            add [rax], al
01115ef62e: 0000            add [rax], al

Dynamic Variables

Another frequently used feature is dynamic variables, implemented by the namespaces vocabulary. For example, the definition of the print word looks for the current value of the output-stream variable and then calls stream-print on it:

( scratchpad ) \ print see
USING: namespaces ;
IN: io
: print ( str -- ) output-stream get stream-print ; inline

The optimized output inlines the implementation of get:

( scratchpad ) [ "Hello, world" print ] optimized.
[
    "Hello, world" \ output-stream 0 context-object assoc-stack
    stream-print
]

You can inspect the machine code generated, seeing references to the factor words that are being called:

( scratchpad ) [ "Hello, world" print ] disassemble
011c0c6c40: 4c8d1df9ffffff        lea r11, [rip-0x7]
011c0c6c47: 6820000000            push dword 0x20
011c0c6c4c: 4153                  push r11
011c0c6c4e: 4883ec08              sub rsp, 0x8
011c0c6c52: 4983c618              add r14, 0x18
011c0c6c56: 48b8dbc5a31a01000000  mov rax, 0x11aa3c5db
011c0c6c60: 498946f0              mov [r14-0x10], rax
011c0c6c64: 498b4500              mov rax, [r13+0x0]
011c0c6c68: 488b4040              mov rax, [rax+0x40]
011c0c6c6c: 498906                mov [r14], rax
011c0c6c6f: 48b86c91810e01000000  mov rax, 0x10e81916c
011c0c6c79: 498946f8              mov [r14-0x8], rax
011c0c6c7d: e8de4e36ff            call 0x11b42bb60 (assoc-stack)
011c0c6c82: 4883c418              add rsp, 0x18
011c0c6c86: 488d1d05000000        lea rbx, [rip+0x5]
011c0c6c8d: e94e5264ff            jmp 0x11b70bee0 (stream-print)
011c0c6c92: 0000                  add [rax], al
011c0c6c94: 0000                  add [rax], al
011c0c6c96: 0000                  add [rax], al
011c0c6c98: 0000                  add [rax], al
011c0c6c9a: 0000                  add [rax], al
011c0c6c9c: 0000                  add [rax], al
011c0c6c9e: 0000                  add [rax], al

Wednesday, October 12, 2011

Optimizing 2^x

A great little article was posted last year about optimizing 2^x for doubles (by approximation). The author gets 40+% performance improvements in C#. I wondered if we could get similar improvements in Factor.

The basic idea is to use the fact that a number, ab+c can be factored into ab * ac. If we separate a floating point number into two components: an integer and a fractional part, we can show that:

2n = 2integer(n) * 2fractional(n).

We are going to approximate this value by using a lookup table to compute the fractional part (within a specified precision). For example, to compute within a 0.001 precision, we need 1000 lookup values, essentially performing this calculation:

2n = ( 1 << Int(n) ) * Table[ (int) ( Frac(n) * 1000 ) ];

Implementation

So, we need a word that can split a floating point number into those two values:

: float>parts ( x -- float int )
    dup >integer [ - ] keep ; inline

Instead of one table with 1000 values, we will copy the original authors decision to use three lookup tables for additional precision. The following code calculates these lookup tables:

CONSTANT: BITS1 10
CONSTANT: BITS2 $[ BITS1 2 * ]
CONSTANT: BITS3 $[ BITS1 3 * ]

CONSTANT: PRECISION1 $[ 1 BITS1 shift ]
CONSTANT: PRECISION2 $[ 1 BITS2 shift ]
CONSTANT: PRECISION3 $[ 1 BITS3 shift ]

CONSTANT: MASK $[ PRECISION1 1 - ]

CONSTANT: FRAC1 $[ 2 PRECISION1 iota [ PRECISION1 / ^ ] with map ]
CONSTANT: FRAC2 $[ 2 PRECISION1 iota [ PRECISION2 / ^ ] with map ]
CONSTANT: FRAC3 $[ 2 PRECISION1 iota [ PRECISION3 / ^ ] with map ]

The function pow2 looks pretty similar to our original mathematical definition:

: pow2 ( n -- 2^n )
    >float 2^int 2^frac * >float ;

The guts of the implementation is in the 2^int and 2^frac words:

: 2^int ( n -- 2^int frac )
    [ float>parts ] keep 0 >= [ 1 swap shift ] [
        over 0 < [ [ 1 + ] [ 1 - ] bi* ] when
        1 swap neg shift 1.0 swap /
    ] if swap ; inline

: 2^frac ( frac -- 2^frac )
    PRECISION3 * >fixnum
    [ BITS2 neg shift FRAC1 nth-unsafe ]
    [ BITS1 neg shift MASK bitand FRAC2 nth-unsafe ]
    [ MASK bitand FRAC3 nth-unsafe ] tri * * ; inline

Testing

Let's try it and see how well it works for small values:

( scratchpad ) 2 1.5 ^ .
2.82842712474619

( scratchpad ) 1.5 pow2 .
2.82842712474619

It seem's to work, how about larger values:

( scratchpad ) 2 16.3 ^ .
80684.28027297248

( scratchpad ) 16.3 pow2 .
80684.28026255539

The error is clearly detectable, but to test that it really works the way we expect it too, we will need to calculate relative error:

: relative-error ( approx value -- relative-error )
    [ - abs ] keep / ;

Our test case will generate random values, compute 2x using our approximation and verify the relative error is less than 0.000000001 when compared with the correct result:

[ t ] [
    10000 [ -20 20 uniform-random-float ] replicate
    [ [ pow2 ] [ 2 swap ^ ] bi relative-error ] map
    supremum 1e-9 <
] unit-test

And to verify performance, we will benchmark it against the built-in (and more accurate ^ word):

: pow2-test ( seq -- new old )
    [ [ pow2 drop ] [ each ] benchmark ]
    [ 2 swap [ ^ drop ] with [ each ] benchmark ] bi ;

We got almost a 40% performance improvement at the cost of a small loss of precision - not bad!

( scratchpad ) 10000 [ -20 20 uniform-random-float ] replicate
               pow2-test / >float .
0.6293153754782392

The code for this is on my Github.

Sunday, September 25, 2011

Approximating Pi

A few days ago, it was announced on the Wolfram Blog that a 13-year-old had made a record calculating 458 million terms for the continued fraction of pi. In the spirit of that, I thought I would show how to solve a question that sometimes gets asked at interviews:

Given that Pi can be estimated using the function 4 * (1 - 1/3 + 1/5 - 1/7 + ...) with more terms giving greater accuracy, write a function that calculates Pi to an accuracy of 5 decimal places.

Using Factor, we can calculate the nth approximation of pi using vector arithmetic:

: approximate-pi ( n -- approx )
    [1,b] 2 v*n 1 v-n 1 swap n/v
    [ odd? [ neg ] when ] map-index sum 4 * ;

This isn't ideal if we want to try an increasing number of terms (looking for a particularly accuracy), since a lot of the the work would be redone unnecessarily. Instead, we can write a word that adds successive terms until the difference between the previous approximation and the current approximation is less than our requested accuracy.

: next-term ( approx i -- approx' )
    [ 2 * 1 + ] [ odd? [ neg ] when ] bi 4.0 swap / + ; inline

:: find-pi-to ( accuracy -- n approx )
    1 4.0 [
        dup pick next-term [ - ] keep
        swap abs accuracy >= [ 1 + ] 2dip
    ] loop ;

To show its performance, we can time it:

( scratchpad ) [ 0.00001 find-pi-to ] time .
Running time: 0.026030341 seconds

3.141597653564762

An equivalent function in Python might look like this:

def find_pi_to(accuracy):
    i = 1
    approx = 4.0
    while 1:
        term = (2 * i) + 1
        if i % 2 == 1:
            term = -term
        new = approx + 4.0/term
        if abs(new - approx) < accuracy:
            approx = new
            break
        i += 1
        approx = new
    return approx

But, if we time this version (not counting startup or compile time), it takes 0.134 seconds. Doing the math shows that Factor is 5 times faster than Python in this case. Not bad.