Saturday, September 19, 2015

ISBN

Most books are issued a unique International Standard Book Number (ISBN) number. Often different formats of the same book will have different ISBN numbers. On a print book, you can usually find the ISBN on a barcode on the back cover.

Most countries seem to have a national ISBN registration agency. In some countries this is a free service provided by a government agency. In other countries, this is operated by a commercial entity. In the United States, one company (R.R. Bowker LLC) has an apparent monopoly on issuing ISBN numbers which can cost $125 for one (less if you buy in bulk).

The ISBN is 13 digits long if assigned starting in 2007, and 10 digits long if assigned before 2007. Each ISBN contains a check digit which is used for basic error detection We are going to build a few words in Factor to calculate the check digits and validate ISBNs.

We need to turn an ISBN (which might include spaces or dashes) into numeric digits:

: digits ( str -- digits )
    [ digit? ] filter string>digits ;

For ISBN-10, the check digit is the sum of each of 10 digits multiplied by a weight (descending from 10 to 1) modulo 11.

: isbn-10-check ( digits -- n )
    0 [ 10 swap - * + ] reduce-index 11 mod ;

For ISBN-13, the check digit is the sum of each of 13 digits multiplied by a weight (alternating between 1 and 3) modulo 10.

: isbn-13-check ( digits -- n )
    0 [ even? 1 3 ? * + ] reduce-index 10 mod ;

We can validate an ISBN by grabbing the digits and running either the ISBN-10 or ISBN-13 check and verifying that the result is zero.

: valid-isbn? ( str -- ? )
    digits dup length {
        { 10 [ isbn-10-check ] }
        { 13 [ isbn-13-check ] }
    } case 0 = ;

The code (and some tests) for this is on my GitHub.

Friday, September 11, 2015

Pig Latin

Pig Latin is a somewhat ridiculous language game which modifies words in such a funny way that is hard to figure out if you don't know how it works but easy if you do. Using Factor, we will build a converter from English to Pig Latin words.

There are two basic rules we should implement:

  1. For words that begin with consonant sounds, the initial consonant or consonant cluster is moved to the end of the word, and "ay" is added to the end.
{ "igpay" } [ "pig" pig-latin ] unit-test
{ "ananabay" } [ "banana" pig-latin ] unit-test
{ "ashtray" } [ "trash" pig-latin ] unit-test
{ "appyhay" } [ "happy" pig-latin ] unit-test
{ "uckday" } [ "duck" pig-latin ] unit-test
{ "oveglay" } [ "glove" pig-latin ] unit-test
  1. For words that begin with a vowel sounds or silent letter, add "way" to the end.
{ "eggway" } [ "egg" pig-latin ] unit-test
{ "inboxway" } [ "inbox" pig-latin ] unit-test
{ "eightway" } [ "eight" pig-latin ] unit-test

We can implement our two basic rules:

: pig-latin ( str -- str' )
    dup [ "aeiou" member? ] find drop [
        "way" append
    ] [
        cut swap "ay" 3append
    ] if-zero ;

We could improve this by:

  • better handling of words that start with capital vowels or are all consonants
  • reverse the rules to convert Pig Latin back to English
  • variations such as adding "yay" (or "i") instead of "way"
  • different rules like adding "ag" before each vowel ("pagig lagatagin")
  • support language games in other languages

Anyway, this is available on my GitHub.

Sunday, August 30, 2015

Bowling Scores

Today we are going to explore building a bowling score calculator using Factor. In particular, we will be scoring ten-pin bowling.

There are a lot of ways to "golf" this, including this short version in F#, but we will build this in several steps through transformations of the input. The test input is a string representation of the hits, misses, spares, and strikes. The output will be a number which is your total score. We will assume valid inputs and not do much error-checking.

A sample game might look like this:

12X4--3-69/-98/8-8-

Our first transformation is to convert each character to a number of pins that have been knocked down for each ball. Strikes are denoted with X, spares with /, misses with -, and normal hits with a number.

: pin ( last ch -- pin )
    {
        { CHAR: X [ 10 ] }
        { CHAR: / [ 10 over - ] }
        { CHAR: - [ 0 ] }
        [ CHAR: 0 - ]
    } case nip ;

We use this to convert the entire string into a series of pins knocked down for each ball.

: pins ( str -- pins )
    f swap [ pin dup ] { } map-as nip ;

A single frame will be either one ball, if a strike, or two balls. We are going to use cut-slice instead of cut because it will be helpful later.

: frame ( pins -- rest frame )
    dup first 10 = 1 2 ? short cut-slice swap ;

A game is 9 "normal" frames and then a last frame that could have up to three balls in it.

: frames ( pins -- frames )
    9 [ frame ] replicate swap suffix ;

Some frames will trigger a bonus. Strikes add the value of the next two balls. Spares add the value of the next ball. We build this by "un-slicing" the frame and calling sum on the next balls.

: bonus ( frame -- bonus )
    [ seq>> ] [ to>> tail ] [ length 3 swap - ] tri head sum ;

We can score the frames by checking for frames where all ten pins are knocked down (either spares or strikes) and adding their bonus.

: scores ( frames -- scores )
    [ [ sum ] keep over 10 = [ bonus + ] [ drop ] if ] map ;

We can solve the original goal by just adding all the scores:

: bowl ( str -- score )
    pins frames scores sum ;

And write a bunch of unit tests to make sure it works:

{ 0 } [ "---------------------" bowl ] unit-test
{ 11 } [ "------------------X1-" bowl ] unit-test
{ 12 } [ "----------------X1-" bowl ] unit-test
{ 15 } [ "------------------5/5" bowl ] unit-test
{ 20 } [ "11111111111111111111" bowl ] unit-test
{ 20 } [ "5/5-----------------" bowl ] unit-test
{ 20 } [ "------------------5/X" bowl ] unit-test
{ 40 } [ "X5/5----------------" bowl ] unit-test
{ 80 } [ "-8-7714215X6172183-" bowl ] unit-test
{ 83 } [ "12X4--3-69/-98/8-8-" bowl ] unit-test
{ 150 } [ "5/5/5/5/5/5/5/5/5/5/5" bowl ] unit-test
{ 144 } [ "XXX6-3/819-44X6-" bowl ] unit-test
{ 266 } [ "XXXXXXXXX81-" bowl ] unit-test
{ 271 } [ "XXXXXXXXX9/2" bowl ] unit-test
{ 279 } [ "XXXXXXXXXX33" bowl ] unit-test
{ 295 } [ "XXXXXXXXXXX5" bowl ] unit-test
{ 300 } [ "XXXXXXXXXXXX" bowl ] unit-test
{ 100 } [ "-/-/-/-/-/-/-/-/-/-/-" bowl ] unit-test
{ 190 } [ "9/9/9/9/9/9/9/9/9/9/9" bowl ] unit-test

This is available on my GitHub.

Wednesday, August 26, 2015

Haikunator

The Haikunator is a project to provide "Heroku-like memorable random names". These names usually consist of an adjective, a noun, and a random number or token. The original repository is implemented in Ruby, with ports to Go, Javascript, Python, PHP, Elixer, .NET, Java, and Dart.

We will be implementing this in Factor using the qw vocabulary that provides a simple way to make "arrays of strings" using the qw{ syntax.

First, a list of adjectives:

CONSTANT: adjectives qw{
    autumn hidden bitter misty silent empty dry dark summer icy
    delicate quiet white cool spring winter patient twilight
    dawn crimson wispy weathered blue billowing broken cold
    damp falling frosty green long late lingering bold little
    morning muddy old red rough still small sparkling throbbing
    shy wandering withered wild black young holy solitary
    fragrant aged snowy proud floral restless divine polished
    ancient purple lively nameless lucky odd tiny free dry
    yellow orange gentle tight super royal broad steep flat
    square round mute noisy hushy raspy soft shrill rapid sweet
    curly calm jolly fancy plain shinny
}

Next, a list of nouns:

CONSTANT: nouns qw{
    waterfall river breeze moon rain wind sea morning snow lake
    sunset pine shadow leaf dawn glitter forest hill cloud
    meadow sun glade bird brook butterfly bush dew dust field
    fire flower firefly feather grass haze mountain night pond
    darkness snowflake silence sound sky shape surf thunder
    violet water wildflower wave water resonance sun wood dream
    cherry tree fog frost voice paper frog smoke star atom band
    bar base block boat term credit art fashion truth disk
    math unit cell scene heart recipe union limit bread toast
    bonus lab mud mode poetry tooth hall king queen lion tiger
    penguin kiwi cake mouse rice coke hola salad hat
}

We will make a token out of digits:

CONSTANT: token-chars "0123456789"

Finally, a simple haikunate implementation:

: haikunate ( -- str )
    adjectives random
    nouns random
    4 [ token-chars random ] "" replicate-as
    "%s-%s-%s" sprintf ;

We can try it a few times, to see how it works:

IN: scratchpad haikunate .
"odd-water-8344"

IN: scratchpad haikunate .
"flat-tooth-9324"

IN: scratchpad haikunate .
"wandering-lion-8346"

IN: scratchpad haikunate .
"yellow-mud-9780"

IN: scratchpad haikunate .
"patient-unit-4203"

IN: scratchpad haikunate .
"floral-feather-1023"

Some versions of "haikunate" in other languages include features such as:

  • allow customization of the delimiter (dots are popular)
  • allow the token to be specified as a range of possible numbers
  • allow the token to be restricted to a maximum length
  • allow the token to be represented using hex digits
  • allow the token to be represented with custom character sets
  • etc.

This is available on my GitHub.

Monday, August 17, 2015

Random Desktop Background

As a follow-up to my Desktop Background post, I wanted to show how to set your desktop background to random images from various online image sites. We will be downloading a URL to a local file and then setting the desktop picture to that file:

: download-and-set-desktop-picture ( url -- )
    dup "/" split1-last nip cache-file
    [ download-to ] [ set-desktop-picture ] bi ;

Okay, now we need some random images!

Imgur is a huge image-hosting website frequently used on sites like Reddit.

: random-imgur ( -- url )
    "https://imgur.com/random" scrape-html nip
    "image_src" "rel" find-by-attribute-key-value
    first "href" attribute ;

XKCD has some fun comics. Maybe they would look good on the desktop!

: random-xkcd ( -- url )
    "http://dynamic.xkcd.com/random/comic/" http-get nip
    R@ http://imgs\.xkcd\.com/comics/[^\.]+\.(png|jpg)@
    first-match >string ;

WallpaperStock has a bunch of more traditional desktop images. We will scrape their random wallpaper page, find the first wallpaper thumbnail, load that wallpaper page, find the default image page, and then load that to find the image URL.

: random-wallpaperstock ( -- url )
    "http://wallpaperstock.net/random-wallpapers.html"
    scrape-html nip "wallpaper_thumb" find-by-class-between
    "a" find-by-name nip "href" attribute
    "http://wallpaperstock.net" prepend scrape-html nip
    "the_view_link" find-by-id nip "href" attribute
    "http:" prepend scrape-html nip "myImage" find-by-id nip
    "src" attribute "http:" prepend ;

Using this is as easy as:

IN: scratchpad random-imgur
               download-and-set-desktop-picture

IN: scratchpad random-xkcd
               download-and-set-desktop-picture

IN: scratchpad random-wallpaperstock
               download-and-set-desktop-picture

This is available on my GitHub.

Friday, August 14, 2015

Desktop Background

One of the benefits of learning to program is learning how to automate tasks performed with a computer. I thought it might be fun to build a simple vocabulary to allow getting and setting of the desktop background picture. Since Factor makes it pretty easy to build cross-platform vocabularies, we will implement this on Mac OS X, Linux, and Windows.

Our API consists of two words, one that gets the current desktop picture and one that sets a new desktop picture, dispatching based on which operating system we are running on (technically, based on the value of the os variable).

HOOK: get-desktop-picture os ( -- path )

HOOK: set-desktop-picture os ( path -- )

Mac OS X

On Mac OS X, we use AppleScript to ask the Finder what the path to the current desktop picture is, or tell it to set the desktop picture to a specific path.

M: macosx get-desktop-picture
    {
        "osascript" "-e"
        "tell app \"Finder\" to get posix path of (get desktop picture as alias)"
    } utf8 [ readln ] with-process-reader ;

M: macosx set-desktop-picture
    absolute-path
    "tell application \"Finder\" to set desktop picture to POSIX file \"%s\""
    sprintf run-apple-script ;

Windows

On Windows, we use the SystemParametersInfo function to get and set the desktop wallpaper.

CONSTANT: SPI_GETDESKWALLPAPER 0x0073

CONSTANT: SPI_SETDESKWALLPAPER 0x0014

M: windows get-desktop-picture
    SPI_GETDESKWALLPAPER MAX_PATH dup 1 + WCHAR <c-array> [
        0 SystemParametersInfo win32-error<>0
    ] keep alien>native-string ;

M: windows set-desktop-picture
    [ SPI_SETDESKWALLPAPER 0 ] dip utf16n encode
    0 SystemParametersInfo win32-error<>0 ;

Linux

On Linux, which has many different desktops, we are going to assume a GNOME environment. Other window managers have different ways to change the desktop background.

M: linux get-desktop-picture
    {
        "gsettings"
        "get"
        "org.gnome.desktop.background"
        "picture-uri"
    } utf8 [ readln ] with-process-reader
    "'file://" ?head drop "'" ?tail drop ;

M: linux set-desktop-picture
    {
        "gsettings"
        "set"
        "org.gnome.desktop.background"
        "picture-uri"
    } swap absolute-path "file://" prepend suffix try-process ;

This is available on my GitHub.

Thursday, August 6, 2015

Automated Reasoning

There was a post about Automated Reasoning in F#, Scala, Haskell, C++, and Julia that uses a simple algorithm from John Harrison's book Handbook of Practical Logic and Automated Reasoning to simplify this equation:

e = (1 + (0 * x)) * 3) + 12

Factor has support for ML-style pattern matching and I thought it would be fun to contribute a simple solution using the match vocabulary.

We want to define a few types of expressions:

TUPLE: Var s ;
TUPLE: Const n ;
TUPLE: Add x y ;
TUPLE: Mul x y ;
Note: we could have made this simpler by assuming integers are constants and strings are variables rather than define the Const and Var tuples, but I wanted to keep this close to the code in the original blog post.

To be able to pattern match, we need to define some match variables:

MATCH-VARS: ?x ?y ;

We want a way to do a single simplification of an expression:

: simplify1 ( expr -- expr' )
    {
        { T{ Add f T{ Const f 0 } ?x } [ ?x ] }
        { T{ Add f ?x T{ Const f 0 } } [ ?x ] }
        { T{ Mul f ?x T{ Const f 1 } } [ ?x ] }
        { T{ Mul f T{ Const f 1 } ?x } [ ?x ] }
        { T{ Mul f ?x T{ Const f 0 } } [ T{ Const f 0 } ] }
        { T{ Mul f T{ Const f 0 } ?x } [ T{ Const f 0 } ] }
        { T{ Add f T{ Const f ?x } T{ Const f ?y } }
                                       [ ?x ?y + Const boa ] }
        { T{ Mul f T{ Const f ?x } T{ Const f ?y } }
                                       [ ?x ?y * Const boa ] }
        [ ]
    } match-cond ;

We have a way to recursively simplify some expressions:

: simplify ( expr -- expr' )
    {
        { T{ Add f ?x ?y } [ ?x ?y [ simplify ] bi@ Add boa ] }
        { T{ Mul f ?x ?y } [ ?x ?y [ simplify ] bi@ Mul boa ] }
        [ ]
    } match-cond simplify1 ;

Finally, we have a word that tries to simplify a value to a constant:

: simplify-value ( expr -- str )
    simplify {
        { T{ Const f ?x } [ ?x ] }
        [ drop "Could not be simplified to a Constant." ]
    } match-cond ;

To check that it works, we can write a unit test that simplifies the original expression above:

{ 15 } [
    T{ Add f
        T{ Mul f
            T{ Add f
                T{ Const f 1 }
                T{ Mul f
                    T{ Const f 0 }
                    T{ Var f "x" } } }
            T{ Const f 3 } }
        T{ Const f 12 } }
    simplify-value
] unit-test

That's cool, but wouldn't it be better if we could work on quotations directly? Let's make a word that converts a quotation to an expression:

: >expr ( quot -- expr )
    [
        {
            { [ dup string? ] [ '[ _ Var boa ] ] }
            { [ dup integer? ] [ '[ _ Const boa ] ] }
            { [ dup \ + = ] [ drop [ Add boa ] ] }
            { [ dup \ * = ] [ drop [ Mul boa ] ] }
        } cond
    ] map concat call( -- expr ) ;

Now that we have that, our test case is a lot simpler:

{ 15 } [
    [ "x" 0 * 1 + 3 * 12 + ] >expr simplify-value
] unit-test

The code for this is on my GitHub.

Note: this takes advantage of a small feature that I added to the match-cond word to provide a way to easily handle a fall-through pattern like the cond word.

Friday, June 19, 2015

Bit Test

Factor has a bit? generic that is used to test an integer to see if a particular bit is set. When operating on fixnum integers, this is implemented by the fixnum-bit? word:

: fixnum-bit? ( x n -- ? )
    { fixnum fixnum } declare
    dup 0 >= [ neg shift even? not ] [ 2drop f ] if 

Many CPU architectures have instructions for performing a Bit Test. On x86, the BT instruction is available. Below, we are going to implement a compiler intrinsic in the hopes of speeding up this operation in Factor.

Intrinsic

In cpu.architecture, add a generic %bit-test based on our CPU:

HOOK: %bit-test cpu ( dst src1 src2 temp -- )

In cpu.x86, implement %bit-test on x86 (returning a boolean using a temporary register and the CMOVB instruction to check the carry flag which holds the result of the bit test):

M:: x86 %bit-test ( dst src1 src2 temp -- )
    src1 src2 BT
    dst temp \ CMOVB (%boolean) ;

In compiler.cfg.instructions, add a ##bit-test instruction:

FOLDABLE-INSN: ##bit-test
def: dst/tagged-rep
use: src1/int-rep src2/int-rep
temp: temp/int-rep ;

In compiler.codegen, we link the ##bit-test instruction with the %bit-test word.

CODEGEN: ##bit-test %bit-test

In compiler.cfg.intrinsics, enable replacing fixnum-bit? with the %bit-test intrinsic:

: enable-bit-test ( -- )
    {
        { fixnum-bit? [ drop [ ^^bit-test ] binary-op ] }
    } enable-intrinsics ;

Disassemble

The old assembly using fixnum-bit? looked like this:

000000010ea00250: mov [rip-0x1ff256], eax
000000010ea00256: sub rsp, 0x8
000000010ea0025a: call 0x10eedf3b0 (integer>fixnum-strict)
000000010ea0025f: mov rax, [r14]
000000010ea00262: cmp rax, 0x0
000000010ea00266: jl 0x10ea002a7 (M\ fixnum bit? + 0x57)
000000010ea0026c: neg rax
000000010ea0026f: mov [r14], rax
000000010ea00272: call 0x10ebec060 (fixnum-shift)
000000010ea00277: mov rax, [r14]
000000010ea0027a: test rax, 0x10
000000010ea00281: mov rax, 0x1
000000010ea0028b: mov rbx, 0x1010eff4c
000000010ea00295: cmovnz rax, rbx
000000010ea00299: mov [r14], rax
000000010ea0029c: mov [rip-0x1ff2a2], eax
000000010ea002a2: add rsp, 0x8
000000010ea002a6: ret
000000010ea002a7: sub r14, 0x8
000000010ea002ab: mov qword [r14], 0x1
000000010ea002b2: mov [rip-0x1ff2b8], eax
000000010ea002b8: add rsp, 0x8
000000010ea002bc: ret

The new assembly looks like this with BT:

000000010e656ec0: mov [rip-0x293ec6], eax
000000010e656ec6: sub rsp, 0x8
000000010e656eca: call 0x10e467ea0 (integer>fixnum-strict)
000000010e656ecf: mov rax, [r14]
000000010e656ed2: mov rbx, [r14-0x8]
000000010e656ed6: sar rbx, 0x4
000000010e656eda: sar rax, 0x4
000000010e656ede: bt rbx, rax
000000010e656ee2: mov rbx, 0x1
000000010e656eec: mov rcx, 0x1018f169c
000000010e656ef6: cmovb rbx, rcx
000000010e656efa: sub r14, 0x8
000000010e656efe: mov [r14], rbx
000000010e656f01: mov [rip-0x293f07], eax
000000010e656f07: add rsp, 0x8
000000010e656f0b: ret

Performance

Turns out that this speeds up fixnum-bit? by a lot. Using a simple test case:

: bench-fixnum-bit ( x n -- ? )
    { fixnum fixnum } declare
    [ 0 100,000,000 ] 2dip
    '[ _ _ bit? [ 1 + ] when ] times ;

Our old version was decently fast:

IN: scratchpad gc [ 0b101010101 0 bench-fixnum-bit ] time
Running time: 0.821433838 seconds

But our new version is much faster!

IN: scratchpad gc [ 0b101010101 0 bench-fixnum-bit ] time
Running time: 0.239439108 seconds

This has been committed to the development version and is available now.

Wednesday, June 17, 2015

Compressed Sieve of Eratosthenes

In my previous post, we used "2-3-5-7" wheel factorization to produce a faster Sieve of Eratosthenes with a compression factor of 2 (by storing only the odd numbers). Samuel Tardieu reminded me that the sieve variation used by Factor uses less memory with a higher compression factor of 3.75.

I thought would be a fun implementation to show below.

Version 6

The "2-3-5" wheel with the first 30 numbers that are not divisible by 2, 3, or 5.

CONSTANT: wheel-2-3-5 $[
    30 [1,b] [
        { 2 3 5 } [ divisor? ] with any? not
    ] B{ } filter-as
]
Note: We use 30 here because the "{ 2 3 5 } product" cycle is 30.

If you look at the wheel, you see that there are 8 candidates in every 30 numbers:

IN: scratchpad wheel-2-3-5 .
{ 1 7 11 13 17 19 23 29 }

The number 8 is interesting because it suggests that we can use a single byte to store the 8 candidates in each block of 30 numbers, using a bitmask for possible primes or f:

CONSTANT: masks $[
    30 [0,b) [
        wheel-2-3-5 index [ 7 swap - 2^ ] [ f ] if*
    ] map
]

We can index into this array to find the byte/mask for every number:

: byte-mask ( n -- byte mask/f )
    30 /mod masks nth ;

Using the byte/mask we can check if a number is marked in the sieve byte-array.

:: marked? ( n sieve -- ? )
    n byte-mask :> ( byte mask )
    mask [ byte sieve nth mask bitand zero? not ] [ f ] if ;

And we can unmark a number in the sieve in a similar fashion:

:: unmark ( n sieve -- )
    n byte-mask :> ( byte mask )
    mask [ byte sieve [ mask bitnot bitand ] change-nth ] when ;

Using that, we can unmark multiples of prime numbers:

:: unmark-multiples ( i upper sieve -- )
    i sq upper i 2 * <range> [ sieve unmark ] each ;

That's all we need to build our sieve, starting from a byte-array initialized with all numbers marked and then progressively unmarking multiples of primes:

:: sieve6 ( n -- sieve )
    n 30 /i 1 + [ 0xff ] B{ } replicate-as :> sieve
    sieve length 30 * 1 - :> upper
    3 upper sqrt 2 <range> [| i |
        i sieve marked? [
            i upper sieve unmark-multiples
        ] when
    ] each sieve ;

As you can see, storing prime numbers up to 10 million takes only 333 KB!

IN: scratchpad 10,000,000 sieve6 length .
333334

This approach is used by the math.primes vocabulary to quickly check primality of any number less than 9 million using only 300 KB of memory.

Neat!

Monday, June 8, 2015

Genuine Sieve of Eratosthenes

Iain Gray posted on the mailing list about a paper called the Genuine Sieve of Eratosthenes by Melissa O'Neill, a Professor of Computer Science at Harvey Mudd College.

It begins with a discussion about some clever looking Haskell code that is just trial division and not the Sieve of Eratosthenes. At the end of the paper is an interesting discussion of performance improvements that I thought might be fun to implement in Factor as a followup to the three versions that I posted about recently.

Version 4

We are going to use wheel factorization as a way of reducing the search space by ignoring some multiples of small prime numbers. It is a bit similar to ignoring all even numbers in our previous "Version 3". In this case, we want to ignore all multiples of the first four prime numbers.

To calculate the "2-3-5-7" wheel, we start with 11 and filter the next 210 numbers that are not divisible by 2, 3, 5, or 7.

CONSTANT: wheel-2-3-5-7 $[
    11 dup 210 + [a,b] [
        { 2 3 5 7 } [ divisor? ] with any? not
    ] B{ } filter-as differences
]
Note: We use 210 here because the "{ 2 3 5 7 } product" cycle is 210.

Next, we want a way to iterate across all the numbers that might be prime, calling a quotation on each number that is prime.

:: each-prime ( upto sieve quot -- )
    11 upto '[ dup _ <= ] [
        wheel-2-3-5-7 [
            over dup 2/ sieve nth [ drop ] quot if
            +
        ] each
    ] while drop ; inline

For each prime found, we will want to mark all odd multiples as not prime.

:: mark-multiples ( i upto sieve -- )
    i sq upto i 2 * <range> [| j |
        t j 2/ sieve set-nth
    ] each ; inline

We will use a bit-array that is sized in 210-bit blocks, so the number of bits needed is:

: sieve-bits ( n -- bits )
    210 /i 1 + 210 * 2/ 6 + ; inline

Finally, we calculate our sieve, first for 11 and then for all the primes above 11.

:: sieve4 ( n -- primes )
    n sieve-bits <bit-array> :> sieve
    t 0 sieve set-nth
    t 4 sieve set-nth
    n sqrt sieve [ n sieve mark-multiples ] each-prime
    V{ 2 3 5 7 } clone :> primes
    n sieve [
        dup n <= [ primes push ] [ drop ] if
    ] each-prime primes ;

Calculating prime numbers up to 10 million is getting even faster.

IN: scratchpad gc [ 10,000,000 sieve4 ] time
Running time: 0.08523477 seconds

Version 5

If you use the compiler.tree.debugger to inspect the optimized output, you can see some dynamic dispatch with branches for fixed-width (fixnum) and arbitrary-precision (bignum) integers.

It would be nice to perform fixnum specialization or improve type propagation in the compiler, but for now we are going to write some ugly code to force fixnum math for performance.

We make a version of each-prime that inlines the iteration:

:: each-prime2 ( upto sieve quot -- )
    11 upto >fixnum '[ dup _ <= ] [
        wheel-2-3-5-7 [
            over dup 2/ sieve nth-unsafe [ drop ] quot if
            fixnum+fast
        ] each
    ] while drop ; inline

And similarly for mark-multiples:

:: mark-multiples2 ( i upto sieve -- )
    i 2 fixnum*fast :> step
    i i fixnum*fast upto >fixnum '[ dup _ <= ] [
        t over 2/ sieve set-nth-unsafe
        step fixnum+fast
    ] while drop ; inline

And use them in our sieve computation:

:: sieve5 ( n -- primes )
    n sieve-bits <bit-array> :> sieve
    t 0 sieve set-nth
    t 4 sieve set-nth
    n sqrt sieve [ n sieve mark-multiples2 ] each-prime2
    V{ 2 3 5 7 } clone :> primes
    n sieve [
        dup n <= [ primes push ] [ drop ] if
    ] each-prime2 primes ;

Calculating prime numbers up to 10 million is super fast!

IN: scratchpad gc [ 10,000,000 sieve5 ] time
Running time: 0.033914378 seconds

We can compute all the 11,078,937 prime numbers up to 200 million in about a second!

IN: scratchpad gc [ 200,000,000 sieve5 ] time
Running time: 0.970876855 seconds

And all the 50,847,534 prime numbers up to 1 billion in about six seconds!

IN: scratchpad gc [ 1,000,000,000 sieve5 ] time
Running time: 6.141224805 seconds

Not quite primegen or primesieve, but not bad for a laptop!

This is available in the math.primes.erato.fast vocabulary on our nightly builds!

Sunday, June 7, 2015

Sieve of Eratosthenes

I noticed that the Crystal Programming Language homepage has an example of the Sieve of Eratosthenes algorithm for finding prime numbers. I thought it would be fun to show a few simple variations of that algorithm in Factor.

The algorithm works by iteratively marking as not prime the multiples of each prime, starting with the multiples of 2. In pseudocode, it looks like this:

Input: an integer n > 1
 
Let A be an array of Boolean values, indexed by integers 2 to n,
initially all set to true.
 
 for i = 2, 3, 4, ..., not exceeding n:
  if A[i] is true:
    for j = i2, i2+i, i2+2i, i2+3i, ..., not exceeding n :
      A[j] := false
 
Output: all i such that A[i] is true.

We will be looking at the performance and memory usage to calculate the 664,579 prime numbers between 1 and 10,000,000.

Version 1

Our first version, is a straightforward implementation from the pseudocode:

:: sieve1 ( n -- primes )
    n 1 + t <array> :> sieve
    f 0 sieve set-nth
    f 1 sieve set-nth

    2 n sqrt [a,b] [| i |
        i sieve nth [
            i sq n i <range> [| j |
                f j sieve set-nth
            ] each
        ] when
    ] each

    sieve [ drop ] selector [ each-index ] dip ;

We can show that it works for the first few prime numbers:

IN: scratchpad 25 sieve1 .
V{ 2 3 5 7 11 13 17 19 23 }

Calculating prime numbers up to 10 million takes less than half a second:

IN: scratchpad [ 10,000,000 sieve1 ] time
Running time: 0.408065579 seconds

The memory used is basically an array with 10 million elements in it, each element of the array is a boolean (64 bits on my laptop). Total memory used is 80 MB.

Version 2

Storing booleans can be more memory-efficient using bit-arrays, where each boolean is stored as a single bit. The logic is the same, however primes are marked by false rather than true.

:: sieve2 ( n -- primes )
    n 1 + <bit-array> :> sieve
    t 0 sieve set-nth
    t 1 sieve set-nth

    2 n sqrt [a,b] [| i |
        i sieve nth [
            i sq n i <range> [| j |
                t j sieve set-nth
            ] each
        ] unless
    ] each

    sieve [ drop not ] selector [ each-index ] dip ;

Calculating prime numbers up to 10 million is faster:

IN: scratchpad [ 10,000,000 sieve2 ] time
Running time: 0.241894524 seconds

The memory used is now one bit per number, so only 10 million bits or 1.25 MB.

Version 3

Since we know that even numbers (except 2) are not prime, we can only work with odd numbers and start our search from 3. For each number found, we can check only the odd multiples by marking off i2 + k*i for only even k.

:: sieve3 ( n -- sieve )
    n dup odd? [ 1 + ] when 2/ <bit-array> :> sieve
    t 0 sieve set-nth

    3 n sqrt 2 <range> [| i |
        i 2/ sieve nth [
            i sq n i 2 * <range> [| j |
                t j 2/ sieve set-nth
            ] each
        ] unless
    ] each

    V{ 2 } clone :> primes

    sieve [
        swap [ drop ] [ 2 * 1 + primes push ] if
    ] each-index

    primes ;
Note: the convenience of selector is not available because we want to initialize the list of primes with 2. Instead, we just do it manually.

Calculating prime numbers up to 10 million is much faster!

IN: scratchpad [ 10,000,000 sieve3 ] time
Running time: 0.110387127 seconds

As for memory usage, it is storing one bit for each odd number, so only 5 million bits or 625 KB.

Monday, June 1, 2015

SEND + MORE = MONEY?

There's a classic verbal arithmetic puzzle that was published in 1924 by Henry Dudeney, where you solve the equation:

    S E N D
+   M O R E
-----------
= M O N E Y
Note: Each letter is a unique digit and S and M can't be zero.

A few days ago, I noticed a blog post with solutions in Perl 6 that references another blog post with a solution in Haskell. I thought I'd show a solution in Factor using the backtrack vocabulary that provides something like John McCarthy's amb ambiguous operator.

First, we need a list of available digits, just the numbers 0 through 9:

CONSTANT: digits { 0 1 2 3 4 5 6 7 8 9 }

Next, we turn a sequence of digits into a number (e.g., { 1 2 3 4 } becomes 1234):

: >number ( seq -- n ) 0 [ [ 10 * ] dip + ] reduce ;

We can then implement our solver, by choosing digits while progressively restricting the possible values for future digits using the ones we've chosen so far (using local variables to store the digits):

:: send-more-money ( -- )
    [
        digits { 0 } diff amb-lazy :> s
        digits { s } diff amb-lazy :> e
        digits { s e } diff amb-lazy :> n
        digits { s e n } diff amb-lazy :> d
        { s e n d } >number :> send

        digits { 0 s e n d } diff amb-lazy :> m
        digits { s e n d m } diff amb-lazy :> o
        digits { s e n d m o } diff amb-lazy :> r
        { m o r e } >number :> more

        digits { s e n d m o r } diff amb-lazy :> y
        { m o n e y } >number :> money

        send more + money = [
            send more money "   %s\n+  %s\n= %s\n" printf
        ] when

    ] amb-all ;
Note: We search all solutions using amb-all (even though there is only one). In this case, it is effectively an iterative search, which we could implement without backtracking. If we wanted the first solution, we could use if-amb.

So, what's the answer? Let's see!

IN: scratchpad send-more-money
   9567
+  1085
= 10652

Neat! And it's fast too -- solving this puzzle in about 2.5 seconds on my laptop.

The code for this is on my GitHub.

Tuesday, May 5, 2015

File Monitor

Factor has a cross-platform file-system change monitor which supports detecting changes to file names, attributes and contents under a specified directory.

There is some minor platform differences between Mac OS X, Windows, and Linux which may be worth looking at if you are building on top of the io.monitors vocabulary. I was curious about what kind of events are generated for various test-cases and built a small utility to experiment with it.

Some code to monitor for changed paths recursively in a directory and print each one out:

: watch-loop ( monitor -- )
    dup next-change path>> print flush watch-loop ;

: watch-directory ( path -- )
    [ t [ watch-loop ] with-monitor ] with-monitors ;

I've committed this as the file-monitor tool (with support for an optional command-line argument to specify which directory to monitor as well as printing the change descriptors). You can run it very simply:

factor -run=file-monitor [path]

An example session on Linux, monitoring some simple changes to files in /tmp:

$ factor -run=file-monitor /tmp &
Monitoring /tmp

$ touch /tmp/a
{ +add-file+ } /tmp/a
{ +add-file+ } /tmp/a
{ +modify-file+ } /tmp/a

$ echo "test" > /tmp/a
{ +modify-file+ } /tmp/a

$ rm /tmp/a
{ +remove-file+ } /tmp/a

Friday, May 1, 2015

File Server

Python has a neat feature that lets you serve files from the current directory.

# Python 2
python2 -m SimpleHTTPServer

# Python 3
python3 -m http.server

I always thought this was a quick and useful way to share files on a local network. Given that Factor has a HTTP Server, we should be able to implement this!

We already have support for serving static content and serving CGI scripts, so we can very simply implement a script to create and launch a HTTP server for the current directory (or the one specified on the command-line), logging HTTP connections to stdout.

This is available in the file-server vocabulary, now you can:

factor -run=file-server [--cgi] [path]

Currently, this defaults to serving files on port 8080 from all available network interfaces. In the future, it would be nice to add the ability to specify port and network interfaces to bind.

Wednesday, April 29, 2015

Burrows-Wheeler Transform

The Burrows–Wheeler transform is a reversible method of rearranging text used to improve the performance of compression algorithms, such as bzip2.

We will implement transform, bwt, and inverse transform, ibwt, in both Python and Factor. First with a slow and simple algorithm, and then second with a faster version.

Version 1

We start with the pseudocode suggested in the Wikipedia article:

function BWT (string s)
   append an 'EOF' character to s
   create a table, rows are all possible rotations of s
   sort rows alphabetically
   return (last column of the table)

In Python, this might look like:

def bwt(s):
    s = s + '\0'
    n = len(s)
    m = sorted(s[i:] + s[:i] for i in range(n))
    return ''.join(x[-1] for x in m)

In Factor, using all-rotations, it might look like this:

: bwt ( seq -- seq' )
    0 suffix all-rotations natural-sort [ last ] map ;

The pseudocode to perform the inverse transform:

function inverseBWT (string s)
   create empty table 
       
   repeat length(s) times
       // first insert creates first column
       insert s as a column of table before first column
       sort rows of the table alphabetically
   return (row that ends with the 'EOF' character)

In Python, this might look like:

def ibwt(s):
    n = len(s)
    m = [''] * n
    for _ in range(n):
        m = sorted(s[i] + m[i] for i in range(n))
    return [x for x in m if x.endswith('\0')][0][:-1]

In Factor, we could implement it like this:

: ibwt ( seq -- seq' )
    [ length [ "" <array> ] keep ] keep
    '[ _ [ prefix ] 2map natural-sort ] times
    [ last 0 = ] find nip but-last ;

Unfortunately, this is very slow, with most of the performance loss in the invert transform.

Version 2

Another way to increase the speed of BWT inverse is to use an algorithm that returns an index into the sorted rotations along with the transform.

In Python, it looks like this:

def bwt(s):
    n = len(s)
    m = sorted(s[i:] + s[:i] for i in range(n))
    return m.index(s), ''.join(x[-1] for x in m)

In Factor, it might look like this:

: bwt ( seq -- i seq' )
    dup all-rotations natural-sort
    [ index ] [ [ last ] map ] bi ;

In Python, the inverse transform looks like this:

def ibwt(k, s):
    def row(k):
        permutation = sorted((t, i) for i, t in enumerate(s))
        for _ in s:
            t, k = permutation[k]
            yield t
    return ''.join(row(k))

In Factor, that roughly translates to:

: ibwt ( i seq -- seq' )
    [ length ] [ <enum> sort-values ] bi
    '[ _ nth first2 ] replicate nip ;

An improved version 2 is available in the development version. In particular, it uses rotated virtual sequences for increased performance and returns transformations that match the type of the input sequence.

Tuesday, April 28, 2015

Writing MIDI Files

Previously, I wrote about Reading MIDI Files using Factor.

Now, we are going to create a writer for MIDI files in less than 180 lines of additional code.

Variable-Length Quantity

To write a variable-length integer, we first "reverse" it, tagging the 8th bit of each additional byte. Then, we write each byte out to the output-stream.

: write-number ( n -- )
    [ 0x7f bitand ] keep

    [ -7 shift dup zero? ] [
        [ 8 shift ] dip
        [ 0x7f bitand 0x80 bitor + ] keep
    ] until drop

    [ [ -8 shift ] [ 7 bit? ] bi ]
    [ dup 0xff bitand write1 ] do while drop ;
Note: there is probably a cleaner way to do this. Patches are welcome! ☺

Text

Strings are encoded in UTF-8, prefixed with their encoded length in bytes (as a variable-length quantity).

: write-string ( str -- )
    utf8 encode [ length write-number ] [ write ] bi ;

Writing Events

The three types of events will each have to be handled differently. To do this, we will make a generic method that is given the previous status byte (to enable "running status" for MIDI events) and returns the new status byte.

GENERIC: write-event ( prev-status event -- status )

First, we write MIDI events, implementing the "running status".

: write-status ( prev-status status -- )
    dup 0xf0 < [
        [ = ] keep swap [ drop ] [ write1 ] if
    ] [
        nip write1
    ] if ;

: write-channel ( prev-status value status quot -- status )
    [
        swap [
            "channel" of + [ write-status ] keep
        ] keep
    ] dip call ; inline

M: midi-event write-event
    [ delta>> write-number ] [ value>> ] [ name>> ] tri {

        { "note-off" [
            0x80 [
                [ "note" of write1 ]
                [ "velocity" of write1 ] bi
            ] write-channel ] }
        { "note-on" [
            0x90 [
                [ "note" of write1 ]
                [ "velocity" of write1 ] bi
            ] write-channel ] }
        { "polytouch" [
            0xa0 [
                [ "note" of write1 ]
                [ "value" of write1 ] bi
            ] write-channel ] }
        { "control-change" [
            0xb0 [
                [ "control" of write1 ]
                [ "value" of write1 ] bi
            ] write-channel ] }
        { "program-change" [
            0xc0 [ "program" of write1 ] write-channel ] }
        { "aftertouch" [
            0xd0 [ "value" of write1 ] write-channel ] }
        { "pitchwheel" [
            0xe0 [
                "pitch" of min-pitchwheel -
                [ 0x7f bitand write1 ]
                [ -7 shift write1 ] bi
            ] write-channel ] }

        ! system common messages
        { "sysex" [
            [ drop 0xf0 dup write1 ] dip
            write 0xf7 write1 ] }
        { "quarter-made" [
            [ drop 0xf1 dup write1 ] dip
            [ "frame-type" of 4 shift ]
            [ "frame-value" of + ] bi write1 ] }
        { "songpos" [
            [ drop 0xf2 dup write1 ] dip
            [ 0x7f bitand write1 ]
            [ -7 shift write1 ] bi ] }
        { "song-select" [
            [ drop 0xf3 dup write1 ] dip write1 ] }
        { "tune-request" [ 2drop 0xf6 dup write1 ] }

        ! real-time messages
        { "clock" [ 2drop 0xf8 dup write1 ] }
        { "start" [ 2drop 0xfa dup write1 ] }
        { "continue" [ 2drop 0xfb dup write1 ] }
        { "stop" [ 2drop 0xfc dup write1 ] }
        { "active-sensing" [ 2drop 0xfe dup write1 ] }
        { "reset" [ 2drop 0xff dup write1 ] }
    } case ;

Next, we write meta events:

M: meta-event write-event
    [ delta>> write-number ] [ value>> ] [ name>> ] tri 
    0xff write1 {
        { "sequence-number" [
            B{ 0x00 0x02 } write 2 >be write ] }
        { "text" [ 0x01 write1 write-string ] }
        { "copyright" [ 0x02 write1 write-string ] }
        { "track-name" [ 0x03 write1 write-string ] }
        { "instrument-name" [ 0x04 write1 write-string ] }
        { "lyrics" [ 0x05 write1 write-string ] }
        { "marker" [ 0x06 write1 write-string ] }
        { "cue-point" [ 0x07 write1 write-string ] }
        { "device-name" [ 0x09 write1 write-string ] }
        { "channel-prefix" [ B{ 0x20 0x01 } write write1 ] }
        { "midi-port" [ B{ 0x21 0x01 } write write1 ] }
        { "end-of-track" [ B{ 0x2f 0x00 } write drop ] }
        { "set-tempo" [ B{ 0x51 0x03 } write 3 >be write ] }
        { "smpte-offset" [
            B{ 0x54 0x05 } write {
                [ "frame-rate" of 6 shift ]
                [ "hours" of + write1 ]
                [ "minutes" of write1 ]
                [ "seconds" of write1 ]
                [ "frames" of write1 ]
                [ "subframes" of write1 ]
            } cleave ] }
        { "time-signature" [
            B{ 0x58 0x04 } write {
                [ "numerator" of write1 ]
                [ "denominator" of 2 /i write1 ]
                [ "clocks-per-tick" of write1 ]
                [ "notated-32nd-notes-per-beat" of write1 ]
            } cleave ] }
        { "key-signature" [
            B{ 0x59 0x02 } write
            key-signatures value-at write ] }
        { "sequencer-specific" [
            0x7f write1
            [ length write-number ] [ write ] bi ] }
    } case drop f ;

Finally, we write system-exclusive events:

M: sysex-event write-event
    drop
    [ delta>> write-number ]
    [ type>> write1 ]
    [ bytes>> write ] tri f ;

Writing a MIDI header and tracks, generically as "chunks":

GENERIC: write-chunk ( chunk -- )

M: midi-header write-chunk
    $[ "MThd" >byte-array ] write
    $[ 6 4 >be ] write
    [ format>> ] [ #chunks>> ] [ division>> ] tri
    [ 2 >be write ] tri@ ;

M: midi-track write-chunk
    $[ "MTrk" >byte-array ] write
    binary [
        events>> f swap [ write-event ] each drop
    ] with-byte-writer
    [ length 4 >be write ] [ write ] bi ;

Finally, words to write MIDI objects, either to a byte-array, or to a file.

: write-midi ( midi -- )
    [ header>> write-chunk ]
    [ chunks>> [ write-chunk ] each ] bi ;

: midi> ( midi -- byte-array )
    binary [ write-midi ] with-byte-writer ;

: midi>file ( midi path -- )
    binary [ write-midi ] with-file-writer ;

This is available now in the midi vocabulary.