Random Numbers Between 0 and n-1

Phil Pilgrim (PhiPi)

I was doing some AVR programming and need a quick, short random number generator that produced values between zero and a limit - 1. The usual way to do this is to generate a full-sized random number, divide it by the limit, and use the remainder. But this was too complicated, and I didn't have room for a division routine anyway. So I started to think about how LFSRs are written, and it dawned on me that each step in the LFSR could produce one bit of a number that I could use to multiply the limit by to produce a fractional product (i.e. the product's upper bits). From there it became just a matter of inserting a conditional add and a shift into the LFSR routine.

So I wrote a version for the Propeller, based on Chip's LFSR in the Spin byte code interpreter. It required two additional instructions:

[b]CON[/b]

  [b]_clkmode[/b]      = [b]xtal1[/b] + [b]pll16x[/b]
  [b]_xinfreq[/b]      = 5_000_000

[b]VAR[/b]

  [b]long[/b]          arg[noparse][[/noparse]*2]

[b]OBJ[/b]

  sio   :       "FullDuplexSerial"

[b]PUB[/b] Start

  sio.start(31, 30, 0, 9600)
  [b]cognew[/b](@random, @arg)
  [b]repeat[/b] 1000
    sio.dec(limit_rand(100_000))
    sio.tx(13)

[b]PUB[/b] limit_rand(limit)

'Return a uniformly-distributed random number between 0 and limit - 1.

  arg := limit
  [b]repeat[/b] [b]while[/b] arg
  [b]return[/b] arg[noparse][[/noparse]*1]

[b]DAT[/b]

              [b]org[/b]       0
random        [b]mov[/b]       xaddr,[b]par[/b]               'Get address of input value.
              [b]mov[/b]       yaddr,[b]par[/b]               'Get address of output value.
              [b]add[/b]       yaddr,#4

:waitarg      [b]rdlong[/b]    x,xaddr [b]wz[/b]              'Is input value non-zero?
        [b]if_z[/b]  [b]jmp[/b]       #:waitarg               '  No:  Try again.

              [b]mov[/b]       y,#0                    'Initialize result.
              [b]mov[/b]       bitcnt,#32              'Initialize bit count.
                                                                  
:randlp       [b]test[/b]      rand,#%10111 [b]wc[/b]         'LFSR: Get parity of tap values.
              [b]rcr[/b]       rand,#1                 'Shift odd parity in though top.
        [b]if_c[/b]  [b]add[/b]       y,x [b]wc[/b]                  'Add input to product if carry (like multiply). ***
              [b]rcr[/b]       y,#1                    'Shift carry into result from top.              ***
              [b]djnz[/b]      bitcnt,#:randlp         'Back for another bit.

              [b]wrlong[/b]    y,yaddr                 'Save result in hub.
              [b]mov[/b]       x,#0                    'Zero input.
              [b]wrlong[/b]    x,xaddr                 'Save zero in hub to signal done.
              [b]jmp[/b]       #:waitarg               'Dum de dum...             

rand          [b]long[/b]      $1234_5678              'Random seed.

xaddr         [b]res[/b]       1                       'Hub address of input.
yaddr         [b]res[/b]       1                       'Hub address of result.
x             [b]res[/b]       1                       'Input value.
y             [b]res[/b]       1                       'Result value.
bitcnt        [b]res[/b]       1                       'Bit counter.
             
                                                '*** These two statements added to regular LFSR for limiting.
  

To give it a quick test, I ran the above code and plotted each output value against the one following it in a scatter chart to see if there were any obvious correlations. There don't seem to be any, but more rigorous tests would not be inappropriate:



This method suffers the same shortcoming that the remainder method does, though: it maps a domain of 2³² - 1 values into a range of n values. If n doesn't evenly divide the domain size, the relative frequency of some result values will be higher than others by one.

-Phil

mpark

Another PhiPi stroke of genius!

ekosminsky

PhiPi,
I really like the code you created. I'm trying to create a game that requires an array filled with random integers, so your limit_rand method works perfectly. However, i was trying to make it so that the seed value changed so the game wouldn't have the same array of "random" integers by setting the seed value: rand long, to CNT, but this doesn't seem to be working. Here is part of my code:

PUB simonGame

  Listen.start(31, 30, 0, 9600)
  waitcnt(cnt + 1440000000)
  cognew(@random, @arg) 
  Listen.dec(cnt)
  repeat 30
     Listen.dec(limitRand(10))
     Listen.tx(13)

Listen is an object of FullDuplexSerial. I also have the line:

rand long cnt

in the DAT block.

Suggestions?

Phil Pilgrim (PhiPi)

You would have to do this explicitly, either in Spin (rand := cnt) before the cognew or in assembly (mov rand,cnt). I'm not sure this would give you the desired result, though, since cnt may always have the same value at startup. You may want to use Chip's RealRandom object (in the OBEX) to produce a seed — or as a complete replacement for my code.

-Phil

ekosminsky

You would think that CNT would have the same value at startup and when it creates the array, but since it increments so quickly, there are often one or two digits inconsistent every run, so the value is much more random. Thank you.

Also, I'm fairly new to spin. What is the OBEX and how can i implement an object of the RealRandom class?

Phil Pilgrim (PhiPi)

OBEX == Propeller Object Exchange.

See the manual under OBJ for general info and Chip's RealRandom object for implementation specifics.

-Phil

northcove

Phil, just curious...for your application why wouldn't the Spin ? operator suffice?

Peter Nachtwey

I am sure the method above will work for most people but for the purist:
http://www.codeproject.com/KB/recipes/SimpleRNG.aspx
0 to (2^32)-1 can just map into 0 to (1-1/(2^32))

Phil Pilgrim (PhiPi)

northcove,

Yes, it would, so long as speed was not a factor. The thrust of my effort was providing the quickest possible 0..n-1 pseudorandom generator.

-Phil

northcove

Ah yes. I would be surprised if the built-in implementation was performed in fewer cycles than yours but as you know those guys are pretty clever.