XProfan

Titre de la fenêtre "Algorithm 488: Gauss-Random Pseudozufallsgenerator"
' Fortran-source et Urheberrechtsträger: Algorithm 488 of collected algorithms
' (C) ACM: https://www.acm.org/ , Algorithm appeared dans comm. acm, vol. 17, no. 12, p. 704.
' Veröffentlicht dans https://www.netlib.org/toms/
' Test for Migration to XProfan11.2a, 2014-10 by P.Specht, Wien; sans chacun Gewähr!
Fenêtre Style 24:Fenêtre 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-Tour (initialisation ausgelagert comme GRandInit)
'
' Except on le first call grand returns a pseudo-random number having a gaussian
' (i.e.normal) distribution with zero mean and unit standard deviation.
' Thus, le density is  f(x) = exp(-0.5*x**2)/sqrt(2.0*pi). le first call
' initializes grand and returns zero. The paramètre n is dummy.
' gRand calls a function rebord, and il is assumed that successive calls to rebord(0)
' give independent pseudo-random numbers distributed uniformly on (0,1), possibly
' including 0 (but not 1). le method used quoi suggested by de neumann, and
' improved by forsythe, ahrens, dieter and brent.
' on le average there sont 1.37746 calls of rebord for each call of grand.
' Warning - dimension and data statements below sont machine-dependent.
' Dimension of d must être at least le number of bits dans le fraction of a
' floating-point number. Thus, on most machines le data statement below
' can être truncated.
' si le integral of sqrt(2.0/pi)*exp(-0.5*x**2) à partir de
' 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![&Boucle]=val(d$[&Boucle])'::imprimer & Loop,format$("%g",d![& Loop])

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

    Whileloop 200000:cnt&=&Boucle

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

        si abs(Grnd!)<0.002

            locate 2,2:Imprimer 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 &Boucle,yh&-h&[&Boucle]

            endwhile

        endif

    Endwhile

    beep
    locate 2,2:Imprimer 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 si leading bit of u is one.
        g10:
        u!=u!+u!
        cas u!<1:goto "g20"
        u!=u!-1
        inc i&
        a!=a!-d![i&]
        goto "g10"
        ' forme w uniform on 0 < w < d(i+1) à partir de u.
        g20:
        w! = d![i&+1]*u!
        ' forme 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 si v! < u!
        cas v!<u!: goto "g40"
        ' generate random v.
        v! = rnd()
        ' loop si u .gt. v.
        cas u!>v!:goto "g30"
        ' reject w and forme a new uniform u à partir de v and u.
        u! = (v!-u!)/(1-u!)
        goto "g20"
        ' forme new u (to être used on next call) à partir de u and v.
        g40:
        u! = (u!-v!)/(1-v!)
        ' use first bit of u for sign, return normal variate.
        u!=u!+u!
        cas u!<1:goto "g50"
        u! = u! - 1
        grand! = w!-a!
        return grand!
        g50:
        grand! = a! - w!
        return grand!

    endproc