Advanced Galton board simulator

[ParadigmShifter]

Probably advanced as in advanced maths though

You probably want to whack the emulator and processor speed up to max though after typing RUN. It's quite relaxing in a way I guess.

I don't think the Speccy RND function is very random I was expecting it to be much smoother than it is...




Bonus points if you know what it is doing.

Extra bonus points if you can make it run fast.

[C.Born]

a nummeric array should be 0 by default, that saves 3 x 265 steps
use a max varname length of 2 characters
start your prog with RANDOMIZE or RANDOMIZE 0

and i expected a triangle one

Spoiler

https://upload.wikimedia.org/wikipedia/ ... on_box.jpg

changed line a little and tryed to make a "run-save" time registry, bit to much, now you need a time routine !!
run 900

Code: Select all

EDITED
  30 LET z0=INT z0: IF z0<1 OR z0>256 THEN  GO TO 40
  35 LET a(z0)=a(z0)+1: IF a(z0)=177 THEN  GO TO 1200
  38 PLOT z0-1,a(z0)-1
  40 LET z1=INT z1: IF z1<1 OR z1>256 THEN  GO TO 50
  45 LET a(z1)=a(z1)+1: IF a(z1)=177 THEN  GO TO 1200
  48 PLOT z1-1,a(z1)-1
  50 GO SUB 230: GO TO 30
 220 LET ep=1e-37: LET tp=2*PI
 230 LET u1=RND: IF u1<=ep THEN  GO TO 230
 240 LET u2=RND: LET ma=si*SQR (-2*LN u1): LET z0=ma*COS (tp*u2)+mu: LET z1=ma*SIN (tp*u2)+mu: RETURN 
 900 CLEAR USR "a"-2: BORDER 7: CLS : FOR f=0 TO 168: POKE USR "a"-1+f,0: NEXT f: POKE 23672,0: POKE 23674,0: POKE 23673,0: POKE 23672,0
1000 RANDOMIZE 0: BORDER 4: DIM a(256)
1020 LET mu=128: LET si=128/3: GO SUB 220: GO TO 30
1200 BORDER 2: LET u=USR "a"-1: LET a=u+1+3*PEEK u: POKE u,1+PEEK u: POKE a+1,PEEK 23672: POKE a+2,PEEK 23673: POKE a+3,PEEK 23674: PAUSE 1e3: RUN 1e3
1234 CLEAR : SAVE "galtonbord" LINE 900
2345 REM paradigm shifter 2023

[C.Born]

to late to edit again and again so, i stop here, you make a timer yr self
its the same but working (on 16k aswell?)

Code: Select all

  30 LET z0=INT z0: IF z0<1 OR z0>256 THEN  GO TO 40
  35 LET a(z0)=a(z0)+1: IF a(z0)=177 THEN  GO TO 1200
  38 PLOT z0-1,a(z0)-1
  40 LET z1=INT z1: IF z1<1 OR z1>256 THEN  GO TO 50
  45 LET a(z1)=a(z1)+1: IF a(z1)=177 THEN  GO TO 1200
  48 PLOT z1-1,a(z1)-1
  50 GO SUB 230: GO TO 30
 220 LET ep=1e-37: LET tp=2*PI
 230 LET u1=RND: IF u1<=ep THEN  GO TO 230
 240 LET u2=RND: LET ma=si*SQR (-2*LN u1): LET z0=ma*COS (tp*u2)+mu: LET z1=ma*SIN (tp*u2)+mu: RETURN 
 900 CLEAR USR "a"-2: BORDER 7: CLS : LET f=USR "a"-1: FOR f=f TO f+21*8: POKE f,0: NEXT f
1000 RANDOMIZE 0: BORDER 4: DIM a(256)
1010 POKE 23672,0: POKE 23674,0: POKE 23673,0: POKE 23672,0
1020 LET mu=128: LET si=128/3: GO SUB 220: GO TO 30
1200 BORDER 2: LET u=USR "a"-1: LET a=u+1+3*PEEK u: POKE u,1+PEEK u: POKE a,PEEK 23672: POKE a+1,PEEK 23673: POKE a+2,PEEK 23674
1230 FOR f=31 TO 0 STEP -1: PRINT #0; PAPER 3;AT 0,0;TAB f; PAPER 4;" ": PAUSE 25: IF INKEY$="" THEN  NEXT f
1234 RUN 1e3
2000 CLEAR : SAVE "galtonbord" LINE 900
2345 REM paradigm shifter 2023

[flange]

Long time lurker but this caught my fancy so thought I'd dip a toe in the water.

For the ZX81 but should be transferable to the Spectrum, reasonably fast especially with ZX81 double speed BASIC ROM.

Code: Select all

  10 DIM A(10)
  20 FOR X=1 TO 7
  30 LET A(X)=0
  40 NEXT X
  50 PRINT AT 0,5;"BINOMIAL DISTRIBUTION"
  60 PRINT AT 5,15;"*"
  70 PRINT AT 7,14;"* *"
  80 PRINT AT 9,13;"* * *"
  90 PRINT AT 11,12;"* * * *"
 100 PRINT AT 13,11;"* * * * *"
 110 PRINT AT 15,10;"* * * * * *"
 120 PRINT AT 17,9;"-   -   -   -"
 130 PRINT AT 18,11;"-   -   -"
 200 PRINT AT 3,15;"O"
 210 LET X=15
 220 LET Y=3
 230 LET XP=X
 240 GOSUB 400
 250 IF Y=16 THEN GOTO 490
 260 LET R=RND
 270 IF R<0.5 THEN LET X=X+1
 280 IF R>=0.5 THEN LET X=X-1
 290 GOSUB 400
 300 GOTO 230
 400 PRINT AT Y,XP;" "
 410 LET Y=Y+1
 420 PRINT AT Y,X;"O"
 430 RETURN 
 490 PRINT AT 16,X;" "
 500 LET Z=(X-7)/2
 510 LET A(Z)=A(Z)+1
 520 LET P=17
 530 IF Z/2=INT (Z/2) THEN LET P=18
 540 LET Q=X
 550 IF A(Z)>=10 THEN LET Q=Q-1
 560 PRINT AT P,Q;A(Z)
 570 IF A(Z)=50 THEN INPUT Z$
 580 GOTO 200

[ParadigmShifter]

I was definitely expecting a smoother bell curve so either my algorithm implementation is incorrect or the random numbers are bad (I tried to use a better random routine with 49 bytes of ASM code but that didn't make much difference and I have to call it at least twice per pair of samples since it only generates 8 bits at a time so it ended up slower I think).

My algorithm is supposed to be an implementation of this:

https://en.wikipedia.org/wiki/Box%E2%80 ... _transform

Based on the C++ implementation

Code: Select all

#include <cmath>
#include <limits>
#include <random>
#include <utility>

//"mu" is the mean of the distribution, and "sigma" is the standard deviation.
std::pair<double, double> generateGaussianNoise(double mu, double sigma)
{
    constexpr double epsilon = std::numeric_limits<double>::epsilon();
    constexpr double two_pi = 2.0 * M_PI;

    //initialize the random uniform number generator (runif) in a range 0 to 1
    static std::mt19937 rng(std::random_device{}()); // Standard mersenne_twister_engine seeded with rd()
    static std::uniform_real_distribution<> runif(0.0, 1.0);

    //create two random numbers, make sure u1 is greater than epsilon
    double u1, u2;
    do
    {
        u1 = runif(rng);
    }
    while (u1 <= epsilon);
    u2 = runif(rng);

    //compute z0 and z1
    auto mag = sigma * sqrt(-2.0 * log(u1));
    auto z0  = mag * cos(two_pi * u2) + mu;
    auto z1  = mag * sin(two_pi * u2) + mu;

    return std::make_pair(z0, z1);
}

I did check that the epsilon is good enough so that log(epsilon) does not cause a domain error, not sure if that can happen since I think the range for speccy rand is [0, 65535/65536]?

I should cache 2*pi*u2 since it is used twice as well.

Wondering why it is so jaggedy though I was expecting a nice smooth curve :/

I chose sigma = 128 (middle of screen) and standard deviation so that edge of the screen was +/- 3 std deviations away.

This is one project that is much easier in BASIC than ASM cos of the floating point maths and LN and trig.

I also have a Poisson distribution sampler (it's slow though and breaks - because it is adding teeny weeny floats to a larger float so loses accuracy and may never return an answer - which is not good - if you need loads of iterations).

I have a C# version of an exponential distribution sampler as well (which I actually used in production code for weather [rainfall] simulation).

Obviously a Normal distribution sampler is a super handy piece of code, last time I used one in production code (to give reaction times etc. for AI) it used a big lookup table, so the B-M transform is a big improvement on that, glad I found it (even though it didn't actually get used).

I also converted a golf game (to PS1) which used a Normal distribution to add random noise to CPU player shots (otherwise they got a hole in one every par 3 lol). That also used a lookup table approach.

[ParadigmShifter]

Long time lurker but this caught my fancy so thought I'd dip a toe in the water.

For the ZX81 but should be transferable to the Spectrum, reasonably fast especially with ZX81 double speed BASIC ROM.

Code: Select all

  10 DIM A(10)
  20 FOR X=1 TO 7
  30 LET A(X)=0
  40 NEXT X
  50 PRINT AT 0,5;"BINOMIAL DISTRIBUTION"
  60 PRINT AT 5,15;"*"
  70 PRINT AT 7,14;"* *"
  80 PRINT AT 9,13;"* * *"
  90 PRINT AT 11,12;"* * * *"
 100 PRINT AT 13,11;"* * * * *"
 110 PRINT AT 15,10;"* * * * * *"
 120 PRINT AT 17,9;"-   -   -   -"
 130 PRINT AT 18,11;"-   -   -"
 200 PRINT AT 3,15;"O"
 210 LET X=15
 220 LET Y=3
 230 LET XP=X
 240 GOSUB 400
 250 IF Y=16 THEN GOTO 490
 260 LET R=RND
 270 IF R<0.5 THEN LET X=X+1
 280 IF R>=0.5 THEN LET X=X-1
 290 GOSUB 400
 300 GOTO 230
 400 PRINT AT Y,XP;" "
 410 LET Y=Y+1
 420 PRINT AT Y,X;"O"
 430 RETURN 
 490 PRINT AT 16,X;" "
 500 LET Z=(X-7)/2
 510 LET A(Z)=A(Z)+1
 520 LET P=17
 530 IF Z/2=INT (Z/2) THEN LET P=18
 540 LET Q=X
 550 IF A(Z)>=10 THEN LET Q=Q-1
 560 PRINT AT P,Q;A(Z)
 570 IF A(Z)=50 THEN INPUT Z$
 580 GOTO 200

That's a good method but does not scale well for large numbers of pegs, better then to use the Normal approximation (which is very accurate once you use enough pegs). Summing dice rolls to get a binomial approximation soon becomes very inefficient to do runtime.

[C.Born]

RND on zx is not very strong, its random sequence of ROM peeks and so it depends on the randomness off the rom wich aint, its a program and has looads of repeats in code and many call adresses in lower memory, the rom itself.
hence a bad picture. i tried myself in asm some things with the R register but thats not very strong to.

oops, timings

Spoiler

Code: Select all

  30 LET z=INT z: IF z<1 OR z>256 THEN  GO TO 40
  35 LET a(z)=a(z)+1: IF a(z)=177 THEN  GO TO 1200
  38 PLOT z-1,a(z)-1
  40 LET y=INT y: IF y<1 OR y>256 THEN  GO TO 50
  45 LET a(y)=a(y)+1: IF a(y)=177 THEN  GO TO 1200
  48 PLOT y-1,a(y)-1
  50 GO SUB 230: GO TO 30
 220 LET e=1e-37: LET t=2*PI
 230 LET u=RND: IF u<=e THEN  GO TO 230
 240 LET v=RND: LET n=s*SQR (-2*LN u): LET z=n*COS (t*v)+m: LET y=n*SIN (t*v)+m: RETURN 
 900 CLEAR USR "a"-2: BORDER 7: CLS : LET f=USR "a"-1: FOR f=f TO f+21*8: POKE f,0: NEXT f
1000 RANDOMIZE 0: BORDER 4: INPUT "": DIM a(256): LET cp=0
1010 POKE 23672,0: POKE 23674,0: POKE 23673,0: POKE 23672,0
1020 LET m=128: LET s=128/3: GO SUB 220: GO TO 30
1200 BORDER 2: LET u=USR "a"-1: LET a=u+1+3*PEEK u: POKE u,1+PEEK u: POKE a,PEEK 23672: POKE a+1,PEEK 23673: POKE a+2,PEEK 23674
1220 LET fr=PEEK a+256*PEEK (a+1)+65535*PEEK (a+2)
1221 PRINT #0;AT 0,0; PAPER 3;fr;TAB 31;" "
1222 LET s=fr/50: LET se=INT s: LET fr=fr-(se*50)
1223 LET m=se/60: LET mi=INT m: LET se=se-(mi*60)
1224 LET h=mi/60: LET ho=INT h: LET mi=mi-(ho*60)
1225 LET d=ho/24: LET dy=INT d: LET ho=ho-(dy*24)
1226 PRINT #0; PAPER 2;AT 0,16;d;"d";ho;"h";mi;"m";se;"s";fr;"f": REM max 3.2 days
1230 FOR f=31 TO 0 STEP -1: PRINT #0; OVER 1; PAPER 3;AT 1,0;TAB f; PAPER 4;" ": PAUSE 25: IF INKEY$="" THEN  NEXT f: IF cp THEN  COPY 
1234 BORDER 4: RUN 1e3
2000 CLEAR : SAVE d*"galtonbord" LINE 900
2345 REM paradigm shifter 2023

[ParadigmShifter]

It doesn't use ROM peeks it uses an LCG. Formula for it is in the manual, it only generates ints in range 0-65535 and divides by 65536 though I believe.

https://worldofspectrum.net/ZXBasicManu ... hap11.html

Looks like formula for next random number is

RND = ((75 * (seed + 1) - 1)&65535)/65536

LCGs are usually pretty bad though.

This is better but I think it was slower since I have to call it twice (at least).

Code: Select all

; 8-bit Complementary-Multiply-With-Carry (CMWC) random number generator.
; Created by Patrik Rak in 2012, and revised in 2014/2015,
; with optimization contribution from Einar Saukas and Alan Albrecht.
; See http://www.worldofspectrum.org/forums/showthread.php?t=39632

rnd:
	IF RAND_ALWAYS_1
	ld a, 1
	ret
	ENDIF
		ld   hl, table

        ld   a, (hl) ; i = ( i & 7 ) + 1
        and  7
        inc  a
        ld   (hl), a

        inc  l      ; hl = &cy

        ld   d, h    ; de = &q[i]
        add  a, l
        ld   e, a

        ld   a, (de) ; y = q[i]
        ld   b, a
        ld   c, a
        ld   a, (hl) ; ba = 256 * y + cy

        sub  c      ; ba = 255 * y + cy
        jr   nc, $+3
        dec  b

        sub  c      ; ba = 254 * y + cy
        jr   nc, $+3
        dec  b

        sub  c      ; ba = 253 * y + cy
        jr   nc, $+3
        dec  b

        ld   (hl), b ; cy = ba >> 8, x = ba & 255
        cpl         ; x = (b-1) - x = -x - 1 = ~x + 1 - 1 = ~x
        ld   (de), a ; q[i] = x

        ret

randSeed:
table   db   0, 0, 82, 97, 120, 111, 102, 116, 20, 15

        if   (table/256)-((table+9)/256)
        error "whole table must be within single 256 byte block"
        endif

And here is a POKE version of that (although it ran a few iterations before I dumped the memory and made it into a series of POKEs



Although there's code there to randomly plot stuff to the screen as well, which can be removed. Don't think it's relocatable either.

EDIT: Using R register is very bad since it only goes up to 127 before resetting, it's ok to use it to seed or jumble the seed a bit in ASM once in a while.

[C.Born]

seed is filled with RANDOMIZE number, with a 0 it with set some number from the ROM afaik
i have to check the rom for that

https://github.com/ZXSpectrumVault/rom- ... trum48.asm

Code: Select all

; ------------------------
; Handle RANDOMIZE command
; ------------------------
; This command sets the SEED for the RND function to a fixed value.
; With the parameter zero, a random start point is used depending on
; how long the computer has been switched on.

;; RANDOMIZE
L1E4F:  CALL    L1E99           ; routine FIND-INT2 puts parameter in BC.
        LD      A,B             ; test this
        OR      C               ; for zero.
        JR      NZ,L1E5A        ; forward to RAND-1 if not zero.

        LD      BC,($5C78)      ; use the lower two bytes at FRAMES1.

;; RAND-1
L1E5A:  LD      ($5C76),BC      ; place in SEED system variable.
        RET                     ; return to STMT-RET

it is the FRAMES counter after all, no roms

[ParadigmShifter]

Yeah it starts at a fixed seed and if you RANDOMIZE or RANDOMIZE 0 it sets it from FRAMES. You always get the same sequence if you don't do that or set seed manually to a non-zero value.

My guess is that the floating point random number is losing precision bits when it goes through LN, SQR and then SIN or COS.

I never tried it on a PC since I don't know how to draw pixels to the screen easily without using a massive framework (and even then I'm not really a graphics guy).

I did verify the mean and unbiased standard deviation estimate of the algorithm came out correct for a large number of samples... was going to do that on the speccy but writing something to collect more than 1 statistic is painful in basic with no local variables or structs

I wasn't too worried about precision for the algorithm for how I intended to use it (which was to do simulations where extreme accuracy was not required).

Would like to see how it looks on a PC and whether double precision is required or not to get good results though!

[C.Born]

had to check better
occurence
255=1242 times
154=1 time
good chance!! good that its the framecounter

Code: Select all

  20 DIM a(256): FOR f=0 TO 16383: LET a(PEEK f+1)=a(PEEK f+1)+1: NEXT f: FOR f=1 TO 256: PRINT f-1;"=";a(f): NEXT f


0=535
1=243
2=193
3=149
4=144
5=108
6=138
7=94
8=127
9=98
10=63
11=68
12=67
13=110
14=94
15=115
16=156
17=81
18=46
19=48
20=33
21=53
22=96
23=45
24=196
25=103
26=40
27=72
28=76
29=18
30=63
31=73
32=361
33=139
34=131
35=272
36=62
37=46
38=21
39=23
40=259
41=52
42=166
43=125
44=45
45=77
46=18
47=37
48=172
49=80
50=38
51=36
52=57
53=42
54=94
55=68
56=207
57=7
58=82
59=42
60=84
61=53
62=132
63=33
64=53
65=50
66=101
67=39
68=84
69=54
70=52
71=53
72=19
73=45
74=11
75=44
76=20
77=44
78=85
79=106
80=38
81=21
82=69
83=57
84=44
85=40
86=47
87=39
88=13
89=27
90=17
91=37
92=412
93=37
94=42
95=34
96=16
97=51
98=20
99=17
100=11
101=59
102=26
103=29
104=22
105=36
106=8
107=15
108=29
109=14
110=67
111=53
112=31
113=37
114=58
115=38
116=61
117=23
118=24
119=70
120=78
121=53
122=39
123=25
124=36
125=21
126=132
127=17
128=46
129=13
130=19
131=6
132=6
133=5
134=26
135=11
136=14
137=17
138=38
139=7
140=6
141=13
142=10
143=18
144=20
145=21
146=13
147=9
148=13
149=11
150=14
151=7
152=10
153=21
154=1
155=8
156=5
157=3
158=6
159=21
160=25
161=34
162=19
163=8
164=11
165=6
166=10
167=78
168=16
169=17
170=11
171=19
172=10
173=11
174=30
175=50
176=34
177=29
178=25
179=15
180=12
181=16
182=28
183=4
184=27
185=15
186=10
187=5
188=7
189=6
190=25
191=14
192=74
193=105
194=41
195=91
196=32
197=101
198=43
199=7
200=51
201=206
202=31
203=286
204=22
205=591
206=36
207=46
208=30
209=68
210=22
211=23
212=31
213=80
214=37
215=29
216=35
217=132
218=25
219=16
220=12
221=83
222=12
223=61
224=41
225=135
226=20
227=35
228=12
229=151
230=85
231=97
232=15
233=23
234=5
235=149
236=3
237=162
238=49
239=90
240=15
241=85
242=19
243=15
244=27
245=55
246=31
247=21
248=24
249=15
250=27
251=43
252=13
253=219
254=281
255=1242

[ParadigmShifter]

Yeah the ROM is not a good source for uniform random data, although there are patches where it isn't too bad if you are doing something that does not require nice statistical properties for your random numbers (so an arcade game or something).

If you are doing a simulation you want to make sure your RNG has a uniform distribution (which the speccy RND does have, but it has too much correlation which ends up with visible artifacts like diagonal bands etc.).

EDIT: Responding to deleted post which showed histogram of the distribution of bytes in the ROM

EDIT: There's a massive block of contiguous 255s in the ROM somewhere as well, avoid that

EDIT2:

https://skoolkid.github.io/rom/dec/asm/14446.html

1170 bytes of 255 starts at address 14446. Think that got replaced in the 128K ROM (people used to use it for an ISR vector table but that broke in 128K mode).

EDIT3: Looks like you have an off-by-one error in your code, 0 is PEEKable and 16384 is screen RAM

[ParadigmShifter]

Changed all the variables to single letters and did the increment after the plot to save some subtractions. Also multiplied v (was u2) by 2*PI so it only does that once.

Seems a lot faster at max CPU speed anyway




This is how jaggedy it looks though



EDIT: Brackets for SIN and COS can go as well, can just use k*SIN v+m, k*COS v+m and precedence handles it properly

EDIT2: Then I did

Code: Select all

230 LET u=RND: IF u<=e THEN GO TO 230
240 LET u=-2*LN u: LET v=RND*t: LET k=s*SQR u: LET z=k*COS v+m: LET w=k*SIN v+m: RETURN

EDIT3: I'll remove the GO SUB next and rejink the main loop a bit

EDIT4: Painless and fits on 1 page now

[ParadigmShifter]

Added a range scaler (r, and q=1/r) to see if that makes it smoother (currently trying 32) which should tell me whether there's some kind of bias after the calculations or not.



I'll post an image when it finishes lol

EDIT: First I did a cheeky LET r = r * 176 to avoid 2 multiplies in the mainloop after it calculates q

EDIT2: Pretty sure the seed will wrap around during this run though (if I'm lucky it will wrap around after an odd number of in range samples though so the period might be doubled?). I think the period is only 65536 (since that's the range of the seed anyway). Probably need more precision random numbers (at least 24 bits).

EDIT3: Results, dunno if it is a random precision/cycle length issue or cumulative error/loss of precision from the maths stuff



EDIT4: I expect SIN and COS get quantised a lot and lose a ton of precision if passed arguments way outside [-2PI, 2PI] which may happen I guess

[C.Born]

they have bigger PI here
https://rosettacode.org/wiki/Pi

[ParadigmShifter]

10 digits of precision is easy enough for PI (so a double these days), no point being more accurate than that.

There is an algorithm which will produce the hexadecimal digits of PI though

https://en.wikipedia.org/wiki/Bailey%E2 ... fe_formula

ZX Basic uses PI/2 to calculate PI I think... and looks like it stores it in a compressed form to make it extra slow for no reason, considering they wasted those 1K+ bytes of 255, that was because they took it from the ZX81 ROM I think

https://skoolkid.github.io/rom/dec/asm/12997.html

atan(1) = PI/4 anyway which is what most people use if PI is not available

so

const double PI = 4 * atan(1); // accidentally had int before lol, so PI = 3 then, which makes stuff easier but sadly less accurate

[ParadigmShifter]

Improving the RNG seems to make it better, this is scaled down 8 (so 8 marbles per bin for 1 pixel)



Code, can pick between 16 and 32 bits of precision by REMing out the one you don't want




RNG code which ended up at address 37943 this time

Code: Select all

basic_rnd16:
	call rnd
	ld c, a
	push bc
	call rnd
	pop bc
	ld b, a
	ret

basic_rnd:
	call rnd
	ld b, 0
	ld c, a
	ret

; 8-bit Complementary-Multiply-With-Carry (CMWC) random number generator.
; Created by Patrik Rak in 2012, and revised in 2014/2015,
; with optimization contribution from Einar Saukas and Alan Albrecht.
; See http://www.worldofspectrum.org/forums/showthread.php?t=39632

rnd:
	IF RAND_ALWAYS_1
	ld a, 1
	ret
	ENDIF
		ld   hl, table

        ld   a, (hl) ; i = ( i & 7 ) + 1
        and  7
        inc  a
        ld   (hl), a

        inc  l      ; hl = &cy

        ld   d, h    ; de = &q[i]
        add  a, l
        ld   e, a

        ld   a, (de) ; y = q[i]
        ld   b, a
        ld   c, a
        ld   a, (hl) ; ba = 256 * y + cy

        sub  c      ; ba = 255 * y + cy
        jr   nc, $+3
        dec  b

        sub  c      ; ba = 254 * y + cy
        jr   nc, $+3
        dec  b

        sub  c      ; ba = 253 * y + cy
        jr   nc, $+3
        dec  b

        ld   (hl), b ; cy = ba >> 8, x = ba & 255
        cpl         ; x = (b-1) - x = -x - 1 = ~x + 1 - 1 = ~x
        ld   (de), a ; q[i] = x

        ret

randSeed:
table   db   0, 0, 82, 97, 120, 111, 102, 116, 20, 15

        if   (table/256)-((table+9)/256)
        error "whole table must be within single 256 byte block"
        endif

I'll do another run with a larger bin size I expect it still won't be super smooth in which case it's probably the precision of the floating point maths (although speccy has 4 bytes of mantissa so better precision than a 32 bit float but not as good as a double)

EDIT: Well at least there's a better version of RND for you all anyway Probably slower than built in RND

EDIT2: If you want to do a RANDOMIZE though you have to POKE the 3 bytes of the SYSVAR FRAMES into any of the addresses following the last 201 which is where the seed table is stored.

[ParadigmShifter]

Scale 32, pretty smooth then considering the number of bins I think



So I think it is mainly the RND function which needed to be better

This is what it should converge towards





That function isn't very useful unless you want to draw the graph (or prove the Central Limit Theorem of course)... it's the area under the graph that you want which is a non-standard function (erf(x) in maths libraries I think). Non-standard = not expressible exactly as a formula in terms of standard operations (+, -, *, /, nth roots, exponentiation, logarithms, trigonometry).

EDIT: Originally had 1/(SQR PI), should be 1/(SQR(2*PI)). Looks the same though to me lol Working out the max and scaling cancelled that out. Fixed the listing and the output anyway

[ParadigmShifter]

So after doing an overnight run I think with the better version of RND everything is good so the maths functions are fine.

Used a scale of 256 (so 256 marbles per bin for 1 pixel). Took hours to run even at max emulation CPU speed multiplier

[C.Born]

for better or worse,:

my "regr + frames implementation" attempt:

Code: Select all

ld hl,(23672)
ld a,r     ; fetch %0xxxxxxx byte
rrca       ; move bit 0 in bit 7 , sets 'gab' at bit 6 %x0xxxxxx
xor h
ld h,a
ld a,r
rrca
xor l
ld l,a
ld (23670),hl   ; seems ignored with out loading BC with HL ?????????????????????????????
push hl
pop bc
ret