Thursday, May 8, 2014

Problem 84 - Red: Monopoly, bugs in a beta, RNG issues, oh my

Continuing in my trend that I began with Zimbu, I decided to use another "fledgling" language for problem 84. This language is Red. Red is based on the REBOL language, that I used earlier in this challenge. Therefore it has a weird love for square brackets, but otherwise a very nice, readable, easy to learn syntax. However, the language is currently at version 0.4.2. This caused...some issues. Notably, many things are in the language spec that have not been implemented. Some examples of issues I ran into:

Horrible bug with array indexing. If you have a list L, you can access the first element with L/1. Indeed, you can even do fancy stuff like L/(3 + 4 - 6). However, you cannot do L/x, where x is any variable. The solution? L/(x + 0) does work, though it took me quite a while to find.

Switch statements are nice, and the syntax for them is quite simple in Red. However, though switch statements have been implemented in the language, switch statements with a default clause have not been implemented...great...this explains my switch cases that do nothing basically for a lot of inputs.

I decided that I would use Red for this problem, which involves running simulations of Monopoly. In order to simulate Monopoly, one must roll dice. In order to roll dice, one must generate random numbers. I started with the bare knowledge I have of Lehmer's algorithm / Linear Congruential generators for forming random numbers (Using a sequence z_n = a^n + b mod m) for some parameters a,b, m. However, everywhere I looked online gave good parameters that would have required integers larger than 32 bits (or at least unsigned 32 bit integers), which Red does not support. Then, I thought I would implement the suggestion of a math professor, which was to take the ones digit of cos(n^2) for some n, and this generates a pretty random sequence. However, this would be just about impossible: not only does Red not have a library cosine function, Red does not yet support floating point arithmetic. So, I went back to the drawing board, until I randomly guessed LCG parameters that gave the right answer. This eventually led to me coming across the right answer by sheer luck. Of course, I felt bad about that, so I went back to the drawing board to get more legitimately random numbers. These were my conditions for my random numbers:

1) Consistency across number of trials: the answer should converge relatively uniformly, I shouldn't get the right results after 1000 trials, and then have very wrong results after 100000, that would be a bad sign.

2) Consistency across input problems: When I first found a bad LCG that gave the right answer to the euler problem, it did not give the correct answer in the class 6-sided dice monopoly case, so it was clearly illegitimate.

3) Consistency in theory: I did not allow myself to guess parameters. I went with large primes for m, 0 for b (thus I used Lehmer's algorithm, not a general LCG), and for a, I chose the primitive root of m, as this is supposed to be a good choice.

Faced with these conditions, I tried to look for more help from the internet, maybe finding some magical small LCG. But no, I found a solution in a somewhat obvious observation in a text file somewhere on the internet (ftp://ftp.ietf.org/rfc/rfc1750.txt): if you have two kind of random, uncorrelated sequences, and you xor them together, you get a theoretically more random sequence. So, I took three large primes with different primitive roots (3, 5, 17), and I started them at different points in the sequence (to decrease any chance for correlation), and I xored the results of the three LCGs together. These gave random numbers that led to the correct results in a pretty consistent manner. I will admit that this is probably not the world's best sequence of random numbers, but it did pretty well, and I am fairly happy with the theory behind it: for example, adding another "correlated" LCG (one with the same a), never helped my results, which supports the theory behind all this. Anyway, after all of these difficulties that led to many insights, my solution is below. Runs in 14s on my machine (I am doing 200000 trials). Also on the time: Red is apparently supposed to  be able to be compiled, but every Red program that I have written that can be legally interpreted is treated as an "Invalid Red Program" by the compiler, so I gave up. Code is here:
lastrand: 25
lastrand2: 17
lastrand3: 59049
rand: func [ 
    return: [integer!] 
] [
    lastrand: (lastrand * 5) // 163847
    lastrand2: (lastrand2 * 17) // 8191
    lastrand3: (lastrand3 * 3) // 104743
    lastrand xor lastrand2 xor lastrand3
]
 
squares: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]

pos: 0
turn: 0
ccn: 0
chn: 0
while [turn < 200000] [
    goagain: true
    dubcount: 0
    while [goagain] [
        d1: 1 + ((rand) % 4)
        d2: 1 + ((rand) % 4)
        either d1 = d2 [
            goagain: true
            dubcount: dubcount + 1
        ] [
            goagain: false
        ]
        either dubcount = 3 [
            pos: 11
            goagain: false
        ] [
            pos: ((pos + d1 + d2) % 40)
            switch pos [
                0 [pos: 40]
                31 [pos: 11
                    goagain: false]
                3 18 34 [
                    switch ccn [
                        0 [pos: 1]
                        1 [pos: 11
                           goagain: false]
                        2 3 4 5 6 7 8 9 10 11 12 13 14 15 [pos: pos]
                    ]
                    ccn: (ccn + 1) % 16
                ]
                8 23 37 [
                    switch chn [
                        0 [pos: 1]
                        1 [pos: 11
                           goagain: false]
                        2 [pos: 12]
                        3 [pos: 25]
                        4 [pos: 40]
                        5 [pos: 6]
                        6 7 [switch pos [
                                    8 [pos: 16]
                                    23 [pos: 26]
                                    37 [pos: 6]
                            ]]
                        8 [switch pos [
                                    8 37 [pos: 13]
                                    23 [pos: 29]
                            ]]
                        9 [pos: pos - 3]
                        10 11 12 13 14 15 [pos: pos]
                    ]
                    chn: (chn + 1) % 16
                ]
                1 2 4 5 6 7 9 10 11 12 13 14 15 16 17 19 20 [pos: pos]
                21 22 24 25 26 27 28 29 30 32 33 34 35 36 38 39 40 [pos: pos]
            ]
        ]
    ]
    squares/(pos + 0): squares/(pos + 0) + 1
    turn: turn + 1  
]
caps: [1000000 0 0 0]
ans: [0 0 0]
i: 1
while [i < 4][
    j: 1
    while [j < 41][
        if all [squares/(j + 0) < caps/(i + 0) squares/(j + 0) > caps/(i + 1)][
            caps/(i + 1): squares/(j + 0)
            ans/(i + 0): j - 1
        ]
        j: j + 1
    ]
    i: i + 1
]
print squares
print ans

No comments:

Post a Comment