Smallissimo

Saturday, April 7, 2012

Meteor contest - Part 4 the SmallInteger

In previous step, we were at about 1,000,000 loops for only 10,000 productive.

And we have about a factor 100 to gain compared to C++.

I'm not sure we can be that clever and avoid any false solution.

Yes, sure, we cannot reach speed of C++ yet with our VM, even the Cog one.

But we can do something to reduce a single loop cost.

One thing to notice is that the 1st part generating all pieces positions do a lot of useless work. Until pieces reach the north edge, we repeat the same bit pattern every two rows. That's a clue indicating something is wrong.

Moreover, in 32 bit Squeak image, we only have 30 bits left for representing a positive SmallInteger.

That's exactly 6 rows. Then we'll start memory intensive LargePositiveInteger allocations, and more involving computations too.

Now the idea: instead of shifting the pieces position north, why not shift the board south every time the pair of bottom rows is full ?

Since we fill east to west and from south to north, this is doable and we then have to generate the possible pieces positions only for the first pair of row. For example

can be shifted two rows south

Since a piece can span over 5 rows, and we consider filling on row 1 and 2, the northern two rows are useless and the board requires only 6 rows, that is 30 bits, a positive SmallInteger...

Good. For the first row, we can use a specially crafted southPossiblePositions.

But how to handle north rows above row 6 and prevent them to spread out of north edge?

Well, we don't have to care of north. We can use the 180° symmetry and as soon as 6 rows is full, reverse the whole board and start at south again.

Of course, we still have to shift the solution for filling the board. To avoid creating a LargeInteger, we'll store the rowOffset shift in a new ShootoutMeteorPiece object, and use it in possible positions, and in solutions.

Object subclass: #ShootoutMeteorPiece
    instanceVariableNames: 'mask row'
    classVariableNames: ''
    poolDictionaries: ''
    category: 'Shootout'

mask: aPieceMask
mask := aPieceMask

forRow: rowOffset
row := rowOffset

mask
^mask

We also delegate whether the piece fits on the board and how to fill the compact string solution (with implicit knowledge of our board rotation trick above 6th row).

fitOnBoard: aBoardMask
^0 == (aBoardMask bitAnd: mask)

fillSolution: aString ncol: ncol withColor: c
    | offset |
    ^ row >= 6
        ifTrue:
            [offset := aString size + 1 - (row - 6 * ncol).
            mask bitsDo: [:k | aString at: offset - k put: c]]
        ifFalse:
            [offset := row * ncol.
            mask bitsDo: [:k | aString at: offset + k put: c]]

ShootoutMeteorBoard>>boardStringWithPieces: pArray
    | board |
    board := String new: ncell.
    1 to: pArray size do: [:i | | c |
        c := '0123456789*' at: i.
        (pArray at: i) fillSolution: board ncol: ncol withColor: c].
    ^board

One subtle thing is that we already handled the symmetry by removing half orientations of piece #6. If we rotate the board, we have to consider the rotated orientations of piece #6. Hmm, this will be a bit technical to keep this optimization but we'll do.

First, we must stop generating pieces positions above twoRows:

ShootoutMeteorBoard>>placesFor: aPieceMask do: aBlock
   | westMask eastMask cellNumber |
   eastMask := self shiftSEmost: aPieceMask.

   [[westMask := eastMask.
   [(cellNumber := westMask lowBit) > twoRows ifTrue: [^self].
   (self hasEastOrWestIsland: westMask) ifFalse: [aBlock value: westMask].
   self canShiftW: westMask] whileTrue: [westMask := self shiftW: westMask].
   self canShiftNE: eastMask] whileTrue: [eastMask := self shiftNE: eastMask].
   self canShiftNW: eastMask] whileTrue: [eastMask := self shiftNW: eastMask]

Let now start the massacre of our previous code.

ShootoutMeteorBoard>>initializePossiblePositions
   | positionsPerPiecePerCell |
   positionsPerPiecePerCell := pieces collect: [:aPieceMask |
        | possible iRot |
        possible := (Array new: twoRows) collect: [:freeCell | Array new: 12 withAll: (ShootoutMeteorPiece new mask: 0)].
        iRot := 0.
        self rotationsOf: aPieceMask do: [:rotated |
            iRot := iRot + 1.
            self placesFor: rotated do: [:shifted |
                (possible at: shifted lowBit) at: iRot put: (ShootoutMeteorPiece new mask:
                    ((self hasEastOrWestIsland: shifted) ifTrue: [0] ifFalse: [shifted]))]].
        possible].
   positionsPerPiece := (1 to: 4) collect: [:rowStatus |
        (1 to: twoRows - (rowStatus \\ 2 * ncol)) collect: [:cellNumber |
            (1 to: pieces size)
                collect: [:pieceNumber |
                    | possible |
                    possible := (positionsPerPiecePerCell at: pieceNumber) at: cellNumber.
                     ((pieceNumber = 6)
                        ifTrue: [possible copyFrom: rowStatus // 3 * 6 + 1 to: rowStatus // 3 * 6 + 6]
                        ifFalse: [possible])
                            reject: [:aPiece | aPiece mask = 0
                                or: [rowStatus odd and: [self hasSouthIsland: aPiece mask]]]]]].

Here the rowStatus will be 1 on first row, 2 from row 2 to 6, 3 at row 10 (after we just rotated the board), 4 from rows 9 down to 7.

Notice also that we prefer to handle the cellNumber index first, before the pieceIndex, as that fits better our solve message.

Of course, since we moved the 6th piece trick above, we remove it from generation of rotations

rotationsOf: aPieceMask do: aBlock
    | next |
    aBlock value: (next := aPieceMask); value: (self flip: next).
    5 timesRepeat: [aBlock value: (next := self rotate: next); value: (self flip: next)]

Now that we lost the northEdge, why bother filling toward south ?

We will have to search islands only at south

hasSouthIsland: aMask
    ^(self hasInsetZero: (southEdge bitAnd: aMask))
    or: [(self hasCornerIsland: aMask edge: southEdge edge: eastEdge)
    or: [(self hasCornerIsland: aMask edge: southEdge edge: westEdge)]]

A bit more clean-ups to stop masks at sixRows:

fromString: aString
    | rawString |
    rawString := aString reject: [:e | e isSeparator].
    ncell := rawString size.
    ncol := (aString readStream upTo: Character cr) count: [:e | e isSeparator not].
    twoRows := ncol * 2.
    sixRows := ncol * 6.
    self initializeRowColMasks.
    pieces := rawString asSet asSortedArray collect: [:char |
        self shiftSEmost:
            (rawString inject: 0 into: [:pmask :c | pmask * 2 + (c = char ifTrue: [1] ifFalse: [0])])].
    self initializePossiblePositions.

initializeRowColMasks
    southEdge := 1 << ncol - 1.
    southToNorthMasks := (0 to: 5) collect: [:i | southEdge << (ncol * i)].
    eastEdge := 1<<sixRows-1/southEdge.
    eastToWestMasks := (0 to: ncol - 1) collect: [:i | eastEdge << i].
    westEdge := eastToWestMasks last.
    oddRowsMask := 1<<sixRows-1/(1<<twoRows-1)*southEdge.
    evenRowsMask := oddRowsMask << ncol.
    northWestMask := westEdge bitAnd: evenRowsMask.
    northEastMask := eastEdge bitAnd: oddRowsMask.
   southWestMask := southEdge bitOr: (westEdge bitAnd: evenRowsMask).
   southEastMask := southEdge bitOr: (eastEdge bitAnd: oddRowsMask).

printMask: aPieceMask
^String new: ncol * 2 + 1 * 6 streamContents: [:aStream | self printMask: aPieceMask on: aStream]

printMask: aPieceMask on: aStream
    6 to: 1 by: -1 do: [:irow |
        | mask |
        irow odd ifTrue: [aStream space].
        mask := 1 << (southToNorthMasks at: irow) highBit.
        1 to: ncol do: [:icol |
            aStream
                nextPut: ((aPieceMask bitAnd: (mask := mask >> 1)) = 0
                    ifTrue: [$0]
                    ifFalse: [$1]);
                space].
        aStream cr]

We are ready to write our new solving algorithm:

searchPuzzlesWithColorMask: colorMask boardMask: bMask rowOffset: rowOff pieces: pArray ifFound: solutionBlock
    | nextFreeCell possibles colorBit iRow boardMask rowStatus |
    colorMask = 0 ifTrue: [ ^solutionBlock value: (self boardStringWithPieces: pieces) ].
    loopCount := loopCount + 1.
    boardMask := bMask.
    iRow := rowOff.
    [(nextFreeCell := (boardMask + 1) lowBit) > twoRows]
        whileTrue:
            [boardMask := (iRow := iRow + 2) = 6
                ifTrue: [boardMask bitReverse: sixRows]
                ifFalse: [ boardMask >> twoRows]].
    rowStatus :=
        (iRow < 6 ifTrue: [0] ifFalse: [2]) +
        (((iRow = 0 or: [iRow = 6]) and: [nextFreeCell <= ncol]) ifTrue: [1] ifFalse: [2]).
    possibles := (positionsPerPiece at: rowStatus) at: nextFreeCell.
    colorBit := 1.
    1 to: pieces size do: [:pieceNumber |
        (colorMask bitAnd: colorBit) = 0
            ifFalse: [(possibles at: pieceNumber) do: [:aPiece |
                (aPiece fitOnBoard: boardMask)
                    ifTrue:
                        [pieces at: pieceNumber put: (aPiece forRow: iRow).
                        self searchPuzzlesWithColorMask: colorMask - colorBit boardMask: boardMask + aPiece mask rowOffset: iRow pieces: pArray ifFound: solutionBlock]]].
        colorBit := colorBit * 2].
    ^nil

solvedPuzzleDo: solutionBlock
    loopCount := 0.
    self searchPuzzlesWithColorMask: 1 << pieces size - 1 boardMask: 0 rowOffset: 0 pieces: pieces copy ifFound: solutionBlock.
    ^loopCount

That is a bit more complex right now, that's the tribute to optimizations. Especially the rowStatus management cries for a refactoring. But right now, we just focus on CPU efficiency, not source sustainability.

Note that bitReverse: may be absent from a Pharo image. Just pick it in Squeak.

What is interesting is to spy what occurs at bit reversal

Just before the first reversal:

0 0 0 0 0
0 0 0 1 0
0 0 0 0 1
0 1 1 0 1
1 1 1 1 1
1 1 1 1 1
And just after:

1 1 1 1 1
1 1 1 1 1
1 0 1 1 0
1 0 0 0 0
0 1 0 0 0
0 0 0 0 0
The now two north-most rows are filled with ones, and it's on purpose that we did not shift before reversing, to keep a barrier toward north.

What about the performance?

MessageTally spyOn: [ShootoutMeteorBoard solveDefault].
->1218297 loops

That's a 20% more loops than Part3, but that's possible, we didn't use the same order for filling the board, especially after the rotation we started back at southEast in a rather less constrained area, and we delayed detection of bad boards.

But now, the cost is about 1.7 seconds instead of 7.

Great! This would place Smalltalk at a not so ridiculous rank among those shootout languages.

And we still have room to reduce the number of combinations.

The code is at http://ss3.gemstone.com/ss/Shootout/Shootout.blog-nice.5.mcz

Friday, April 6, 2012

Meteor Contest Part 3 - reducing combinations

So our naive algorithm costs quite a lot. The analysis of MessageTally tells it is not in the initialization phase, but rather in the solving phase. Let's measure how many iterations we performed:

solvedPuzzleDo: solutionBlock
^self searchPuzzlesWithPermutation: (1 to: pieces size) asArray rank: 1 mask: 0 pieces: pieces copy ifFound: solutionBlock

searchPuzzlesWithPermutation: perms rank: i mask: boardMask pieces: pArray ifFound: solutionBlock
    | index boardLowBit count |
    count := 1.
    i > perms size
        ifTrue:
            [solutionBlock value: (self boardStringWithPieces: pArray).
            ^count].
    boardLowBit := (boardMask + 1) lowBit.
    i to: perms size do: [:j |
        index := perms at: j.
        perms swap: i with: j.
        ((positionsPerPiece at: index) at: boardLowBit) do: [:rotMask |
            0 = (boardMask bitAnd: rotMask)
                ifTrue:
                    [pArray at: index put: rotMask.
                    count := count + (self searchPuzzlesWithPermutation: perms rank: i + 1 mask: boardMask + rotMask pieces: pArray ifFound: solutionBlock)]].
        perms swap: i with: j].
    ^count

ShootoutMeteorBoard default solvedPuzzleDo: [:e | ]

-> 3096168

Since 10 loops are necessary for producing a solution, (but 11 are counted) and since there are 2098 solutions, that's 23078 productive loops...

First, we perform a minor refactoring: we handle permutations with a color bit mask, beacause it is more in line with the style of our algorithm:

searchPuzzlesWithColorMask: colorMask boardMask: boardMask pieces: pArray ifFound: solutionBlock
    | boardLowBit colorBit count |
    count := 1.
    colorMask = 0
        ifTrue:
            [solutionBlock value: (self boardStringWithPieces: pArray).
            ^count].
    boardLowBit := (boardMask + 1) lowBit.
    colorBit := 1.
    1 to: pArray size do: [:colorRank |
        (colorMask bitAnd: colorBit) = 0
            ifFalse: [
                ((positionsPerPiece at: colorRank) at: boardLowBit) do: [:rotMask |
                    0 = (boardMask bitAnd: rotMask)
                        ifTrue:
                            [pArray at: colorRank put: rotMask.
                            count := count + (self searchPuzzlesWithColorMask: colorMask - colorBit boardMask: boardMask + rotMask pieces: pArray ifFound: solutionBlock)]]].
        colorBit := colorBit * 2].
    ^count

solvedPuzzleDo: solutionBlock
^self searchPuzzlesWithColorMask: 1 << pieces size - 1 boardMask: 0 pieces: pieces copy ifFound: solutionBlock

Now let's look an individual piece position:

| board |

board printMask: (((board := ShootoutMeteorBoard default) instVarNamed: 'pieces') at: 5)

...
0 0 0 0 0
0 0 1 0 0
0 0 0 1 1
0 0 1 0 1

Ah yes, some are always invalid, because they generate a hole that we will never be able to fill. How could we detect this ?

For south row, this is easy, if there is one or more zero holes between ones, then there is a hole we cannot fill. Translated into bit tricks:

hasInsetZero: aMask
    | allOnes |
    allOnes := aMask bitOr: aMask - 1.
    ^(allOnes bitAnd: allOnes + 1) > 0

The first operation replaces trailing zeroes with ones:

  110011000 mask          0111000

| 110010111 mask - 1    | 0110111

  110011111               0111111

The second one replace trailing ones with zeroes

  110011111 allOnes       0111111 

& 110100000 allOnes + 1 & 1000000

110000000 0000000

If the result is different from zero, then yes, there is an inset zero...

hasSouthOrNorthIsland: aMask
^(self hasInsetZero: (southEdge bitAnd: aMask))

or: [self hasInsetZero: (northEdge bitAnd: aMask)]

What about east and west edges ?

If we bitAnd the piece mask with eastEdge, then we are looking for some bit pattern like

00001

00000

00001

How to remove the four west most columns? Well we don't, we fill them with 1s multiplying with 11111 (this happens to be our southEdge)

11111

00000

11111

Magically, we can use the hasInsetZero: trick on above bit pattern

hasEastOrWestIsland: aMask
^ (self hasInsetZero: southEdge * (eastEdge bitAnd: aMask))

or: [self hasInsetZero: southEdge * (westEdge bitAnd: aMask)]

There is one more case to eliminate, which is the hole in a corner. But a piece could fit in that corner if it has exactly 5 cell holes. In fact, a single piece can generate up to 7 holes:

...
0 0 0 0 0
0 0 0 1 1
0 0 1 0 0
0 0 1 0 0
0 1 0 0 0

One solution is to fill the mask with ones and count the filled bits

fill: aMask startingAt: pos count: countBlock
    | filled |
    (aMask bitAnd: pos) = 0 ifFalse: [^aMask].
    countBlock value.
    filled := aMask + pos.
    (self canShiftE: pos)
        ifTrue: [filled := self fill: filled startingAt: (self shiftE: pos) count: countBlock].
    (self canShiftNE: pos)
        ifTrue: [filled := self fill: filled startingAt: (self shiftNE: pos) count: countBlock].
    (self canShiftNW: pos)
        ifTrue: [filled := self fill: filled startingAt: (self shiftNW: pos) count: countBlock].
    (self canShiftW: pos)
        ifTrue: [filled := self fill: filled startingAt: (self shiftW: pos) count: countBlock].
    (self canShiftSW: pos)
        ifTrue: [filled := self fill: filled startingAt: (self shiftSW: pos) count: countBlock].
    (self canShiftSE: pos)
       ifTrue: [filled := self fill: filled startingAt: (self shiftSE: pos) count: countBlock].
    ^filled

It could be interesting to factor above code with a loop on directions as proposed in Part1, but we'll let that as an exercise. We could also change the countBlock into a fillCount instance variable, we are state-full after all. We won't bother right now.

fill: aMask startingAt: pos
    | count |
    count := 0.
    self fill: aMask startingAt: pos count: [count := count + 1].
    ^count

hasCornerIsland: aMask edge: verticalEdge edge: horizontalEdge
    | corner |
    ^(aMask bitAnd: verticalEdge) > 0
        and: [(aMask bitAnd: horizontalEdge) > 0
        and: [(aMask bitAnd: (corner := verticalEdge bitAnd: horizontalEdge)) = 0
        and: [(self fill: aMask startingAt: corner) \\ 5 > 0]]]

Let's detect these simple cases altogether and remove them from possible positions:

hasIsland1: aMask
    ^(self hasEastOrWestIsland: aMask)
    or: [(self hasCornerIsland: aMask edge: southEdge edge: eastEdge)
    or: [(self hasCornerIsland: aMask edge: southEdge edge: westEdge)
    or: [(self hasCornerIsland: aMask edge: northEdge edge: eastEdge)
    or: [(self hasCornerIsland: aMask edge: northEdge edge: westEdge)
    or: [self hasSouthOrNorthIsland: aMask]]]]]

Again, some tests will be duplicated, but we don't care yet, this is not a big CPU contributor.

initializePossiblePositions
    positionsPerPiece := pieces collect: [:aPiece |
        | possible |
        possible := (Array new: ncell) collect: [:lowBit | Set new: 12].
        self rotationsOf: aPiece do: [:rotated |
            self placesFor: rotated do: [:shifted |
                (self hasIsland1: shifted) ifFalse: [(possible at: shifted lowBit) add: shifted]]].
        possible collect: [:e | e asArray]].

MessageTally spyOn: [ShootoutMeteorBoard solveDefault].

Still gives the good result, and now take a bit less than 13s.

ShootoutMeteorBoard default solvedPuzzleDo: [:e | ]
-> 1999211
OK, we saved about 33% loops, and gained about 33% speed with this simple trick.

Now we can also see that turning the board by 180° gives us another possible solution. Great! we can add the pieces into a Set and stop as soon as we reach the quota of solutions.

Turning 180° is just reversing bit pattern, or reversing the compact string solution:

ShootoutMeteorBoard class>>solveDefaultUntilFound: n
    | board solutions |
    solutions := Set new: n.
    board := ShootoutMeteorBoard default.
    ^[:print |
        board solvedPuzzleDo:
            [:aString |
            (solutions add: aString; add: aString reversed; size) >= n ifTrue: [^print value] ].
        print value] value:
            [solutions := solutions sorted.
            String streamContents:
                [:outputStream |
                outputStream print: solutions size; nextPutAll: ' solutions found'; cr; cr.
                board printSolution: solutions first on: outputStream.
                outputStream cr.
                board printSolution: solutions last on: outputStream]]

Deceivingly, we don't get a factor 2 speed up, only 2 seconds...
We transform the loopCount into an instance variable and count loops again:

-> 1664525 loops

Hmm case of unlucky ordering of pieces?

This is at http://ss3.gemstone.com/ss/Shootout/Shootout.blog-nice.3.mcz

But we don't really need a Set... We can omit the emission of those symmetric solutions by removing 3 rotations (and flipped rotations) of a single piece:
(ShootoutMeteorBoard default instVarNamed: 'positionsPerPiece')
collect: [:e | e detectSum: [:e2 | e2 size]]
-> #(195 221 231 228 219 273 231 263 183 228)

The piece of rank 6 seems more interesting

rotationsOf: aPieceMask do: aBlock
    | next |
    aBlock value: (next := aPieceMask); value: (self flip: next).
    ((aPieceMask = (pieces at: 6)) ifTrue: [2] ifFalse: [5])
        timesRepeat: [aBlock value: (next := self rotate: next); value: (self flip: next)]

Then we add the symmetric solution again like in our Set variant:

ShootoutMeteorBoard class>>solveDefault
    ^String streamContents: [:outputStream |
        | board count minSolution maxSolution |
        count := 0.
        minSolution := String new: 50 withAll: $9.
        maxSolution := String new: 50 withAll: $0.
        (board := ShootoutMeteorBoard default) solvedPuzzleDo:
            [:direct |
                {direct. direct reversed} do: [:aString |
                    count := count + 1.
                    aString < minSolution ifTrue: [minSolution := aString].
                    aString > maxSolution ifTrue: [maxSolution := aString]]. ].
        outputStream print: board loopCount; nextPutAll: ' loops'; cr;cr.
        outputStream print: count; nextPutAll: ' solutions found'; cr; cr.
        board printSolution: minSolution on: outputStream.
        outputStream cr.
        board printSolution: maxSolution on: outputStream]

MessageTally spyOn: [ShootoutMeteorBoard solveDefault].

-> 1081569 loops

And a bit less than 7s.

This is at http://ss3.gemstone.com/ss/Shootout/Shootout.blog-nice.4.mcz

We have already saved 66% computations with easy steps.

In next iteration, we'll try harder.

Meteor contest Part 2 - the naive solution

To continue our Meteor contest, we now have to enumerate all possible pavements...

First, a method was missing in previous post:

Integer>>bitsDo: aBlock
    | mask |
    self < 0 ifTrue: [^self error: 'Cannot enumerate bits of a negative integer'].
    mask := self.
    [mask = 0]
       whileFalse:
           [aBlock value: mask lowBit.
           mask := mask bitAnd: mask - 1]

The trick here is that mask bitAnd: mask - 1 sets the lowest 1-bit to 0, thus having as many loops as the bitCount. The rest is obvious.

  10110011000 mask 

& 10110010111 mask - 1

10110010000

For the case of LargeInteger, it is more efficient to avoid repeated allocation of LargeInteger:

LargePositiveInteger>>bitsDo: aBlock
     | mask offset |
   1 to: self digitLength do: [:iByte |
       offset := iByte - 1 << 3.
       mask := (self digitAt: iByte).
       [mask = 0]
           whileFalse:
               [aBlock value: mask lowBit + offset.
               mask := mask bitAnd: mask - 1]]

One possibility is to generate all permutations of colors like squeak does:

searchPuzzlesWithPermutation: perms rank: i
    | index boardLowBit |
    i > perms size ifTrue: [ ^perms ].
    i to: perms size do: [:j |
       perms swap: i with: j.
       self searchPuzzlesWithPermutation: perms rank: i + 1]].
       perms swap: i with: j].
    ^nil

Then, we need to pass the occupied boardMask and isolate the next free cell in it, which happens to be the lowest bit at 0. Using a bit trick
nextFreeCell := (boardMask + 1) lowBit.

Or we could have used (boardMask bitXor: boardMask + 1) highBit
Indeed:

  10110010111 boardMask 

^ 10110011000 boardMask + 1

00000001111

We then have to enumerate possible orientations of i-th color piece which fit into this board,

((positionsPerPiece at: colorIndex) at: nextFreeCell) do: [:rotMask |

If the place is free

0 = (boardMask bitAnd: rotMask)

then record the rotated and shifted pieceMask into pArray (indexed by color rank)

pArray at: index put: rotMask.

and continue with permutations of higher ranks

self searchPuzzlesWithPermutation: perms rank: i + 1 mask: boardMask + rotMask pieces: pArray ifFound: solutionBlock.

When found, the string-i-fied solution is passed as argument to a block of code.
i > perms size ifTrue: [ ^solutionBlock value: (self boardStringWithPieces: pArray) ].

Putting it all together:

searchPuzzlesWithPermutation: perms rank: i mask: boardMask pieces: pArray ifFound: solutionBlock
    | colorIndex boardLowBit |
    i > perms size ifTrue: [ ^solutionBlock value: (self boardStringWithPieces: pArray) ].
   nextFreeCell := (boardMask + 1) lowBit.
    i to: perms size do: [:j |
        colorIndex := perms at: j.
        perms swap: i with: j.
        ((positionsPerPiece at: colorIndex) at: nextFreeCell) do: [:rotMask |
            0 = (boardMask bitAnd: rotMask)
                ifTrue:
                    [pArray at: index put: rotMask.
                    self searchPuzzlesWithPermutation: perms rank: i + 1 mask: boardMask + rotMask pieces: pArray ifFound: solutionBlock]].
        perms swap: i with: j].
    ^nil

At upper level, we initiate the search loop with:

solvedPuzzleDo: solutionBlock
self searchPuzzlesWithPermutation: (1 to: pieces size) asArray rank: 1 mask: 0 pieces: pieces copy ifFound: solutionBlock

We then have to print the solutions, in the compact form:

or in the nicer board form:

printSolution: aString on: aStream
    | src |
    src := aString readStream.
    [src atEnd]
        whileFalse:
            [1 to: ncol do: [:j | aStream nextPut: src next; space].
            aStream cr.
            1 to: ncol do: [:j | aStream space; nextPut: src next].
            aStream cr]

and we are done:

ShootoutMeteorBoard class>>solveDefault
   ^String streamContents: [:outputStream |
       | board count minSolution maxSolution |
       count := 0.
       minSolution := String new: 50 withAll: $9.
       maxSolution := String new: 50 withAll: $0.
       (board := ShootoutMeteorBoard default) solvedPuzzleDo:
           [:aString |
               count := count + 1.
               aString < minSolution ifTrue: [minSolution := aString].
               aString > maxSolution ifTrue: [maxSolution := aString]. ].
       outputStream print: count; nextPutAll: ' solutions found'; cr; cr.
       board printSolution: minSolution on: outputStream.
       outputStream cr.
       board printSolution: maxSolution on: outputStream]

This works great, but rather slooowwwwly.
MessageTally spyOn: [ShootoutMeteorBoard solveDefault].
2098 solutions found

0 0 0 0 1
2 2 2 0 1
2 6 6 1 1
2 6 1 5 5
8 6 5 5 5
8 6 3 3 3
4 8 8 9 3
4 4 8 9 3
4 7 4 7 9
7 7 7 9 9

9 9 9 9 8
9 6 6 8 5
6 6 8 8 5
6 8 2 5 5
7 7 7 2 5
7 4 7 2 0
1 4 2 2 0
1 4 4 0 3
1 4 0 0 3
1 1 3 3 3

Almost 20 seconds on my mac mini with a Cog VM...
The version is at http://ss3.gemstone.com/ss/Shootout/Shootout.blog-nice.1.mcz

In next steps we are going to break this work

Playing with Meteor contest in Smalltalk - Part 1

The shootout language benchmarks can be a source of fun too.

It will provide us this little exercise: the meteor contest.

The game consists in finding all the possible pavements of a 10x5 hexagonal cells board with 10 pieces of 5 hexagons.

The obvious idea is to represent the board as a bit mask, where 1 means occupied cells, 0 means free cells.

Also each piece can be represented as a bit mask of occupied cell.

I decided to initialize the pieces from the string representation of a board:

ShootoutMeteorBoard class>>default
^self basicNew fromString:
'0 0 0 0 1
2 2 2 0 1
2 6 6 1 1
2 6 1 5 5
8 6 5 5 5
8 6 3 3 3
4 8 8 9 3
4 4 8 9 3
4 7 4 7 9
7 7 7 9 9'

Where each digit represent a different color.

The least significant bit is located at SE corner.

The pieces can be rotated in 6 different directions (E, NE, NW, W, SW, SE) and also flipped, giving 12 different orientations.

Here is our trivial ideas:

Generate all possible positions of a piece on the board.
Then fill the board starting at south east, and progressing along the bit rank, toward W, then N directions.
Try this for all different permutations of the colors, so as to generate all the possible solutions.

This is a rather naive algorithm, which won't be efficient because it will explore too many combinations. But let's do it, we will then see how to reduce the costs.

Lets see how to shift toward different directions : a piece cannot be shifted if it is already on an edge, so we have to test that. Let's just do it with bit masks. We modify a bit fromString:
ncell := rawString size.

ncol := (aString readStream upTo: Character cr) count: [:e | e isSeparator not].
nrow := ncell / ncol.

twoRows := ncol * 2.
self initializeRowColMasks.

And

ShootoutMeteorBoard>>initializeRowColMasks
    southEdge := 1 << ncol - 1.
    southToNorthMasks := (0 to: nrow - 1) collect: [:i | southEdge << (ncol * i)].
    northEdge := southToNorthMasks last.
    eastEdge := 1<<ncell-1/southEdge.
    eastToWestMasks := (0 to: ncol - 1) collect: [:i | eastEdge << i].
    westEdge := eastToWestMasks last.
    oddRowsMask := 1<<ncell-1/(1<<twoRows-1)*southEdge.
    evenRowsMask := oddRowsMask << ncol.

Ah, good old bit tricks, remember the division I used in variant of Sieve of Eratosthenes?

If not, you can read old pages of this blog.

You see, 2r000010000100001 * 2r11111 -> 2r111111111111111, so just reverse the operation and you'll see that the division is the a way to create a bit pattern with one bit set every n bit.

Just for debug purpose, printing a mask won't be luxury

We just have to enumerate all bits from HSB to LSB.

ShootoutMeteorBoard>>printMask: aPieceMask on: aStream
| mask |
mask := 1 << ncell.

    nrow to: 1 by: -1 do: [:irow |
        irow odd ifTrue: [aStream space].
        1 to: ncol do: [:icol |
            aStream
                nextPut: ((aPieceMask bitAnd: (mask := mask >> 1)) = 0
                    ifTrue: [$0]
                    ifFalse: [$1]);
                space].
        aStream cr]

ShootoutMeteorBoard>>printMask: aPieceMask
^String new: ncol * 2 + 1 * nrow streamContents: [:aStream | self printMask: aPieceMask on: aStream]

Lets try, ShootoutMeteorBoard default printMask: (ShootoutMeteorBoard default instVarNamed: 'oddRowsMask')

Good, in Smalltalk indices starts at 1, the the 1st row is odd (this is the south one).

But why do we need masks for odd and even rows ? Lets see with color 7

...
0 0 0 0 0
0 7 0 7 0
7 7 7 0 0

There is one cell on the W edge, but since it is on an odd row this does not prevent us to shift toward NW.

...
7 0 7 0 0
7 7 7 0 0
0 0 0 0 0

Now there is a 7 on an even west edge, we cannot shift toward NE anymore...

Indeed, NE east shift shifts odd rows by 5, but even rows are shifted by 1 more bit. Thus:

shiftNW: aPieceMask
    | evens odds |
    odds := oddRowsMask bitAnd: aPieceMask.
    evens := evenRowsMask bitAnd: aPieceMask.
    ^evens << 1 + odds << ncol

shiftNE: aPieceMask
    | evens odds |
    odds := oddRowsMask bitAnd: aPieceMask.
    evens := evenRowsMask bitAnd: aPieceMask.
    ^odds >> 1 + evens << ncol

shiftSE: aPieceMask
    | evens odds |
    odds := oddRowsMask bitAnd: aPieceMask.
    evens := evenRowsMask bitAnd: aPieceMask.
    ^odds >> 1 + evens >> ncol

shiftSW: aPieceMask
    | evens odds |
    odds := oddRowsMask bitAnd: aPieceMask.
    evens := evenRowsMask bitAnd: aPieceMask.
    ^evens << 1 + odds >> ncol

Shifting toward E, or W is more trivial:

shiftE: aPieceMask
^aPieceMask >> 1

shiftW: aPieceMask
^aPieceMask << 1

But we must first test if we can shift.

canShiftNW: aPieceMask

^(aPieceMask bitAnd: (northEdge bitOr: (westEdge bitAnd: evenRowsMask))) = 0

canShiftNE: aPieceMask

^(aPieceMask bitAnd: (northEdge bitOr: (eastEdge bitAnd: oddRowsMask))) = 0

canShiftSW: aPieceMask

^(aPieceMask bitAnd: (southEdge bitOr: (westEdge bitAnd: evenRowsMask))) = 0

canShiftSE: aPieceMask

^(aPieceMask bitAnd: (southEdge bitOr: (eastEdge bitAnd: oddRowsMask))) = 0

canShiftE: aPieceMask
^(aPieceMask bitAnd: eastEdge) = 0

canShiftW: aPieceMask
^(aPieceMask bitAnd: westEdge) = 0

If we want to enumerate the directions, we can encode them with an integer:

ShootoutMeteorBoard>>initializeShiftMasks
    E := 1.
    NE := 2.
    NW := 3.
    W := 4.
    SW := 5.
    SE := 6.

(canShiftMasks := Array new: 6)

at: E put: eastEdge;

at: NE put: (northEdge bitOr: (eastEdge bitAnd: oddRowsMask)));

at: NW put: (northEdge bitOr: (westEdge bitAnd: evenRowsMask)));

at: W put: westEdge;

at: SW put: (southEdge bitOr: (westEdge bitAnd: evenRowsMask)));

at: SE put: (southEdge bitOr: (eastEdge bitAnd: oddRowsMask))

And now, let's test:

canShift: aPieceMask toward: aDirection

^(aPieceMask bitAnd: (canShiftMasks at: aDirection)) = 0

And if we want to generalize the shifts too,

(oddShifts := Array new: 6)

at: E put: -1;

at: NE put: ncol - 1;

at: NW put: ncol;

at: W put: 1;

at: SW put: 0 - ncol;

at: SE put: -1 - ncol.

(evenShifts := Array new: 6)

at: E put: -1;

at: NE put: ncol;

at: NW put: ncol + 1;

at: W put: 1;

at: SW put: 1 - ncol;

at: SE put: 0 - ncol.

shift: aPieceMask toward: aDirection

^((aPieceMask bitAnd: oddRowsMask) bitShift: (oddShifts at: aDirection)) +

((aPieceMask bitAnd: evenRowsMask) bitShift: (evenShifts at: aDirection))

Now we are going to shift a piece to the south east most position possible.

shiftSEmost: aPieceMask
    | mostSEMask eastColumn lowBit |
    aPieceMask odd ifTrue: [^aPieceMask].
    lowBit := aPieceMask lowBit.
    mostSEMask := aPieceMask >> (lowBit - 1 // (2*ncol) * (2*ncol)).
    (mostSEMask bitAnd: southEdge) = 0
        ifTrue: [mostSEMask := (self canShiftSE: mostSEMask)
            ifTrue: [self shiftSE: mostSEMask]
            ifFalse: [self shiftSW: mostSEMask]].
    eastColumn := eastToWestMasks findFirst: [:e | (e bitAnd: mostSEMask) > 0].
    ^mostSEMask >> (eastColumn - 1)

Does it work ?

| board |
(board := ShootoutMeteorBoard default) printMask: (board shiftSEmost: ((board instVarNamed: 'pieces') at: 8))

...
0 0 0 0 0
0 0 0 0 0
0 0 1 0 1
0 1 1 1 0
OK, fine.

Now, let's flip the piece. we can flip east-west or north-south.

But in both cases, the odd rows becomes even, so we have to shift by ncol in either direction.

flip: aPieceMask
   ^self shiftSEmost: (southToNorthMasks
        inject: 0 into: [:mask :rowMask |
            mask << ncol + ((rowMask bitAnd: aPieceMask) >> (rowMask lowBit - 1))]) >> ncol

flipEW: aPieceMask
   ^self shiftSEmost: (eastToWestMasks
        inject: 0 into: [:mask :columnMask |
            mask << 1 + ((columnMask bitAnd: aPieceMask) >> (columnMask lowBit - 1))]) << ncol

That's it... Here is the N-S version

| board |
(board := ShootoutMeteorBoard default) printMask: (board flip: ((board instVarNamed: 'pieces') at: 8))

...

and the E-W

| board |
(board := ShootoutMeteorBoard default) printMask: (board flipEW: ((board instVarNamed: 'pieces') at: 8))

Now the rotation is a bit more involving with just these shifts...

Here is our algorithm, assuming the piece is shifted SE most:

1) choose the pivot point, invariant during rotation on south row, west most.

rotate: aPieceMask
    | rotatedMask pivot rotatedPivot irow row |
    rotatedMask := 0.
    irow := 1.
    row := aPieceMask bitAnd: (southToNorthMasks at: irow).
    rotatedPivot := pivot := 1 << (row highBit - 1).

2) choose to rotate anti clock wise and sweep the board toward east. if the next cell is toward east, then the next rotated cell will be toward NE.

    [rotatedMask := rotatedMask + rotatedPivot.
    [(row bitAnd: pivot - 1) = 0]
        whileFalse:
            [pivot := self shiftE: pivot.
            rotatedPivot := self shiftNE: rotatedPivot.
            (row bitAnd: pivot) = 0
                ifFalse:
                    [rotatedMask := rotatedMask + rotatedPivot]].

3) Now continue on next row, toward north, until it is empty

(row := aPieceMask bitAnd: (southToNorthMasks at: (irow := irow + 1))) = 0]
whileFalse:

4) Begin scanning NE of pivot if possible
            [(self canShiftNE: pivot)
                ifTrue:
                    [pivot := self shiftNE: pivot.
                    rotatedPivot := self shiftNW: rotatedPivot]
                ifFalse:
                    [pivot := self shiftNW: pivot.
                    rotatedPivot := self shiftW: rotatedPivot].

5) Then scan toward west, until west most is encountered on this row
            [row >= (pivot << 1)]
                whileTrue:
                    [pivot := self shiftW: pivot.

Notice that shifting SW is dangerous, we might rotate a cell out of the board.

We detect such case and move the rotated piece two rows up.

                    (self canShiftSW: rotatedPivot)
                        ifFalse:
                            [rotatedPivot := rotatedPivot << twoRows.
                            rotatedMask := rotatedMask << twoRows.].
                    rotatedPivot := self shiftSW: rotatedPivot]].

6) Now we're done, just shift the result SE most to follow our convention
^self shiftSEmost: rotatedMask

Surely a bit too clever, this will by far be our worse method. Let's just try it:

| board |
(board := ShootoutMeteorBoard default) printMask: (board rotate: ((board instVarNamed: 'pieces') at: 8))

Good. If I tell you that above code was not my first attempt, you'll probably believe me. How great it is to have an interpreted language and a good debugger, it's invaluable for those super fast feedback loops.

Now lets enumerate possible orientations of a piece:

Then shift a piece at every possible position.
We start SE most, then scan toward W while possible, then toward NE if possible, else toward NW, this gonna give us a quite regular code:

placesFor: aPieceMask do: aBlock
    | westMask eastMask cellNumber |
    eastMask := self shiftSEmost: aPieceMask.

    [[westMask := eastMask.
    [aBlock value: westMask.
    self canShiftW: westMask] whileTrue: [westMask := self shiftW: westMask].
    self canShiftNE: eastMask] whileTrue: [eastMask := self shiftNE: eastMask].
    self canShiftNW: eastMask] whileTrue: [eastMask := self shiftNW: eastMask]

We are ready to fulfil our first algorithm step, that is to generate all possible pieces positions, which is the map of possible position for each orientation of each piece

self placesFor: rotated do: [:shifted |

(possible at: shifted lowBit) add: shifted]].
possible collect: [:e | e asArray]].

Our step is finished, let's put it together

    twoRows := ncol * 2.
    nrow := ncell / ncol.
    self initializeRowColMasks.
    pieces := rawString asSet asSortedArray collect: [:char |
        self shiftSEmost:
            (rawString inject: 0 into: [:pmask :c | pmask * 2 + (c = char ifTrue: [1] ifFalse: [0])])].
    self initializePossiblePositions.

Next step will be in next post. If you are impatient, code is located at squeaksource3. But it might be better with explanations, I was rather sparing with comments, following the existing Shootout code style.

Friday, January 27, 2012

Which fraction has these digits ?

Sven Van Caekenberghe posted these nice example on Pharo mailing list:

1/9801 asScaledDecimal: 200.
1/998001 asScaledDecimal: 3000.

This also works with:

1/9801 printShowingDecimalPlaces: 200.

We could conjecture that 1/99980001 asScaledDecimal: 40000 should give a nice serie of digits too.

But if we want to form our own fraction, how can we construct it ?

Let us consider a fraction in interval ]0,1[

In "decimal" notation (I say decimal, but it is valid in any other base b), this fraction print in base b with a head composed of head size digits, and a repeated tail composed of tail size digits 0.head_tail_tail_tail...

For example, b := 10, head := '123' and tail := '40', will lead to (fraction asScaledDecimal: tail size * 3 + head size) = '0.123404040s9'

I once thought we could extend the syntax with such notation to denote infinite repetition:

fraction := 0.123(40).

But we didn't, so what is the value of the fraction?

Let us work with numbers rather than strings:

h := head ifEmpty: [0] ifNotEmpty: [Number readFrom: head base: b].

t := tail ifEmpty: [0] ifNotEmpty: [Number readFrom: tail base: b].

nh := head size.

nt := tail size.

If we had infinite series (lazy one would be preferred), we could define:

positiveIntegers := 1 to: Integer infinity.

infiniteTail := (positiveIntegers collect: [:i | t / (b raisedTo: nt * i)]).

fraction := (infiniteTail + h) / (b raisedTo: nh).

We have no such ready made infinity, nor lazy collections...

But even with such objects, we cannot evaluate this infiniteTail that easily.

Really? If we shift it by tail size digits, we will obtain this number with same infiniteTail:

infiniteTail * (b raisedTo: nt) = (t + infiniteTail).

From this property, it is easy to get:

infiniteTail := t / ((b raisedTo: nt) - 1).

If we want a leading head, we simply have to shift above fraction:

fraction := (t / ((b raisedTo: nt) - 1) + h) / (b raisedTo: nh).

So, let us reconstruct Sven example:

| base tail |
base := 10.
tail := ((0 to: 97) , #(99) inject: '' writeStream into: [:s :n | n printOn: s base: base length: 2 padded: true. s]) contents.
^(Number readFrom: tail base: base) / ((base raisedTo: tail size) - 1)

-> 1/9801

Good. Now, this expression won't lead to nice fraction for any tail...

If we want a fraction with numerator = 1, we need to find a divisor of the denominator 999....99. That's the nice thing with Sven example. I don't know how or where he found it, but I guess it'll be hard to be as elegant. Oh, maybe thanks to the extension to any base you can play with it too (you will need to implement printShowingDecimalPlaces:base:).