XProfan

WindowTitle "Algorithm 488: Gauss-Random Pseudozufallsgenerator"
' Fortran-Quelle und Urheberrechtsträger: Algorithm 488 of collected algorithms
' (C) ACM: https://www.acm.org/ , Algorithm appeared in comm. acm, vol. 17, no. 12, p. 704.
' Veröffentlicht in https://www.netlib.org/toms/
' Test for Migration to XProfan11.2a, 2014-10 by P.Specht, Wien; Ohne jede Gewähr!
WindowStyle 24:Window 0,0-%maxx,%maxy-40:randomize:font 2
var xh&=width(%hwnd)/2:var yh&=height(%hwnd)*9/10
usepen 0,1,0:line 0,yh& - 2*xh&,yh&:line xh&,0 - xh&,2*yh&
'
' Function GRand() 'Gauss-Rnd (Initialisierung ausgelagert als GRandInit)
'
' Except on the first call grand returns a pseudo-random number having a gaussian
' (i.e.normal) distribution with zero mean and unit standard deviation.
' Thus, the density is  f(x) = exp(-0.5*x**2)/sqrt(2.0*pi). the first call
' initializes grand and returns zero. The parameter n is dummy.
' gRand calls a function rand, and it is assumed that successive calls to rand(0)
' give independent pseudo-random numbers distributed uniformly on (0,1), possibly
' including 0 (but not 1). the method used was suggested by von neumann, and
' improved by forsythe, ahrens, dieter and brent.
' on the average there are 1.37746 calls of rand for each call of grand.
' Warning - dimension and data statements below are machine-dependent.
' Dimension of d must be at least the number of bits in the fraction of a
' floating-point number. Thus, on most machines the data statement below
' can be truncated.
' if the integral of sqrt(2.0/pi)*exp(-0.5*x**2) from
' a(i) to infinity is 2^(-i), then d(i) = a(i) - a(i-1).
GRandInit:
declare d$[],d![],u!:d$[]=explode( \
"0,0.674489750,0.475859630,0.383771164,0.328611323,0.291142827,0.263684322,"+\
"0.242508452,0.225567444,0.211634166,0.199924267,0.189910758,0.181225181,"+\
"0.173601400,0.166841909,0.160796729,0.155349717,0.150409384,0.145902577,"+\
"0.141770033,0.137963174,0.134441762,0.131172150,0.128125965,0.125279090,"+\
"0.122610883,0.120103560,0.117741707,0.115511892,0.113402349,0.111402720,"+\
"0.109503852,0.107697617,"+\
"0.105976772,0.104334841,0.102766012,0.101265052,0.099827234,0.098448282,"+\
"0.097124309,0.095851778,0.094627461,0.093448407,0.092311909,0.091215482,"+\
"0.090156838,0.089133867,0.088144619,0.087187293,0.086260215,0.085361834,"+\
"0.084490706,0.083645487,0.082824924,0.082027847,0.081253162,0.080499844,"+\
"0.079766932,0.079053527,0.078358781,0.077681899" , "," )

whileloop 0,60:d![&Loop]=val(d$[&Loop])'::print &Loop,format$("%g",d![&Loop])

    endwhile:clear d$[]'::waitinput
    ' end of machine-dependent statements, but:
    ' u must be preserved between calls!
    GLOCKENKURVE_DARSTELLEN:
    Declare cnt&,grnd!,idx&,h&[2*xh&],diehöllezufriert&

    Whileloop 200000:cnt&=&Loop

        grnd!=GRand()
        idx&=xh&+(xh&*grnd!/5)
        h&[idx&]=h&[idx&]+1

        if abs(Grnd!)<0.002

            locate 2,2:Print cnt&,tab(10);format$("%g",grnd!);"    ";:moveto 0,yh&
            usepen 0,1,rgb(rnd(255),rnd(255),h&[idx&])'cnt&,0,h&[idx&])

            whileloop 0,2*xh&:lineto &Loop,yh&-h&[&Loop]

            endwhile

        endif

    Endwhile

    beep
    locate 2,2:Print cnt&,tab(10);format$("%g",grnd!);"    ";
    waitinput
    end

    proc GRand

        declare a!,i&,v!,w!,grand!
        ' initialize displacement a and counter i.
        a! = 0.0
        i& = 0
        ' increment counter and displacement if leading bit of u is one.
        g10:
        u!=u!+u!
        case u!<1:goto "g20"
        u!=u!-1
        inc i&
        a!=a!-d![i&]
        goto "g10"
        ' form w uniform on 0 < w < d(i+1) from u.
        g20:
        w! = d![i&+1]*u!
        ' form v = 0.5*((w-a)**2 - a**2). note that 0 < v < log(2).
        v! = w!*(0.5*w!-a!)
        ' generate new uniform u.
        g30:
        u! = rnd()
        ' accept w as a random sample if v! < u!
        case v!<u!: goto "g40"
        ' generate random v.
        v! = rnd()
        ' loop if u .gt. v.
        case u!>v!:goto "g30"
        ' reject w and form a new uniform u from v and u.
        u! = (v!-u!)/(1-u!)
        goto "g20"
        ' form new u (to be used on next call) from u and v.
        g40:
        u! = (u!-v!)/(1-v!)
        ' use first bit of u for sign, return normal variate.
        u!=u!+u!
        case u!<1:goto "g50"
        u! = u! - 1
        grand! = w!-a!
        return grand!
        g50:
        grand! = a! - w!
        return grand!

    endproc