\ http://projecteuler.net/index.php?section=problems&id=213

\ Problem:

\ A 30x30 grid of squares contains 900 fleas, initially one flea per
\ square.  When a bell is rung, each flea jumps to an adjacent square
\ at random (usually 4 possibilities, except for fleas on the edge of
\ the grid or at the corners).

\ What is the expected number of unoccupied squares after 50 rings of
\ the bell? Give your answer rounded to six decimal places.

: f+! ( r addr -- )
    dup f@ f+ f! ;

30 constant xlen
30 constant ylen

0 value nextstate
0 value currentstate

: genstate ( r -- addr )
    \ generate a grid with all squares containing r
    xlen ylen * floats dup allocate throw ( u addr )
    swap 0 ?do ( addr )
        fdup dup i + f!
    1 floats +loop
    fdrop ;

: up ( r addr -- )
    \ adds r to the field upwards from addr
    xlen floats - f+! ;

: down ( r addr -- )
    \ adds r to the field downwards from addr
    xlen floats + f+! ;

: left ( r addr -- )
    \ adds r to the field leftwards from addr
    1 floats - f+! ;

: right ( r addr -- )
    \ adds r to the field rightwards from addr
    1 floats + f+! ;

: middle ( r addr -- )
    \ r is the number of fleas in a square, addr identifies the next
    \ state of this square.  Updates the adjacent squares.
    0.25e f*
    fdup dup up
    fdup dup down
    fdup dup left
    fdup dup right
    fdrop drop ;

: top-edge ( r addr -- )
    \ r is the number of fleas in a square, addr identifies the next
    \ state of this square.  Updates the adjacent squares.
    3e f/
    fdup dup down
    fdup dup left
    fdup dup right
    fdrop drop ;

: bottom-edge ( r addr -- )
    \ r is the number of fleas in a square, addr identifies the next
    \ state of this square.  Updates the adjacent squares.
    3e f/
    fdup dup up
    fdup dup left
    fdup dup right
    fdrop drop ;

: left-edge ( r addr -- )
    \ r is the number of fleas in a square, addr identifies the next
    \ state of this square.  Updates the adjacent squares.
    3e f/
    fdup dup up
    fdup dup down
    fdup dup right
    fdrop drop ;

: right-edge ( r addr -- )
    \ r is the number of fleas in a square, addr identifies the next
    \ state of this square.  Updates the adjacent squares.
    3e f/
    fdup dup up
    fdup dup down
    fdup dup left
    fdrop drop ;

: tl-corner ( r addr -- )
    \ r is the number of fleas in a square, addr identifies the next
    \ state of this square.  Updates the adjacent squares.
    0.5e f*
    fdup dup down
    fdup dup right
    fdrop drop ;

: tr-corner ( r addr -- )
    \ r is the number of fleas in a square, addr identifies the next
    \ state of this square.  Updates the adjacent squares.
    0.5e f*
    fdup dup down
    fdup dup left
    fdrop drop ;

: bl-corner ( r addr -- )
    \ r is the number of fleas in a square, addr identifies the next
    \ state of this square.  Updates the adjacent squares.
    0.5e f*
    fdup dup up
    fdup dup right
    fdrop drop ;

: br-corner ( r addr -- )
    \ r is the number of fleas in a square, addr identifies the next
    \ state of this square.  Updates the adjacent squares.
    0.5e f*
    fdup dup up
    fdup dup left
    fdrop drop ;

: p ( addr1 -- addr2 r addr3 )
    \ addr1 identifies a square in currentstate, r is its contents,
    \ addr2 is the next square, and addr3 is the square in nextstate
    \ corresponding to addr1.
    dup f@ dup float+ swap currentstate - nextstate + ;

: ring ( -- )
    \ ring the bell (see problem statement)
    0e genstate to nextstate
    currentstate ( addr )
    p tl-corner xlen 2 ?do p top-edge loop p tr-corner
    ylen 2 ?do
        p left-edge xlen 2 ?do p middle loop p right-edge
    loop
    p bl-corner xlen 2 ?do p bottom-edge loop p br-corner
    drop
    currentstate nextstate to currentstate free throw ;

: combine-probs ( acc dist -- acc )
    \ combine the probability distribution of a flea with the others up to now
    xlen ylen * floats 0 ?do
        1e dup i + f@ f- over i + dup f@ f* f!
    1 floats +loop
    drop ;

: flea ( u -- )
    \ compute the probability distribution of flea u in currentstate
    0e genstate to currentstate
    1e currentstate swap floats + f!
    50 0 ?do ring loop ;

: flea-circus ( -- r )
    \ r is the number of expected empty fields.  Compute the
    \ probabilities for all squares by combining the probabilities of
    \ all fleas, then sum up these probabilities.
    1e genstate
    xlen ylen * 0 ?do
        i flea currentstate combine-probs
        currentstate free throw
    loop
    0e xlen ylen * 0 ?do
        dup f@ f+ float+
    loop ;

flea-circus f. cr