We're on ANOTHER random number generation technique today. Let's work on the Lagged Fibonacci Sequence, a fun way to use the Fibonacci Sequence on itself to generate pseudo-random numbers.
One of our long term goals is to add a bunch of tools to our toolbox. Making reusable components makes doing more complicated things easier. You can rifle through your toolbox in hopes of an "Oh, I already have something we can use for this part!".
What's the Fibonacci sequence again?
It's a simple pattern that appears in a kinds of places in nature. For example, tree branches follow a Fibonacci pattern. Who knew trees were good at math?
In a series, the next value is the sum of the previous two.
So if we start with 0 and 1, next is... 1. Then...
Let's do this on a VIC-20
The lowly VIC 20 doesn't get much love since it's often overshadowed by the Commodore 64. What most people don't know is that for raw calculations, it's actually faster than a 64. Yes there are benchmarks later, calm down.
We need to get two values and add them to get a third. Lets try to use GOSUBs to organize our code so we can add stuff to our toolbox. You never know when you might need a Fibonacci series generator in BASIC. Could happen.
10 F=0:FOR C=0 TO 20 20 GOSUB 100 30 PRINT F; 40 NEXT 50 END 100 REM *** FIBONOCCI *** 110 IF C = 0 THEN F=0:RETURN 120 F1=0:F2=1 130 FOR I=2 TO C 140 F=F1+F2:F1=F2:F2=F 150 NEXT 160 RETURN
Used most famously in FreeCiv, but even Oracle databases use it in some places, the Lagged Fibonacci Generator is another way to generate pseudo random numbers with okayish randomness properties.
Directly from Wikipedia:
For the generator to achieve this maximum period, the polynomial:
y = xk + xj + 1
must be primitive over the integers mod 2. Values of j and k satisfying this constraint have been published in the literature.
Clear as mud? It'll be easier to visualize in the code.
We're going to use
k=3 j=7 for this example to keep this short. Feel free to use our example here and expand it to various keys and test your result.
10 DIMS(8):S$="8675309":F=0:MAX=255:J=3:K=7 15 FORI=1TO7:S(I)=VAL(MID$(S$,I,1)):NEXT 20 GOSUB 100 100 REM *** LAGGED FIB *** 110 F=S(J)+S(K):IFF>255THENF=F-MAX+1 120 FORI=1TO6:S(I)=S(I+1):NEXT:S(7)=F 130 RETURN
If we run this 10 times we get some random numbers:
Monte Carlo Method to test randomness
We'll do a bunch of iterations over this subroutine in a second. First we should talk about how random this may or may not be.
As you might expect, there are many ways to test if a series is random. The Monte Carlo of Pi method is an excellent way to inefficiently estimate Pi, but another use is to test for randomness. The closer your results are to pi, the more random your series is.
We will take our random numbers and scale them between 0.0 and 1.0, then using two values at a time and calculate the square root of x2+y2 . If this value is less than or equal to 1, we place in the circle (with a radius of 1), otherwise it is outside the circle. The estimation of PI is then four times the number of points in the circle (M) divided by the total number of points (N).
We'll add it as a GOSUB starting at line 300. Lines 345-390 do a little printing and formatting to help see what's happening.
1 C=100 5 TI$="000000" 10 J=3:K=7:DIMS(8):S$="8675309":F=0:MAX=255:IN=0:OUT=0 15 FORI=1TO7:S(I)=VAL(MID$(S$,I,1)):NEXT 20 DIMV(C):FOR N=0TOE 30 GOSUB 100 50 V(N)=F:PRINTV(N); 60 NEXT 65 GOSUB 300 70 PRINT:PRINT"JIFFIES:"TI 80 END 100 REM *** LAGGED FIB *** 110 F=S(J)+S(K):IFF>255THENF=F-MAX+1 120 FORI=1TO6:S(I)=S(I+1):NEXT:S(7)=F 130 RETURN 200 PRINT"->";:FORP=1TO7:PRINTS(P);:NEXT:PRINT: RETURN 300 REM *** Monte Carlo *** 305 PRINT:FORN=0TOC-1 STEP2 310 X=(V(N)-MAX/2)/(MAX/2) 320 Y=(V(N+1)-MAX/2)/(MAX/2) 330 Z=SQR(X^2+Y^2) 340 X$=STR$(X):Y$=STR$(Y):Z$=STR$(Z) 345 PRINTLEFT$(X$,4);:PRINTLEFT$(Y$,4);:PRINTLEFT$(Z$,4) 350 IF Z<1 THENIN=IN+1 360 IF Z>1 THENOUT=OUT+1 365 NEXT 370 PRINT:PRINT"IN:"IN"OUT:"OUT 380 PRINT:PRINT"RESULT:"4*IN/(IN+OUT) 390 RETURN
Embracing the constraints
We're doing this on an unexpanded VIC-20 on purpose. Partly because it's cool and partly because it almost immediately shows us a problem with our implementation above. Constraints can be a good thing.
This implementation so far reflects my modern approach to solving this problem where memory is cheap and plentiful.
My daily driver computer is a Mac with 3.2x1010 Bytes of RAM (32GB).
The VIC 20 has 3583 Bytes. Slightly less.
On my Mac, I don't think much of making an array of 50k values. On a VIC-20, even a 100 value array is pretty large. Remember, we have 3583 Bytes to work with. A DIM(100) takes up something like 500 Bytes on a Commodore plus the values in the DIM.
So what do we do about this if we want to calculate this to 50k pairs?
We have limited RAM, but "unlimited" CPU time if we're willing to wait on a 1 Mhz CPU. So let's rethink this:
1 C=100 5 TI$="000000" 10 J=3:K=7:DIMS(8):S$="8675309":F=0:MAX=255:IN=0:OUT=0 15 FORI=1TO7:S(I)=VAL(MID$(S$,I,1)):NEXT 20 FOR N=0TOC 30 GOSUB 100 40 FX=F:PRINT"NEXT PAIR:";:PRINTFX; 45 GOSUB 100 50 FY=F:PRINTFY; 60 GOSUB 300 65 GOSUB400 70 NEXT 75 PRINT:PRINT"JIFFIES:"TI 80 END 100 REM *** LAGGED FIB *** 110 F=S(J)+S(K):IFF>255THENF=F-MAX+1 120 FORI=1TO6:S(I)=S(I+1):NEXT:S(7)=F 130 RETURN 200 PRINT"->";:FORP=1TO7:PRINTS(P);:NEXT:PRINT: RETURN 300 REM *** MONTE CARLO *** 305 PRINT 310 X=(FX-MAX/2)/(MAX/2) 320 Y=(FY-MAX/2)/(MAX/2) 330 Z=SQR(X^2+Y^2) 340 X$=STR$(X):Y$=STR$(Y):Z$=STR$(Z) 345 PRINT"X Y Z:";:PRINTLEFT$(X$,4);:PRINTLEFT$(Y$,4);:PRINTLEFT$(Z$,4); 350 IF Z<1 THENIN=IN+1 360 IF Z>1 THENOUT=OUT+1 365 RETURN 400 PRINT:PRINT"IN:"IN"OUT:"OUT"TOTAL:"IN+OUT 410 PRINT"RESULT:"4*IN/(IN+OUT):PRINT:PRINT 420 RETURN
Much better. Now we're limited by the CPU speed and wont run out of RAM.
We can't do these without benchmarking on different systems can we?
Lets skip the printing along the way and just loop through a bunch of them:
Monte Carlo Visual
It's fun to visualize the dots inside and outside the circle. If we expand our constraints just a little and add the Super Expander Cartridge to our VIC 20, doing a visual is pretty easy.
I don't want to be called lazy, especially on the internet where words mean the difference between life and death.
Or is it "the angry internet people don't matter". Anyway, it's one of those...
If we want to do any large calculations of this or use it in a game we need to implement it in Assembly anyway.
Let's make a routine that's as close to the BASIC version as we can. It's some simple 16 bit addition which we've done before:
fib ldx #0 lda #0 jsr $ddcd lda #32 jsr $ffd2 ldx f2 lda f2+1 jsr $ddcd lda #32 jsr $ffd2 loop clc lda f1 adc f2 sta f0 lda f1+1 adc f2+1 sta f0+1 lda f2 sta f1 lda f2+1 sta f1+1 lda f0 sta f2 lda f0+1 sta f2+1 print ldx f0 lda f0+1 jsr $ddcd lda #32 jsr $ffd2 inc c ldx c cpx #23 bne loop rts f0 .byte 0,0 f1 .byte 0,0 f2 .byte 1,0 c .byte 0
We're going for as close to the BASIC implementation as possible here as well. This is about understanding, not breaking speed or efficiency records.
j .byte 2 k .byte 6 s .byte 8,6,7,5,3,0,9 f .byte 0 fx .byte 0 fy .byte 0 max = 255 slen = 6 ;zero based zp1 = $fb decout = $ddcd ; change to $bdcd on C64 lagfib .block loop ldx j ldy k lda s,x clc adc s,y cmp #max bcc ahead sec sbc #max ahead sta f ldx #0 loop2 lda s+1,x sta s,x inx cpx #slen bne loop2 lda f ldx #slen sta s,x print ldy #0 sty zp1 loop3 lda s,y tax lda #0 jsr decout lda #$20 jsr $ffd2 inc zp1 ldy zp1 cpy #slen+1 bne loop3 lda #$0d jsr $ffd2 done rts .bend
That's it. Now we have an okayish Pseudo Random Number generator that 8 Bit Computers can do fairly easily.
In your spare time it's fun took for places to optimize all of these loops both in the BASIC versions and Assembly. Try them out and see how fast you can make it. I want to hear about your improvements to what we've done here.