; --------------- Alea Test --------------------
;       test your random number generator
; -------- guillaume.tello@orange.fr -----------

; procedures

#label D random
#label E init
#label A test

; variables

#reg 00 counter
#reg 00 average
#reg 11 index
#reg 12 total(x)
#reg 13 total(x)
#reg 14 histogram
#reg 15 seed
#reg 16 number(x)

; reg 01 to reg 10 contains how many numbers
; falled into each of the ten parts 0.0 0.1 0.2 0.3 0.4 ... 0.9

;-----------------
; E to init a test
;-----------------

LBL init
     seed={ @seed CMS }                 ; clear all but the seed
RTN

;--------------------------------------------
; n A to add n tests to the current histogram
;--------------------------------------------

LBL test
    STO counter                         ; store n
    FOR counter
        total(x)+{ SBR _random }        ; accumulate every number
        index={ * $10 + 1 = INT }       ; wich of the 10 parts?
        1 SM* index                     ; one more number in this slice
        total(x)+{ @seed X }          ; accumulates the squares
    LOOP

    counter={ $10 }                     ; 10 for the loop
    X/T                                 ; 10 saved into T
    number(x)={
	      histogram,index={ 0 }       ; uses first zero to clear histogram and index	
            FOR counter                 ; computes the total number of tests
                + RC* counter
            LOOP
            =
              }

    counter={ X/T }                     ; gets back 10 for the loop
    X/T                                 ; and saves it again to compare

    FOR counter                         ; loop through the 10 parts
        index+{ 1 }                     ; index to the correct register
        histogram*{ $10 }               ; shift histogram left
        histogram+{                     ; and adds one more digit
                    RC* index * $50     
                    / @number(x)
                    + $.5 = INT         ; the digit!
                    IF(x>=t)            ; if >= 10
                        9               ; then limit to 9
                    ENDIF   
                  }
                  
    LOOP
    
    RCL histogram R/S                   ; display histogram
                                        ; then R/S display the average value
    average={ @total(x) / @number(x) = }  
        R/S
                                        ; then R/S and you finally get the standard deviation
    { @total(x) / @number(x) - @average x = } SQR

RTN

LBL random
_random:
    seed={ @seed + PI = Y^X 5 = frac }
RTN
