Contributor: RANDALL ELTON DING

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

>I a program I'm currently struggeling with, I need to get a random number
>from a Gaussian distribution. Anybody got any ideas or anybody able to point
>to something which does the job.

This does a pretty good job of generating a gaussian random variable
with mean `a` and standard deviation `d`.
This program also does a graphic plot to demonstrate the function.

First, here is the origional C source if the gaussian function
which I transcribed to beloved pascal..

/* ------------------------------------------------ *
 * gaussian -- generates a gaussian random variable *
 *             with mean a and standard deviation d *
 * ------------------------------------------------ */
 double gaussian(a,d)
 double a,d;
 {
   static double t = 0.0;
   double x,v1,v2,r;
   if (t == 0) {
     do {
       v1 = 2.0 * rnd() - 1.0;
       v2 = 2.0 * rnd() - 1.0;
       r = v1 * v1 + v2 * v2;
     } while (r>=1.0);
     r = sqrt((-2.0*log(r))/r);
     t = v2*r;
     return(a+v1*r*d);
   }
   else {
     x = t;
     t = 0.0;
     return(a+x*d);
   }
 }


* ----------------------------------------------------------------------
* now, the same thing in pascal
* ----------------------------------------------------------------------
*)

{$N+}
program testgaussian;

uses graph,crt;

const
  bgipath = 'e:\bp\bgi';

procedure initbgi;
  var
    errcode,grdriver,grmode: integer;

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



function rnd: double;   { this isn't the best, but it works }
  var                   { returns a random number between 0 and 1 }
    i: integer;
    r: double;

  begin
    r:= 0;
    for i:= 1 to 15 do begin
      r:= r + random(10);
      r:= r/10;
    end;
    rnd:= r;
  end;



function gaussian(a,d: double): double;      { a is mean }
  const                                      { d is standard deviation }
    t: double = 0;   { pascal's equivalent to C's static variable }

  var
    x,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;
      gaussian:= a+v1*r*d;
    end
    else begin
      x:= t;
      t:= 0;
      gaussian:= a+x*d;
    end;
  end;



procedure testplot;
  var
    x,mx,my,y1: word;
    y: array[1..999] of word;
              { ^^^ make this bigger if you have incredible graphics }
  begin
    initbgi;
    mx:= getmaxx+1;
    my:= getmaxy;
    fillchar(y,sizeof(y),#0);
    repeat
      x:= trunc(gaussian(mx/2,50));
      y1:= y[x];
      putpixel(x,my-y1,white);
      y[x]:= y1+1;
    until keypressed;
    closegraph;
  end;



begin
  randomize;
  testplot;
end.