Jump to content
IGNORED

IntyBASIC Special Random Selection?


First Spear

Recommended Posts

Hi all. I have an array with 30 elements, and I want to efficiently select one item from the array that is only in a certain range. Just 1-6,7,12,13,18,19,24,25-30. These elements represent the edge of a box.

 

What is the right pattern to think about to do this, more importantly what is an efficient way to get one of those numbers (or pull at random and exclude the inner numbers)?

 

Thanks.

 

 

post-31694-0-01226800-1549561022.jpg

Link to comment
Share on other sites

I don't know if it is the most efficient way, but the way I would attack this problem would be to construct an array with the border cells and just pick an index in that range at random.

RectBorder:
  Data  1,  2,  3,  4,  5
  Data  6, 12, 18, 24
  Data 30, 29, 28, 27, 26
  Data 25, 19, 13, 7
 

-dZ.

  • Like 1
Link to comment
Share on other sites

I was about to offer pretty much the same suggestion dZ-Jay just gave: Use a lookup table to map your sparse number space to a dense array.

 

You have 18 numbers to choose from, which is slightly inconvenient, since it's not a power of 2. Still, if you do something like this...

.

X = RAND(18)

.

...to pick one of the 18 values, you'll get a nearly uniformly distributed result.

 

4 of the values will get chosen slightly more often than the others. The elements at indexes 0, 4, 9, and 13 will get chosen 15/256 times, while the rest will get chosen 14/256 times. It's a small bias, but it's there.

 

I'm not sure it's worth worrying about.

  • Like 2
Link to comment
Share on other sites

I don't understand the usage here but my thought was the same. a lookup table using just the number that make up the edges.

 

however, The mention of the bias is eye opening. where does the bias come from? I made a (partial) Asteroids type game last year and the color of the asteroids was random but some colors seemed to appear more frequently than others and I couldn't find a fix for it.

Link to comment
Share on other sites

I believe internally random numbers are generated in the range 0-255, i.e. 256 different possibilities. If you ask for a range 0-15, it equals 256/16 = each 16 of the original range is scaled down to your desired range.

 

If you ask for a range 0-17 like in this case, it equals 256/18 = each 14.22222 of the original range is scaled down. The fraction will carry over as an error until it has exceeded 1, which explains the slight bias.

 

Different ranges would distribute the error differently, though uniformly. Usually you'd use a=RANDOM(6)+2 to get a number 2-7. I'm not entirely sure if this (slower) version makes the numbers more evenly distributed:

 

DO : a=RANDOM( 8 ) : LOOP UNTIL a>1

Edited by carlsson
Link to comment
Share on other sites

I believe internally random numbers are generated in the range 0-255, i.e. 256 different possibilities. If you ask for a range 0-15, it equals 256/16 = each 16 of the original range is scaled down to your desired range.

 

If you ask for a range 0-17 like in this case, it equals 256/18 = each 14.22222 of the original range is scaled down. The fraction will carry over as an error until it has exceeded 1, which explains the slight bias.

 

Different ranges would distribute the error differently, though uniformly. Usually you'd use a=RANDOM(6)+2 to get a number 2-7. I'm not entirely sure if this (slower) version makes the numbers more evenly distributed:

 

DO : a=RANDOM( 8 ) : LOOP UNTIL a>1

 

Generating a power-of-two random number, and then looping until it fits within a non-power-of-two range, does ensure that the random number will be uniform provided the random number generator is uniform over its power-of-2 range. However, it does have a non-deterministic run time. The further the range is from a power of two, the more times it will loop. 6 out of 8 isn't too bad—each iteration has a 3/4 chance of being in range, and 1/4 chance of looping. 18/32, on the other hand has a 7/16 chance of looping—nearly 1/2.

 

Another strategy might be to "rotate the bias:"

.

  X = RANDOM(18) + XROT
  IF X > 18 THEN X = X - 18
  IF XROT = 0 THEN XROT = 17 ELSE XROT = XROT - 1

.

The bias is still there, but it will be spread across all values equally across multiple random numbers, as it rotates the number space by 1 each time it's called.

 

This might be slightly faster, and doesn't use an extra variable; however, it only rotates the bias once per frame:

.

X = (18 * ((RANDOM(256) + FRAME) AND 255)) / 256

.

As carlsson says, the bias arises from the fact that RAND and RANDOM start with an 8-bit random number, and then multiply by the random number range. For values that divide evenly into 256 (e.g. any power of 2 from 1 to 256), this naturally gives a uniform result. For other values, you'll get a few values over-represented.

 

I wrote a short C program to determine which values would be over-represented by RAND(18) and RANDOM(18):

.

#include <stdio.h>
  
int h[18];

int main() {
    for (int i = 0; i < 256; i++) {
        int j = (i * 18) / 256;
        h[j]++;
    }

    for (int i = 0; i < 18; i++) {
        printf("%2d: %d\n", i, h[i]);
    }
    return 0;
}

.

It just counts how many times each of the values 0 through 17 would occur by stepping through all of the possible random number generator values 0 through 255. Most of them occur 14 times (as 256/18 = 14.222...). Four of them occur 15 times.

 

If you replace that 18 with 16, you will see all 16 values occur exactly 16 times.

 

Another approach would be to up the math to 15 bits, and reduce the bias to approximately 1821-vs-1820:

.

#TEMP = RANDOM(128) * 256 + RANDOM(256)
X = (#TEMP + (#TEMP / ) / 8 / 256

.

There are a lot of strategies for dealing with the bias. The smaller the random number range relative to the underlying random number, the less the bias. (Again, assuming the underlying random number generator is uniform over its native range. That's 0..255 for RAND and RANDOM; 0..32767 for my last example above.)

Edited by intvnut
  • Like 2
Link to comment
Share on other sites

It is possible I am introducing my own errors here, but I ran this program:

 

 

 MODE 1
  DIM #counter(32)
 
  FOR i=0 TO 31:#counter(i)=0:NEXT
 
  FOR #j=0 TO 1000
    k=RAND(32)
    #counter(k)=#counter(k)+1
    WAIT
  NEXT
 
  FOR i=0 TO 31
    PRINT AT i*5 COLOR $0007,<3>#counter(i)
  NEXT
 
end:
  GOTO end
and got these results:

post-5454-0-61962100-1549623277.gif

 

I replaced the call with RANDOM(32) which doesn't require WAIT which allowed me to extend the loop to 10000 iterations. It gave me these results:

post-5454-0-05062800-1549623359.gif

 

Unless I have myself to blame, the distribution of pseudo-random numbers doesn't seem perfectly uniform.

Link to comment
Share on other sites

OK, so I wrote a C model of IntyBASIC's RANDOM function, and I have to say I'm not very impressed with it. If you hold FRAME constant, it only cycles between 30 values. If you do cycle through all 65536 values of FRAME, it spreads out a bit better but isn't quite uniform, even after many many thousands of random numbers.

 

One option is to skip IntyBASIC's RANDOM and implement your own LFSR-based random number generator, burning a 16-bit variable to do it.

 

Unfortunately, my favorite source for LFSR taps is no longer online: http://www.newwaveinstruments.com/resources/articles/m_sequence_linear_feedback_shift_register_lfsr.htm

 

And archive.org is currently down for maintenance. :-P When it comes back up, you can find nice tables of LFSR taps there.

 

In any case, I can use a known-good polynomial, such as 0xAD52. With that polynomial, you can implement the following LFSR based RAND, if you don't mind blowing a 16-bit variable for a larger seed:

.

NextRand: PROCEDURE
          #SEED = (((#SEED = 0) XOR FRAME) XOR RandTbl(#SEED AND 15)) XOR (#SEED / 16)
          #SEED = RandTbl(#SEED AND 15) XOR (#SEED / 16)
END
RandTbl:  DATA $0000, $4303, $8606, $C505
          DATA $56A9, $15AA, $D0AF, $93AC 
          DATA $AD52, $EE51, $2B54, $6857
          DATA $FBFB, $B8F8, $7DFD, $3EFE

DEF FN MyRand(x) (((x) * ($FF AND #SEED)) / 256)

.

That should spread more uniformly across 0..255, with 0 being slightly under-represented. (255 vs. 256 occurrences across the 65535 element sequence.)

 

Note: I think I found an IntyBASIC bug, or I just don't understand how DEF FN is supposed to work. Paging nanochess! Take a look at this:

.

' Input to IntyBASIC:
DEF FN MyRand(val) (((val) * ($FF AND #SEED)) / 256)
X = MyRand(18)
X = (18 * ($FF AND #SEED)) / 256

' Output from IntyBASIC 1.4.0.  Notice the call to MyRand just zeros X:
    ;[21] X = MyRand(18)
    SRCFILE "/tmp/r.bas",21
    CLRR R0
    MVO R0,var_X
    ;[22] X = (18 * ($FF AND #SEED)) / 256
    SRCFILE "/tmp/r.bas",22
    MVI var_&SEED,R0
    ANDI #255,R0
    MULT R0,R4,18
    SWAP R0
    MVO R0,var_X

.

That looks like a bug to me. Either that, or as I said, I don't actually know how DEF FN works. The call to MyRand(18) just zeros X rather than expanding to the equivalent expression below and actually computing a value.

Edited by intvnut
  • Like 2
Link to comment
Share on other sites

OK, Archive.org came back up. I was able to capture all of their maximal-length LFSR tap tables, and make a PDF of the website. I'll share it here as well.

 

And this goes to the archive.org capture: https://web.archive.org/web/20160313124937/http://www.newwaveinstruments.com/resources/articles/m_sequence_linear_feedback_shift_register_lfsr.htm

 

EDIT: It appears my attachment is bogus, as the various "stages" files just contain HTML complaining about me scraping the site. I'll have to redo my work tonight. :-P

 

EDIT 2: I've downloaded all of the tap files properly now, and updated the ZIP file. If you had downloaded it previously, re-download with this new, larger, more useful file.

LFSR_Taps.zip

Edited by intvnut
  • Like 1
Link to comment
Share on other sites

OK, so I wrote a C model of IntyBASIC's RANDOM function, and I have to say I'm not very impressed with it. If you hold FRAME constant, it only cycles between 30 values. If you do cycle through all 65536 values of FRAME, it spreads out a bit better but isn't quite uniform, even after many many thousands of random numbers.

 

One option is to skip IntyBASIC's RANDOM and implement your own LFSR-based random number generator, burning a 16-bit variable to do it.

 

Unfortunately, my favorite source for LFSR taps is no longer online: http://www.newwaveinstruments.com/resources/articles/m_sequence_linear_feedback_shift_register_lfsr.htm

 

And archive.org is currently down for maintenance. :-P When it comes back up, you can find nice tables of LFSR taps there.

 

In any case, I can use a known-good polynomial, such as 0xAD52. With that polynomial, you can implement the following LFSR based RAND, if you don't mind blowing a 16-bit variable for a larger seed:

.

NextRand: PROCEDURE
          #SEED = (((#SEED = 0) XOR FRAME) XOR RandTbl(#SEED AND 15)) XOR (#SEED / 16)
          #SEED = RandTbl(#SEED AND 15) XOR (#RAND / 16)
END
RandTbl:  DATA $0000, $4303, $8606, $C505
          DATA $56A9, $15AA, $D0AF, $93AC 
          DATA $AD52, $EE51, $2B54, $6857
          DATA $FBFB, $B8F8, $7DFD, $3EFE

DEF FN MyRand(x) (((x) * ($FF AND #SEED)) / 256)

.

That should spread more uniformly across 0..255, with 0 being slightly under-represented. (255 vs. 256 occurrences across the 65535 element sequence.)

 

Note: I think I found an IntyBASIC bug, or I just don't understand how DEF FN is supposed to work. Paging nanochess! Take a look at this:

.

' Input to IntyBASIC:
DEF FN MyRand(val) (((val) * ($FF AND #SEED)) / 256)
X = MyRand(18)
X = (18 * ($FF AND #SEED)) / 256

' Output from IntyBASIC 1.4.0.  Notice the call to MyRand just zeros X:
    ;[21] X = MyRand(18)
    SRCFILE "/tmp/r.bas",21
    CLRR R0
    MVO R0,var_X
    ;[22] X = (18 * ($FF AND #SEED)) / 256
    SRCFILE "/tmp/r.bas",22
    MVI var_&SEED,R0
    ANDI #255,R0
    MULT R0,R4,18
    SWAP R0
    MVO R0,var_X

.

That looks like a bug to me. Either that, or as I said, I don't actually know how DEF FN works. The call to MyRand(18) just zeros X rather than expanding to the equivalent expression below and actually computing a value.

 

I need to check this but DEF FN should include a "equal" sign like:

 

 

DEF FN MyRand(val) = (((val) * ($FF AND #SEED)) / 256)
Link to comment
Share on other sites

 

 

I need to check this but DEF FN should include a "equal" sign like:

DEF FN MyRand(val) = (((val) * ($FF AND #SEED)) / 256)

 

Ah, that did the trick. I didn't notice the missing '=' because I'm so used to #define in C.

 

Why wasn't it a compile error to omit the equals sign, though?

Link to comment
Share on other sites

That's what I'm going to check :dunce:

 

Aha... it did give me some errors:

.

Error: missing = in DEF FN in line 19
Warning: invalid extra characters in line 19
Error: bad syntax for expression in line 20

.

I just didn't notice them because, well, it still generated output. :D

 

I guess I'm spoiled by compilers that refuse to give me output if there's a compile error. It forces me to take notice. ;-)

  • Like 1
Link to comment
Share on other sites

I tried to implement intvnut's code above and got extremely strange results, including that letters F, H, J were printed as part of the routine that prints a three digit number!! First I thought there were gremlins in the compiler, but then I spotted a typo in the NextRand procedure where one occurrence of #SEED was mistyped as the uninitialized variable #RAND. I know that can be avoided by setting Option Explicit, so I have to blame myself for not doing that as a habit.

After fixing that, I'm getting results as per the screenshot above. Possibly I've still made some mistake in how to use the routine, or otherwise that particular LFSR doesn't seem any better than the built-in PRNG, at least not in this implementation/context.

 DEF FN MyRand(x) = (((x) * ($FF AND #SEED)) / 256)
 
  MODE 1
 
  DIM #counter(32)
 
  FOR i=0 TO 31:#counter(i)=0:NEXT
 
  #SEED = 0
 
  FOR #j=0 TO 10000
    k=MyRand(32)
    GOSUB NextRand
    #counter(k)=#counter(k)+1
  NEXT
 
  FOR i=0 TO 31
    PRINT AT i*5 COLOR $0007,<3>#counter(i)
  NEXT
 
end:
  GOTO end
 
NextRand: PROCEDURE
  #SEED = (((#SEED = 0) XOR FRAME) XOR RandTbl(#SEED AND 15)) XOR (#SEED / 16)
  #SEED = RandTbl(#SEED AND 15) XOR (#SEED / 16)
  END
 
RandTbl:  DATA $0000, $4303, $8606, $C505
          DATA $56A9, $15AA, $D0AF, $93AC 
          DATA $AD52, $EE51, $2B54, $6857
          DATA $FBFB, $B8F8, $7DFD, $3EFE


post-5454-0-07675600-1549658731.gif

Link to comment
Share on other sites

I tried to implement intvnut's code above and got extremely strange results, including that letters F, H, J were printed as part of the routine that prints a three digit number!! First I thought there were gremlins in the compiler, but then I spotted a typo in the NextRand procedure where one occurrence of #SEED was mistyped as the uninitialized variable #RAND.

 

Oops, I'll go back and edit my post.

 

 

After fixing that, I'm getting results as per the screenshot above. Possibly I've still made some mistake in how to use the routine, or otherwise that particular LFSR doesn't seem any better than the built-in PRNG, at least not in this implementation/context.

post-5454-0-07675600-1549658731.gif

 

The standard deviation for a binomial distribution is given by sqrt(n * p * (1 - p)).

 

You can view each bin independently, with a 1/32 chance of being chosen, and 31/32 chance of not being chosen. Thus, for 32 bins and 10000 trials, p = 1/32 and n = 10000. The standard deviation is therefore approximately sqrt(10000 * (1/32) * (1 - 1/32)) = 17.4 (approx).

 

The distribution of outcomes is expected to be normal. Therefore, you should expect 68% of the samples to fall within +/- 17.4 of 312.5 (that is, 295.1 to 329.9), and ~95% of the samples to fall within +/- 34.8 of 312.5 (that is, 277.7 to 347.3).

 

That spread looks like it's a little wider than suggested by the binomial distribution. I wonder if it gets better if you remove "XOR FRAME"?

Link to comment
Share on other sites

Without XOR FRAME, I get this distribution:

 

post-5454-0-22280000-1549660308.gif

I put those into a Google Sheet and calculated the population standard deviation, and it came up 16.29. That sounds pretty reasonable—maybe even a shade low—compared to the theoretical 17.4.

Edited by intvnut
Link to comment
Share on other sites

I put the numbers into an Excel sheet. The results for RAND(32) w/ WAIT are skewed due to I only ran that loop for 1/10 of the other tests and multiplied each bin by 10.

attachicon.gifrndtest.xls

 

Now that I'm back at my desk, I was able to load this into gnumeric and take a look. (Previously, I was stuck in a meeting and was limited in what I could do.)

 

The RAND(32) w/ WAIT comes up as 51.85. The theoretical stdev would be 10 * SQRT(1000 * (1/32) * (1 - 1/32)) = 55.02. So, RAND(32) w/ WAIT looks like it's in the right ballpark.

 

Overall, AD52 w/out FRAME seems like the best performer. AD52 vs. RANDOM(32) look similar in performance, at least on this test.

 

In the context of stdev, the small bias introduced by RANDOM(18) is small compared to the standard deviation, and you're unlikely to notice it in practice.

  • Like 1
Link to comment
Share on other sites

On the topic of pseudo-random number generators:

I regularly participate in the yearly Comp.Sys.Sinclair Crap Game Compo. While I usually don't "win" that compo neither, I like it because it is more about making clumsy games with sometimes obscure mechanics based on popular or funny topics, and noone expects you to optimize and polish your code, quite the opposite as you get questioned if you spend too much time developing your games. Also each edition of the CSSCGC tends to have a couple of special challenges, make games on certain subjects or within certain limitations.

One of the challenges in the 2017 edition was to make a Sinclair BASIC game without using any keywords achieved in the E cursor mode. For those of you who don't know about the ZX Spectrum, you are not able to type in BASIC keywords by hand on the 16K and 48K models (though you can on the later Spectrum 128, +2 and +3 models). Instead each key has up to 5 keywords/symbols assigned to it, which you obtain by various combinations of the two Caps Shift and Symbol Shift keys. Each time you press one or both shift keys, the cursor changes to a different letter to tell you which mode you are in.

If you are prohibited from reaching the E cursor mode, the following keywords are no longer available:
ABS, ACS, ASN, ATN, ATTR, BEEP, BIN, BRIGHT, CAT, CHR$, CIRCLE, CLOSE#, CODE, COS, DATA, DEF, FN, ERASE, EXP, FLASH, FN, FORMAT, IN, INK, INKEY$, INT, INVERSE, LEN, LINE, LLIST, LN, LPRINT, MERGE, MOVE, OPEN#, OUT, OVER, PAPER, PEEK, PI, POINT, READ, RESTORE, RND, SCREEN$, SGN, SIN, SQR, STR$, TAB, TAN, USR, VAL$, VAL, VERIFY plus the symbols [ ] ~ | \ { }, copyright and embedded colour control codes.

However you still have access to these keywords:
AND, AT, BORDER, CLEAR, CLS, CONT, COPY, DIM, DRAW, FOR, GOSUB, GOTO, IF, INPUT, LET, LIST, LOAD, NEW, NEXT, NOT, OR, PAUSE, PLOT, POKE, PRINT, RANDOMIZE, REM, RETURN, RUN, SAVE, STEP, STOP, THEN, TO plus the symbols ! " # $ % & ' ( ) * + , . / : ; ? @ < <= <> <= > ^_ £

(rather curious to have access to RANDOMIZE but not be able to use RND to get a value from it...)



Anyway, I decided to make a text based golf game in the good old style of the player typing in which club, strength and angle of the stroke and the game will plot that as ASCII graphics with PRINT AT statements. I had a couple of challenges to tackle:

* Simulate trigonometry without the built-in SIN and COS, which I solved with the help from the BASIC Handbook by David A. Lien
* Simulate the INT function which I did with FOR loops that makes the game extra slow
* Simulate the RND function so the game wouldn't be entirely deterministic

I tried to implement a LFSR but due to Sinclair BASIC only has logical AND/OR instead of bitwise (and no XOR or shift functions at all), I found myself having to re-implement bitwise operators from scratch. While I actually managed to get it to work, it ran super dog slow on the already slow BASIC, almost like I could run the LFSR algorithm by paper and pencil and be done before the computer. That is when I started to look online for other types of PRNG and found a site which describes a routine supposedly present in the C64 version of Defender of the Crown! This routine has a data set of 15 initial values, and then for each random value picked modifies the array. I believe it has an entropy (?) of 1 million or something like that. The best part you only need addition and multiplication, no bitwise operations at all. Of course it uses a bit of RAM but at least when it came to the crap game that was the least of my worries. Perhaps it is less useful on the Intellivision where RAM tends to be expensive.

I have implemented this PRNG as below, and ran the same test. I found that the standard deviation depends a lot on the starting seed (from 12.71 to 21.19), but perhaps that is as expected. There may be more PRNG's of this kind, or with different starting sets of values.

  MODE 1
  DIM #counter(32)
  DIM prng(15)
  SIGNED ry,rx
 
  FOR i=0 TO 31:#counter(i)=0:NEXT
 
  REM Initialize DOTC PRNG
  FOR i=0 TO 14 : prng(i)=prngdef(i) : NEXT
 
  REM Set starting seed
  rx=5:ry=rx-1:k=0 : REM Edit: Or hmm.. perhaps it should be ry=rx+1 to fit better with the algorithm
  
  FOR #j=0 TO 10000
    GOSUB getrand
    k = k AND 31 : REM scale down result to desired range, done with floating point division & multiplication on the ZX Spectrum
    #counter(k)=#counter(k)+1
  NEXT
 
  FOR i=0 TO 31
    PRINT AT i*5 COLOR $0007,<3>#counter(i)
  NEXT
 
end:
  GOTO end
 
getrand: PROCEDURE
  rx=rx-1:if rx<0 then rx=14
  k=prng(ry)+prng(rx):if k>255 then k=k-256
  prng(ry)=k:ry=rx
  END
 
prngdef:
  DATA 241,174,84,90,228,165,57,157,19,96,163,28,71,173,107

rndtest.xls

Edited by carlsson
  • Like 3
Link to comment
Share on other sites

I trimmed the quote below to make it a bit more compact.

 

I tried to implement a LFSR but due to Sinclair BASIC only has logical AND/OR instead of bitwise (and no XOR or shift functions at all), I found myself having to re-implement bitwise operators from scratch. While I actually managed to get it to work, it ran super dog slow on the already slow BASIC, almost like I could run the LFSR algorithm by paper and pencil and be done before the computer. That is when I started to look online for other types of PRNG and found a site which describes a routine supposedly present in the C64 version of Defender of the Crown! This routine has a data set of 15 initial values, and then for each random value picked modifies the array. I believe it has an entropy (?) of 1 million or something like that. The best part you only need addition and multiplication, no bitwise operations at all. Of course it uses a bit of RAM but at least when it came to the crap game that was the least of my worries. Perhaps it is less useful on the Intellivision where RAM tends to be expensive.

I have implemented this PRNG as below, and ran the same test. I found that the standard deviation depends a lot on the starting seed (from 12.71 to 21.19), but perhaps that is as expected. There may be more PRNG's of this kind, or with different starting sets of values.

rx=5:ry=rx-1:k=0 : REM Edit: Or hmm.. perhaps it should be ry=rx+1 to fit better with the algorithm

'... 

getrand: PROCEDURE
  rx=rx-1:if rx<0 then rx=14
  k=prng(ry)+prng(rx):if k>255 then k=k-256
  prng(ry)=k:ry=rx
  END

 

This is a form of lagged Fibonacci random number generator. They can be quite good if you pick the distance between rx and ry and the total array size correctly. The LSB ends up being an LFSR, and the upper bits behave as LFSRs perturbed by the carries (or borrows) of the bits below it. (Knuth covers this in TAOCP Volume 2. There's an online preview here. Click on 3.2.2 and scroll down a couple pages.)

 

I need to check, but I think the lag 1 might be OK, based on the fact that [15,14] provides a maximal length 15-stage LFSR.

I've used lagged Fibonacci generators for years because they're exceptionally cheap computationally. They do need a fair bit of RAM though. If I'm not mistaken, this generator should have a period of 28 * 65,535 - 1= 16,776,959. I need to test it to be sure.

 

One downside of lagged Fibonacci is that it has a lot of inertia, and that may explain your stdev results. In particular, it'll produce pretty bad results for a long time if you initialize it with [1, 0, 0, .... 0]. Also, you need to ensure there's at least one non-zero bit in the LSB of one of the bytes in the initializer, otherwise that bit will remain 0 forever.

 

Other generators, such as WELL, do a better job of "escaping zeroland." It's considerably more expensive though, and needs even more RAM.

 

I did write an implementation of WELL512a in CP-1610 assembly back in 2013; however it's based on a license restricted C implementation, so I haven't made it public. I did get permission to use it in my own Intellivision games, but so far I haven't gotten around to using it.

Edited by intvnut
Link to comment
Share on other sites

I've used lagged Fibonacci generators for years because they're exceptionally cheap computationally. They do need a fair bit of RAM though. If I'm not mistaken, this generator should have a period of 28 * 65,535 - 1= 16,776,959. I need to test it to be sure.

 

OK, I wrote a quick C program to measure the DotC RNG, using an initial state of [1, 0, 0, ... ]. It apparently has a period of 62,912,640, which is larger than I expected.

 

Here's my C code if anyone wants to check it.

#include <assert.h>
#include <stdbool.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

static uint8_t state[15] = { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
static uint8_t orig [15] = { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
static uint8_t rx = 5, ry = 6;

static uint8_t get_rand(void) {
    uint8_t k;
    rx = rx ? rx - 1: 14;
    k = state[rx] + state[ry];
    state[ry] = k;
    ry = rx;
    return k;
}

static bool is_orig(void) {
    return rx == 5 && !memcmp(state, orig, sizeof(state));
}

int main() {
    long long period = 0;
    assert(is_orig());

    do {
        period++;
        get_rand();
    } while (!is_orig());

    printf("Period = %lld\n", period);

    return 0;
}

.

.

 

EDIT: I just changed to use the DotC initializer, and the period seems to be the same, once I fixed a bug in my cycle-detector. 62.9M.

Edited by intvnut
Link to comment
Share on other sites

OK, Archive.org came back up. I was able to capture all of their maximal-length LFSR tap tables, and make a PDF of the website. I'll share it here as well.

 

And this goes to the archive.org capture: https://web.archive.org/web/20160313124937/http://www.newwaveinstruments.com/resources/articles/m_sequence_linear_feedback_shift_register_lfsr.htm

 

EDIT: It appears my attachment is bogus, as the various "stages" files just contain HTML complaining about me scraping the site. I'll have to redo my work tonight. :-P

 

EDIT 2: I've downloaded all of the tap files properly now, and updated the ZIP file. If you had downloaded it previously, re-download with this new, larger, more useful file.

 

Here's an LFSR implementation from my framework. It's actually based on code you, Arnauld, and I worked on in the INTVPROG list several years ago. I guess it could be added to IntyBASIC as an Assembly Language module and called with a macro. I posted it in the forum a few years ago.

 

All the ancillary modules are macro libraries which are part of my framework, and support validation of parameters and error checking; and the routine includes a parameter for deciding whether to unroll the loop for speed.

 

However, none of that is needed really. The core routine is below, although I still recommend that if you use it, you use the fully integrated set of modules, for they include simple, single-line callable macros such as for ranged values, or to seed the PRNG, that can then be invoked from IntyBASIC using a DEF'FN macro.

;* ======================================================================== *;
;*  These routines are placed into the public domain by their author.  All  *;
;*  copyright rights are hereby relinquished on the routines and data in    *;
;*  this file.  -- James Pujals (DZ-Jay), 2014                              *;
;* ======================================================================== *;

PRNG_POLY       EQU     $0000B301                       ; Maximal length sequence (M-Sequence) polynomial:
                                                        ;       x^32 = x^15 + x^13 + x^12 + x^9 + x^8 + x^0
_rnd            STRUCT  0
@@poly_hi       EQU     ((PRNG_POLY SHR 16) AND $FFFF)
@@poly_lo       EQU     ( PRNG_POLY         AND $FFFF)

                ; Need to declare the following variables
@@seed_hi       SYSTEM  1
@@seed_lo       SYSTEM  1
                ENDS

;; ======================================================================== ;;
;;  PM.PRNG.RAND:                                                           ;;
;;  Generates a pseudo-random number of a given magnitude.  The PRNG is a   ;;
;;  Galois implementation of a 32-bit Linear Feedback Shift Register (LFSR).;;
;;                                                                          ;;
;;  The polynomial used is configured via the PRNG_POLY symbol.             ;;
;;                                                                          ;;
;;  There are two entry points to this procedure:                           ;;
;;      PM.PRNG.RAND                 Receives arguments as input record.    ;;
;;                                                                          ;;
;;      PM.PRNG.RAND.1               Receives number of bits in register.   ;;
;;                                                                          ;;
;;  INPUT for PM.PRNG.RAND                                                  ;;
;;      R5      Pointer to invocation record, followed by return address.   ;;
;;              Number of random bits to generate.      1 DECLE             ;;
;;              Mask corresponding to the bits.         1 DECLE             ;;
;;                                                                          ;;
;;  INPUT for PM.PRNG.RAND.1                                                ;;
;;      R0      Number of random bits to generate.                          ;;
;;      R5      Pointer to invocation record, followed by return address.   ;;
;;              Mask corresponding to the bits.         1 DECLE             ;;
;;                                                                          ;;
;;  OUTPUT                                                                  ;;
;;      R0      Pseudo-random bits.                                         ;;
;;      R1      Trashed (high-order word of polynomial).                    ;;
;;      R2      Trashed.                                                    ;;
;;      R3      Trashed (when OPTIMIZED_SP is defined).                     ;;
;; ======================================================================== ;;
PM.PRNG.RAND    PROC
                MVI@    R5,     R0                      ; Get number of bits
@@1:            MOVR    R0,     R2                      ; Make a copy of the number of bits

                MVI     _rnd.seed_lo,           R0      ; \_ Get the current LFSR values
                MVI     _rnd.seed_hi,           R1      ; /

_r.i            QSET    0

            IF (_rnd.poly_hi = 0)
                MVII    #_rnd.poly_lo,          R3      ;
            ENDI

                SLL     R2,     2                       ; \
                SUBI    #@@__update,            R2      ;  |_ Determine entry point into unrolled loop
                NEGR    R2                              ;  |    entry = (exit - (mask * size)) : size = 4
                MOVR    R2,     PC                      ; /

            REPEAT (16)

                ; Shift the LFSR one bit at a time,
                ; for how many bits were requested.
                ; --------------------------------------
                DSLLC   R1, R0, 1

              IF (_rnd.poly_hi = 0)

                ADCR    PC                              ; if (!carry) ...
                XORR    R3,     R0                      ;   lfsr.lo ^= POLY.lo; // (low-order word)

              ELSE

                BC      @@__no_xor[_r.i]                  ; if (!carry) ...
                XORI    #_rnd.poly_lo,          R0      ;   lfsr.lo ^= POLY.lo; // (low-order word)
                XORI    #_rnd.poly_hi,          R1      ;   lfsr.hi ^= POLY.hi; // (high-order word)
@@__no_xor[_r.i]:

              ENDI

_r.i            QSET    (_r.i + 1)

            ENDR

@@__update:     MVO     R0,     _rnd.seed_lo            ; \_ Update the LFSR
                MVO     R1,     _rnd.seed_hi            ; /
                AND@    R5,     R0                      ; Mask the output bits

                JR      R5
                ENDP

The code above unrolls its loop for maximum speed. However, the original module in the ZIP file allows you to select this with an option.

 

It's free to use by anybody who cares or who wants to examine how it works. As always, thanks to Joe Z. and Arnauld C. for their invaluable assistance.

 

-dZ.

Edited by DZ-Jay
Link to comment
Share on other sites

Oh, and I forgot to post, this is the C routine from where I based the CP-1610 implementation. I think I found it in Wikipedia, under Linear-Feedback Shift Register (LFSR):

uint32_t lfsr = 1;
unsigned period = 0;

do {
  /* taps: 32 31 29 1; feedback polynomial: x^32 + x^31 + x^29 + x + 1 */
  lfsr = (lfsr >> 1) ^ (-(lfsr & 1u) ? 0xD0000001u : 0);
  
  ++period;
} while(lfsr != 1u);


uint32_t lfsr2(uint32_t lfsr, uint32_t iters)
{
   while (iters-- != 0)
       lfsr = (lfsr >> 1) ^ (lfsr & 1u ? 0xD0000001u : 0);

   return lfsr;
}

  • Like 1
Link to comment
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.
Note: Your post will require moderator approval before it will be visible.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

Loading...
  • Recently Browsing   0 members

    • No registered users viewing this page.
×
×
  • Create New...