Random-number generation

[Arlet]

The underlying reason that you find the golden ratio in so many places:
https://www.youtube.com/watch?v=sj8Sg8qnjOg

[BigEd]

Very nice video! Continued fractions are very handy, but I didn't learn them in school (or at Uni), even though there's nothing very difficult there. A good example of why it's worthwhile continuing to learn things...

As for the Quasi Random sequences, I recall (re)inventing this idea of adding Phi back in the day, probably early 90s, when I used it to generate addresses (indexes into an array) to probe cache sizes on Sun's various models of Sparc workstation.

[Alex1]

Arlet, I agree with your explanation.

[Arlet]

I just made a github repository to store my code, for easy reference.

For now, it just has my latest 6502 assembly version that was posted here, but I'm planning to add more versions as they pass the appropriate tests, including for other kinds of CPUs, for those interested, and maybe some more 6502 variants for different size/speed vs quality criteria.

[Alex1]

Arlet, try to make a version where there is no pattern

[Arlet]

I can try, assuming "no pattern" means that "no 4-byte sequence should repeat until 4GB of data has been generated". The problem is combining that requirement with high quality random output, but I can try.

[barrym95838]

Arlet wrote:

I can try, assuming "no pattern" means that "no 4-byte sequence should repeat until 4GB of data has been generated". The problem is combining that requirement with high quality random output, but I can try.

I'm no expert on the subject, but I think those goals are mutually exclusive, as you seemed to indicate earlier.

[Arlet]

It should be possible to have a 4GB sequence, of which the first part is high quality. Obviously, if you have already generated 3.9GB, the next 0.1 GB won't be very random, because they are restricted to the sequences that haven't shown up yet. Just like a card counter at a casino would be able to predict the last few remaining cards in the deck.

[Arlet]

Ok. Here's the deal. It's not going to be possible to write a random generator where every 32 bit sliding window produces a unique value. The only way to do that is to use a 32 bit shift register, and then feed in new bits from one side. That's what a LFSR does. I have no doubt that there are more ways than just a LFSR to produce the same kind of sequence, but they aren't going to be a) simple, and b) produce good quality randomness.

Now, what is possible? We could make a random generator that produces 4 billion unique outputs, each 32 bit wide. The problem is that good quality will require a large amount of code & cycles, which I don't find very attractive or interesting.

We could also make a random generator with 32 bit state, producing 4GB worth of single bytes. That's a bit easier, and you can find code below.

I had posted this earlier, and then removed it because I found a bug. Now I'm putting it back because the bug doesn't seem to hurt the actual results. The bug is the presence of the ROL instruction near the end. The whole output function from s -> o was supposed to be a 32x32 bit bijective function, and that ROL breaks that. However, even with a non-bijective function, the actual output byte in the accumulator is still completely fair (every possible byte occurs exactly 2^24 times). I don't understand why.

I am posting this code just as a point of curiosity. It is slower and bigger than my earlier random generators, uses more RAM, and produces lower quality output. Because the 6502 has limited registers, we are forced to use memory for temporary variables, but if you're going to do that anyway, it makes more sense just to use a bigger state. With a bigger state, it's much easier to produce good quality randomness.

Code: Select all

CLC 
LDA #$45
ADC s0
STA s0
ADC s1
STA s1
ADC s2
STA s2
ADC s3
STA s3
EOR s2
ASL   
ADC s1
STA o3
ADC s2
STA o2
ADC s0
ASL  
ADC o2
STA o0
ADC s1
STA o1
ADC o3
EOR o2
ASL  
ROL  
ADC o1
ADC o2
ADC o0
RTS 

s0-s3 is a 32 bit internal state, o0-o3 are temporary variables to calculate the output.

[Alex1]

Interesting algorithm

[Arlet]

It uses the idea that the combination of

Code: Select all

ASL 
ADC o2

is equivalent to:

Code: Select all

ASL
ADC #0
CLC
ADC o2

That ASL + ADC #0 makes a rotate by 1 bit, which helps to mix the LSBs, which are the least random.

[Arlet]

Sneak preview of a new random number generator with 6 bytes of state in s0-s5, and 4 bytes of output in o0-o3. It takes a bit more code than previous designs, but you get 4 output bytes instead of 1, so it's quite efficient in cycles/byte if you need bigger outputs. It passes PractRand <= 256 GB, so quality is quite good. I am still working on improvements, but progress slows down with increasing data set (takes about 45 minutes on my machine to test 256 GB for a single code sequence)

Code: Select all

LDA #$45
CLC
ADC s0
STA s0
ADC s1
STA s1
ADC s2
STA s2
ADC s3
STA s3
ADC s4
STA s4
ADC s5
STA s5
LDA s2
ASL
ADC s5
STA o0
ADC s3
STA o1
ADC s4
STA o3
LDA o2
ASL
ADC o3
STA o2
ADC o0
STA o0
ADC o1
ASL
ADC o0
STA o1
ADC o3
STA o3
ADC o2
ASL
ADC o3
STA o2
ADC o1
STA o1
ADC o0
STA o0
RTS

I have not tested the 6502 code, but with all zero seed, it should produce ec 38 f4 64 9e c3 63 ec d0 01 89 f4 63 67 8c dc as output.

[Alex1]

It's also an interesting algorithm.

But how does Practrand work? I did a test on 256GB and it only does 195 tests. So Practrand does not test all the outputs of an algorithm but only does a few random tests at random?

The FIPS test, on the contrary, analyzes all the data.

Quote:

rng=RNG_stdin8, seed=unknown
length= 256 gigabytes (2^38 bytes), time= 110189 seconds
no anomalies in 195 test result(s)

[Arlet]

Quote:

So Practrand does not test all the outputs of an algorithm but only does a few random tests at random?

It means it does 195 different tests, all running in parallel, and each test gets a copy of the full input data. It prints out a status report every doubling of the data, but then it just continues with additional data to get a more precise answer.

BTW, when I run on 256GB, it says it's doing 372 tests, not 195. Maybe because I run 'stdin' instead of 'stdin8' ?

[Alex1]

I also give you my LCG 64 which dates from 1990,

It generates a random sequence of (2^64)-1 bytes. It's a 64-bit Linear Congruential Generator. Each call to the routine generates an 8-bit byte. There is no short cycle for all seed values.

First a 8-bytes seed must be put at seed:
The seed can vary from 0 to $FFFFFFFFFFFFFFFF.

Each call to the begin: function returns a pseudorandom byte in rndbyte:

First init: must be call one time.
After begin: must be call for each generated byte (rndbyte:) that can be use directly.

The test vector is :

For a seed = $FFFFFFFFFFFFFFFF

the first 50 bytes are :

Quote:

AE 0B D2 B3 F2 5B 8B 4B AD 23 58 5C 68 BB 62 AD 96 07 4F 7C E3 58 54 17 6F 3B EF 9A ED A5 19 A2 C9 6B 47 71 0F 3A 2A 23 38 F5 4F 34 23 9B BF 21 4C 03

Code: Select all

//----------------------------------------------------
//          Pukall LCG 64
//                1990
//----------------------------------------------------

//----------------------------------------------------
//          Main Program
//----------------------------------------------------

init:
    lda     #6
    sta     pos
    rts
begin:
    ldx     pos
    cpx     #6
    bne     rnd
    jsr     start
    ldx     #0
    clc
    php
loop0:
    lda     seed,x
    plp
    adc     dda,x
    php
    sta     seed,x
    inx
    cpx     #8
    bne     loop0
    plp
rnd:
    ldx     pos
    lda     seed,x
    sta     rndbyte
    dex
    cpx     #2
    bne     end2
end:
    lda     #6
    sta     pos
    rts
end2:
    stx     pos
    rts
start:
    lda     #0
    ldx     #15
loop:
    sta     result,x
    dex
    cpx    #7
    bne     loop
    ldy     #64       
part1:
    sty     storey
    ldx     #7
    lsr     seed,x
    php
    dex
loop1:
    plp
    ror     seed,x
    php
    dex
    cpx     #0
    bne     loop1
    plp
    ror     seed
    bcc     part2 
    ldx     #7
    ldy     #0
    clc
    php
loop2:
    inx
    lda     result,x
    plp
    adc     two,y    
    php
    sta     result,x
    iny
    cpy     #8
    bne     loop2
    plp
part2:
    ror     
    php
    ldx     #15         
    sta     result,x
loop3:
    dex
    plp
    ror     result,x
    php
    cpx     #0
    bne     loop3
    plp
    ldy     storey
    dey              
    bne     jump
    ldx     #0
loop4:
    lda     result,x
    sta     seed,x
    inx
    cpx    #8
    bne     loop4
    rts
jump:
    jmp part1

seed: .byte $FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF
rndbyte: .byte $FF
result: .byte $FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF,$FF
two: .byte $2D,$7F,$95,$4C,$2D,$F4,$51,$58
dda: .byte $01,$00,$00,$00,$00,$00,$00,$00
storey: .byte $FF
pos: .byte $FF