vardef normaldeviate =
 save u, v, xa;
 forever:
   u := 1 - uniformdeviate 63/64;
   v := sqrt(8/mexp(256)) * (-1/2 + uniformdeviate 1);
   xa := v / u;
   exitif (xa * xa <= -mlog(u)/64);
 endfor
 xa
enddef;

vardef exponentialdeviate expr mu =
 save u;
 u = 1 - uniformdeviate 1;  % hence 0 < u <= 1 and so you avoid
 -mu * 1/256 mlog(u)        % the danger of calling mlog(0)
enddef;

vardef tand(expr theta) =
 save x, y;
 (x, y) = dir theta;
 if abs(x) < eps: infinity else: y/x fi
enddef;
vardef exp(expr x) = mexp(256x) enddef;
vardef log(expr x) = 1/256 mlog(x) enddef;

% this is defined only for a > 1
vardef gammadeviate(expr a, b) =
 save y, x, v, s, accept; boolean accept;
 s = sqrt(2a - 1);
 forever:
   forever:
     y := tand(uniformdeviate 180);
     exitif abs(y) < 16;
   endfor
   x := s * y + a - 1;
   accept := false;
   if x > 0:
     v := uniformdeviate 1;
     if v <= (1 + y * y) * exp((a - 1) * log(x / (a - 1)) - s * y):
       accept := true;
     fi
   fi
   exitif accept;
 endfor
 x/b
enddef;