program new;
{$DEFINE E_HIDE_ALL}
{$I SimbaExt/SimbaExt.simba}

function __gammavar(alpha,beta:Double): Double;
var
  u1,u2,b,c,v,z,r,ainv:Double;
const
  MAGIC_CONST := 1.0 + Ln(4.5);
begin
  if alpha > 1.0 then
  begin
    ainv := Sqrt(2.0 * alpha - 1.0);
    b := alpha - Ln(4);
    c := alpha + ainv;
    repeat
      u1 := Random();
      if not (0.0000001 < u1) and (u1 < 0.9999999) then
        continue;
      u2 := 1.0 - Random();
      v := Ln(u1 / (1.0 - u1)) / ainv;
      Result := alpha*Exp(v);
      z := u1 * u1 * u2;
      r := b + c * v - Result;
    until (r + MAGIC_CONST - 4.5*z >= 0.0) or (r >= Ln(z));
    Exit(Result * beta);
  end
  else if alpha = 1.0 then
  begin
    u1 := Random();
    while u1 <= 0.0000001 do u1 := Random();
    Exit(-Ln(u1) * beta);
  end else
    RaiseException('Alpha between 0 and 1 is not supported');
end;

function Random(lo,hi:Double; density:Double=14):Double; overload;
var v:Double;
begin
  if density > 1 then
  begin
    v := __gammavar(density,1.0)
    if v = 0 then Exit(lo);
    Result := v / (v + __gammavar(density,1.0));
  end else
    Result := Random();
  Result := Result / 1.0 * (hi-lo) + lo;
end;

function RandomPoint(mid:TPoint; radius:Int32; density:Double=14): TPoint;
begin
  Result.x := Round( Random(mid.x-radius,mid.x+radius, density) );
  Result.y := Round( Random(mid.y-radius,mid.y+radius, density) );
end;


//plot
var
  i:Int32;
  BMP:TRafBitmap;
  pt:TPoint;
begin
  BMP.Create(500,500);

  for i:=0 to 100000 do
  begin
    pt := RandomPoint([250,250],249, 5);  //change the last param (density)
    BMP.Pixel(pt.x,pt.y, 255);
  end;
  BMP.Debug();
end.