[Back to MATH SWAG index] [Back to Main SWAG index] [Original]

    {
    >Now I need a fast way of generating random numbers in Gaussian Distribution

    This is my code for gaussian and uniform distribution.
    After the unit there is a graphic test program.

    From: randyd@alpha2.csd.uwm.edu (Randall Elton Ding)
    }

    unit rndgauss;

    interface

    function rnd: double;                  { returns uniform [0..1] }

    function gauss(a,d: double): double;   { a is mean, d is std deviation }

    implementation

    function rnd: double;
      const
        bias = 1023;

      var
        data: record
                b: byte;
                d: double;
              end;
        x: array[0..8] of byte absolute data;
        e,i,j: word;

      begin
        for i:= 0 to 7 do x[i]:= lo(random(256));
        e:= bias;
        repeat
          j:= 0;
          for i:= 0 to 7 do begin
            j:= (x[i] shl 1) + hi(j);
            x[i]:= lo(j);
          end;
          e:= e-1;
          if (bias-e) mod 8 = 0 then x[0]:= lo(random(256));
        until (x[7] and $10) = $10;
        x[7]:= (x[7] and $0F) or lo(e shl 4);
        x[8]:= lo(e shr 4);
        rnd:= data.d;
      end;

    function gauss(a,d: double): double;
      const
        t: double = 0;

      var
        v1,v2,r: double;

      begin
        if t=0 then begin
          repeat
            v1:= 2*rnd-1;
            v2:= 2*rnd-1;
            r:= v1*v1+v2*v2
          until r<1;
          r:= sqrt((-2*ln(r))/r);
          t:= v2*r;
          gauss:= a+v1*r*d;
        end
        else begin
          gauss:= a+t*d;
          t:= 0;
        end;
      end;

    begin
    end.

    {---------------------- cut ---------------------------------}

    program gaussiantest;

    uses crt,graph,rndgauss;

    const
      bgipath = 'c:\bp\bgi';
      largestx = 999;

    procedure testplot;
      var
        htarry: array[0..largestx] of integer;
        x,y,w,h,m,v: word;

      begin
        fillchar(htarry,sizeof(htarry),#0);
        w:= getmaxx+1;
        h:= getmaxy;
        m:= getmaxx div 2;
        v:= getmaxx div 8;
        while not keypressed do begin
          x:= trunc(gauss(m,v));
          if x<=largestx then begin
            y:= htarry[x];
            if y<=h then begin
              putpixel(x,h-y,white);
              htarry[x]:= y+1;
            end;
          end;
        end;
      end;

    procedure initbgi;
      var errcode,grmode,grdriver: integer;
      begin
        grdriver:= detect;
        initgraph (grdriver,grmode,bgipath);
        errcode:= graphresult;
        if errcode <> grok then begin
          writeln ('Graphics error: ',grapherrormsg (errcode));
          halt (1);
        end;
      end;

    begin
      initbgi;
      testplot;
      closegraph;
    end.

[Back to MATH SWAG index] [Back to Main SWAG index] [Original]