Thursday, April 19, 2012

Twin Primes

The most recent programming challenge from Programming Praxis is to:

Pairs of prime numbers that differ by two are known as twin primes: (3,5), (5,7), (11,13), (17,19), (29,31), (41,43), (59,61), (71,73),... Your task is to write a function that finds the twin primes less than a given input n.

Samuel Tardieu has contributed many improvements to Factor's math.primes vocabulary, which we will be using to solve this puzzle.

We can solve this puzzle in the naive method by computing all prime numbers up to a specified input, and then filtering them for pairs that differ by two:

: twin-primes-upto ( n -- seq )
    primes-upto 2 <clumps> [ first2 - abs 2 = ] filter ;

Another nice word would check to see whether any two numbers are twin primes (using short-circuit combinators to exit early if any of the conditions are not satisfied):

: twin-primes? ( x y -- ? )
    { [ - abs 2 = ] [ nip prime? ] [ drop prime? ] } 2&& ;

The puzzle page suggests a more efficient method for computing twin primes, which might be worth experimenting with...

Thursday, April 5, 2012

Faster Big Ratios!

While working on a post about approximating pi, I noticed that the performance of Factor's large rational numbers was less than stellar.

Specifically, we defined this estimation function to find pi as a rational number:

:: find-pi-to ( accuracy -- n approx )
    1 4 [
        over [ 2 * 1 + ] [ odd? [ neg ] when ] bi
        4 swap / [ + ] keep
        abs accuracy >= [ 1 + ] 2dip
    ] loop ;

Using this for an accuracy of 0.001 was incredibly slow:

IN: scratchpad [ 0.001 find-pi-to ] time
Running time: 1.692090043 seconds

Bignum

In Factor, large integers are stored as arbitrary-precision "bignums". Similar to other languages such as Python and Scheme, Factor stores these numbers as a sequence of large (either 30-bit or 62-bit) digits, and then performs math on this sequence of digits.

It turns out that most of the ratios has both a bignum numerator and bignum denominator. For most basic math operation on ratios, Factor would compute the greatest common divisor to produce a normalized fraction (e.g., 6/4 would become 3/2).

For an accuracy of 0.001 in this example, the denominator had more than 1,700 digits in base 10. As these numbers grew, the performance of Euclid's GCD algorithm began degrading particularly due to the (relatively) expensive calculation of "mod" for bignum's.

Lehmer GCD

I noticed a Python bug report that suggested addressing a similar performance problem for their rational number implementation (in fractions.py) by implementing Lehmer's GCD algorithm.

After implementing a bignum-gcd primitive that used Lehmer's GCD, I created a fast-gcd word that used this for bignum's and the current gcd word for other real numbers.

Performance improvement was impressive! After recompiling with these patches, our original test case takes less than 10% of the time!

IN: scratchpad [ 0.001 find-pi-to ] time
Running time: 0.156713844 seconds

You can find this in the latest development version of Factor.

Friday, March 30, 2012

Mega Millions

With the Mega Millions jackpot over $640 million, you might be asking yourself:

How do I use Factor to win the lottery???

You can find a good article on Forbes about calculating your odds of winning. We can also calculate it with Factor. The lottery works like this: you have to choose 5 "main" numbers (1 through 56) and then pick a "mega" number (1 through 46).

We can use the "n choose k" formula in math.combinatorics to compute the number of ways of picking "5 unordered outcomes from 56 numbers and then 1 of 46 possible mega numbers":

IN: scratchpad 56 5 nCk 46 * .
175711536

Sure enough, thats 175,711,536 possibilities. So if you buy one random lottery ticket, you have a 1 in 175,711,536 chance of winning, or a 0.00000000569114597006% chance. Not that great! What are the other odds?

Since the jackpot is so large, it's got to be worth playing, right? In fact, if you take the total jackpot and the odds of each winning category, you can find the expected value of a ticket:

IN: scratchpad {
                   640000000/175711536
                   250000/3904701
                   10000/689065
                   150/15313
                   150/13781
                   7/306
                   10/844
                   3/141
                   2/75
               } sum "$%.2f" printf
$3.82

Wow, $3.82 of expected value (not counting sharing the jackpot with someone who picked the same numbers, or the fact you can only win the highest prize you qualify for)! But, how do we pick our numbers... well, Factor to the rescue! Let's randomly sample our 5 main numbers and then pick a random mega number and output the result:

IN: scratchpad : mega ( -- seq )
                   56 [1,b] 5 sample natural-sort
                   46 [1,b] random suffix ;

In the spirit of fiver, we can generate these numbers to play:

IN: scratchpad 5 [ mega ] replicate .
{
    { 2 8 25 30 43 29 }
    { 6 15 31 41 47 33 }
    { 11 20 26 29 33 15 }
    { 3 6 19 23 32 15 }
    { 16 25 31 44 50 24 }
}

And tonight we can check the winning numbers and see if it worked!

Friday, March 2, 2012

Next Permutation

Yesterday, I noticed a programming challenge to find the "next greater permutation of digits" for a given number:

Given a number, find the next higher number that uses the same set of digits. For instance, given the number 38276, the next higher number that uses the digits 2 3 6 7 8 is 38627.

While reading the comments, I noticed that some of the C++ solutions used the std::next_permutation function that returns the "lexicographically next greater permutation of elements". Noticing that the math.combinatorics vocabulary lacks a next-permutation word, I thought it would be nice to contribute one to the Factor standard library.

std::next_permutation()

It's a useful place to start to get an overview and example of how the algorithm works. The C++ version is pretty dense and looks like this:

template<typename Iter>
bool next_permutation(Iter first, Iter last)
{
    if (first == last)
        return false;
    Iter i = first;
    ++i;
    if (i == last)
        return false;
    i = last;
    --i;

    for(;;)
    {
        Iter ii = i;
        --i;
        if (*i < *ii)
        {
            Iter j = last;
            while (!(*i < *--j))
            {}
            std::iter_swap(i, j);
            std::reverse(ii, last);
            return true;
        }
        if (i == first)
        {
            std::reverse(first, last);
            return false;
        }
    }
}

Example!

Walking through an example of the algorithm is particularly helpful to understand it:

As an example consider finding the next permutation of:

8342666411

The longest monotonically decreasing tail is 666411, and the corresponding head is 8342.

8342 666411

666411 is, by definition, reverse-ordered, and cannot be increased by permuting its elements. To find the next permutation, we must increase the head; a matter of finding the smallest tail element larger than the head’s final 2.

8342 666411

Walking back from the end of tail, the first element greater than 2 is 4.

8342 666411

Swap the 2 and the 4

8344 666211

Since head has increased, we now have a greater permutation. To reduce to the next permutation, we reverse tail, putting it into increasing order.

8344 112666

Join the head and tail back together. The permutation one greater than 8342666411 is 8344112666.

Implementation

We can use the steps in the example above to help organize our code. First, we find the "cut point" which is the index to the left of the longest monotonic tail:

: cut-point ( seq -- n )
    [ last ] keep [ [ > ] keep swap ] find-last drop nip ;

Next, we want to find the smallest element larger than the element at the "cut point" (searching from the end of the sequence):

: greater-from-last ( n seq -- i )
    [ nip ] [ nth ] 2bi [ > ] curry find-last drop ;

We then need a way to reverse the tail of a sequence (the ! denotes that this modifies the sequence in-place):

: reverse-tail! ( n seq -- seq )
    [ swap 1 + tail-slice reverse! drop ] keep ;

Putting this together gives us a word that looks a bit like our example. I have decided that in the case where the sequence reaches its lexicographic greatest order, we reverse it to its smallest ordering. This allows it to cycle through all possible permutations no matter where you start.

: (next-permutation) ( seq -- seq )
    dup cut-point [
        swap [ greater-from-last ] 2keep
        [ exchange ] [ reverse-tail! nip ] 3bi
    ] [ reverse! ] if* ;

We wrap this with a simple check to make sure the sequence is not empty. Arguably, we could instead check if the length is greater than 1 since a single element sequence has only possible permutation.

: next-permutation ( seq -- seq )
    dup [ ] [ drop (next-permutation) ] if-empty ;

Testing

Factor makes it easy to do unit tests. Here are some of the ones I've used to test this code:

[ "" ] [ "" next-permutation ] unit-test

[ "1" ] [ "1" next-permutation ] unit-test

[ "21" ] [ "12" next-permutation ] unit-test

[ "8344112666" ] [ "8342666411" next-permutation ] unit-test

[ "ABC" "ACB" "BAC" "BCA" "CAB" "CBA" "ABC" ]
[ "ABC" 6 [ dup >string next-permutation ] times ] unit-test

This is available in the latest development branch of Factor.

Bonus

Solving the original programming challenge is as easy as:

IN: scratchpad 38276 number>string next-permutation string>number .
38627

Tuesday, February 21, 2012

TXON

The TXON, also known as "Text Object Notation", is a proposed format for structured data.

Much less popular than other formats such as JSON, XML, or even INI files - I thought it would still be fun to implement encode and decode words in Factor.

An example TXON might look something like this:

Factor:`
    url:`http://factorcode.org`
    development:`Started in 2003`
    license:`Open source (BSD license)`
    influences:`Forth, Lisp, and Smalltalk`
`

Encoding

Since TXON uses "`" characters to delimit values, we need to escape them:

: encode-value ( string -- string' )
    R" `" "\\`" re-replace ;

To implement encoding in a generic way, we dispatch on the type of object being encoded:

GENERIC: >txon ( object -- string )

M: sequence >txon
    [ >txon ] map "\n" join ;

M: assoc >txon
    >alist [
        first2 [ encode-value ] [ >txon ] bi* "%s:`%s`" sprintf
    ] map "\n" join ;

M: string >txon
    encode-value ;

M: number >txon
    number>string >txon ;

Decoding

Although the TXON specification includes an EBNF grammar, I am going to show one way to build a parser from scratch. In the tradition of concatenative languages, we will build our decoder from several smaller words.

For symmetry with the encode-value word, we need a way to unescape the ` characters:

: decode-value ( string -- string' )
    R" \\`" "`" re-replace ;

Since the TXON format is a series of name:`value` pairs, we can parse the name by finding the separator and then decoding the name (which might contain escaped characters):

: parse-name ( string -- remain name )
    ":`" split1 swap decode-value ;

To build a word that finds the first (unescaped) ` character, we will first make a word that looks at adjacent characters, returning true if the second character is an unescaped `:

: `? ( ch1 ch2 -- ? )
    [ CHAR: \ = not ] [ CHAR: ` = ] bi* and ;

By grouping the string into adjacent characters, we can find the first unescaped ` (specially handling the case where the first character is an `):

: (find-`) ( string -- n/f )
    2 clump [ first2 `? ] find drop [ 1 + ] [ f ] if* ;

: find-` ( string -- n/f )
    dup ?first CHAR: ` = [ drop 0 ] [ (find-`) ] if ;

Parsing the value is slightly complicated by the fact that TXON supports values which might themselves be a single value, a sequence of values, or a series of name/value pairs. Basically, that means we need to:

  1. find the first ` character
  2. checks if the previous character is a : (indicating a name/value)
  3. parse all name/values if so, otherwise decode the value(s)

That algorithm can be translated into this code:

DEFER: name/values

: (parse-value) ( string -- values )
    decode-value string-lines dup length 1 = [ first ] when ;

: parse-value ( string -- remain value )
    dup find-` [
        dup 1 - pick ?nth CHAR: : =
        [ drop name/values ] [ cut swap (parse-value) ] if
        [ rest [ blank? ] trim-head ] dip
    ] [ f swap ] if* ;

We want to parse a "name=value" pair, which should be as easy as parsing the name, then the value, then associating into a hashtable:

: (name=value) ( string -- remain term )
    parse-name [ parse-value ] dip associate ;

The string might contain a "name=value" pair, or just a single value:

: name=value ( string -- remain term )
    [ blank? ] trim
    ":`" over subseq? [ (name=value) ] [ f swap ] if ;

We finish by building a word to produce all "name=value" pairs, used in the parse-value word earlier.

: name/values ( string -- remain terms )
    [ dup { [ empty? not ] [ first CHAR: ` = not ] } 1&& ]
    [ name=value ] produce assoc-combine ;

Putting all of that together, we can make a word to parse a TXON string, producing "name=value" pairs until exhausted:

: parse-txon ( string -- objects )
    [ dup empty? not ] [ name=value ] produce nip ;

: txon> ( string -- object )
    parse-txon dup length 1 = [ first ] when ;

Try It

You can try this out in the listener:

IN: scratchpad H{ { "a" "123" } } >txon .
"a:`123`"

IN: scratchpad "a:`123`" txon> .
H{ { "a" "123" } }

Can you improve on this? Maybe by using the peg.ebnf vocabulary to create an EBNF parsing word?

The code for this (and a bunch of tests) are on my Github.

Monday, February 13, 2012

Readability

James O'Beirne wrote a great blog post on why languages matter with some thoughts on predictability, readability, and compactness. In it, he compares some examples of code in dynamic languages such as PHP, Python, and Groovy.

I wanted to compare his simple "readability" examples with Factor, to show why concatenative programming matters.

PHP

This example of some PHP code:

<?php
$x = 1;
$nums = array(10, 20, 30, 40);
$res = 0;

foreach ($nums as $n)
  if ($n > 15)
    $res -= $n*2 + $x;

Groovy

He compares the PHP code favorably to this Groovy code:

def x = 1
def nums = [10, 20, 30, 40]
def res = nums.findAll { it > 15 } 
    .collect { it * 2 + x } 
    .inject(0) {accum, val -> accum - val}

Python

Although James doesn't show a Python example, it could look something like this:

x = 1
nums = [10,20,30,40]
res = 0
for n in nums:
    if n > 15:
        res -= (n*2) + x

We could improve this by using a generator expression (or, equivalently, a list comprehension)):

>>> def foo(x, nums):
...     return -sum(n*2+x for n in nums if n > 15)
... 

>>> foo(1, [10,20,30,40])
-183

Factor

I would argue that a Factor version is pretty readable:

IN: scratchpad { 10 20 30 40 } [ 15 > ] filter
               [ 2 * 1 + ] map sum neg .
-183

If you wanted to factor (ahem) this into a reusable word:

: foo ( x nums -- res )
    [ 15 > ] filter [ 2 * + ] with map sum neg ;

IN: scratchpad 1 { 10 20 30 40 } foo .
-183

Thursday, February 2, 2012

copy

I've used Factor to build several common unix programs including cat, fortune, wc, and others.

Today, I wanted to show how to build the cp ("copy") program using the simple file manipulation words. If we look at the man page, we can see that its usage is two-fold:

  1. Copy several source files to a destination directory
  2. Copy a source file to a destination file or directory

We can make a nice usage string to display if the arguments are not correct:

: usage ( -- )
    "Usage: copy source ... target" print ;

We can implement the first usage, copy-to-dir, by checking to see that the destination is a directory before calling copy-files-into, or printing the usage if it is not:

: copy-to-dir ( args -- )
    dup last file-info directory?
    [ unclip-last copy-files-into ] [ drop usage ] if ;

The second usage, copy-to-file, first checks if the destination exists and is a directory (if so calling our copy-to-dir word), otherwise calling copy-file:

: copy-to-file ( args -- )
    dup last { [ exists? ] [ file-info directory? ] } 1&&
    [ copy-to-dir ] [ first2 copy-file ] if ;

Putting it all together, we can implement our program by checking the number of arguments and assuming the two-argument version is copy-to-file and more arguments are copy-to-dir (anything less gets the usage):

: run-copy ( -- )
    command-line get dup length {
        { [ dup 2 > ] [ drop copy-to-dir  ] }
        { [ dup 2 = ] [ drop copy-to-file ] }
        [ 2drop usage ]
    } cond ;

MAIN: run-copy

The code for this is on my Github.