public static double nextRandom(double k,double theta,Random random){
  final double q1=0.0416666664, q2=0.0208333723, q3=0.0079849875;
  final double q4=0.0015746717, q5=-0.0003349403, q6=0.0003340332;
  final double q7=0.0006053049, q8=-0.0004701849, q9=0.0001710320;
  final double a1=0.333333333, a2=-0.249999949, a3=0.199999867;
  final double a4=-0.166677482, a5=0.142873973, a6=-0.124385581;
  final double a7=0.110368310, a8=-0.112750886, a9=0.104089866;
  final double e1=1.000000000, e2=0.499999994, e3=0.166666848;
  final double e4=0.041664508, e5=0.008345522, e6=0.001353826;
  final double e7=0.000247453;
  if (k < 1.0) {
    final double b=1.0 + 0.36788794412 * k;
    while (true) {
      final double p=b * random.nextDouble();
      if (p <= 1.0) {
        final double gds=Math.exp(Math.log(p) / k);
        if (Math.log(random.nextDouble()) <= -gds) {
          return (gds / theta);
        }
      }
 else {
        final double gds=-Math.log((b - p) / k);
        if (Math.log(random.nextDouble()) <= ((k - 1.0) * Math.log(gds))) {
          return (gds / theta);
        }
      }
    }
  }
 else {
    final double ss, s, d;
    if (k != -1.0) {
      ss=k - 0.5;
      s=Math.sqrt(ss);
      d=5.656854249 - 12.0 * s;
    }
 else {
      ss=0.0;
      s=0.0;
      d=0.0;
    }
    final double v1, v12;
{
      double tv1, tv2, tv12;
      do {
        tv1=2.0 * random.nextDouble() - 1.0;
        tv2=2.0 * random.nextDouble() - 1.0;
        tv12=tv1 * tv1 + tv2 * tv2;
      }
 while (tv12 > 1.0);
      v1=tv1;
      v12=tv12;
    }
    final double b, c, si, q0;
{
      final double t=v1 * Math.sqrt(-2.0 * Math.log(v12) / v12);
      final double x=s + 0.5 * t;
      final double gds=x * x;
      if (t >= 0.0) {
        return (gds / theta);
      }
      final double un=random.nextDouble();
      if (d * un <= t * t * t) {
        return (gds / theta);
      }
      if (k != -1.0) {
        final double r=1.0 / k;
        q0=((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) * r + q3) * r + q2) * r + q1) * r;
        if (k > 3.686) {
          if (k > 13.022) {
            b=1.77;
            si=0.75;
            c=0.1515 / s;
          }
 else {
            b=1.654 + 0.0076 * ss;
            si=1.68 / s + 0.275;
            c=0.062 / s + 0.024;
          }
        }
 else {
          b=0.463 + s - 0.178 * ss;
          si=1.235;
          c=0.195 / s - 0.079 + 0.016 * s;
        }
      }
 else {
        b=0.0;
        c=0.0;
        si=0.0;
        q0=0.0;
      }
      if (x > 0.0) {
        final double v=t / (s + s);
        final double q;
        if (Math.abs(v) > 0.25) {
          q=q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
        }
 else {
          q=q0 + 0.5 * t * t* ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1)* v;
        }
        if (Math.log(1.0 - un) <= q) {
          return (gds / theta);
        }
      }
    }
    while (true) {
      double e, u, sign_u, t;
      do {
        e=-Math.log(random.nextDouble());
        u=random.nextDouble();
        u=u + u - 1.0;
        sign_u=(u > 0) ? 1.0 : -1.0;
        t=b + (e * si) * sign_u;
      }
 while (t <= -0.71874483771719);
      final double v=t / (s + s);
      final double q;
      if (Math.abs(v) > 0.25) {
        q=q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
      }
 else {
        q=q0 + 0.5 * t * t* ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1)* v;
      }
      if (q <= 0.0) {
        continue;
      }
      final double w;
      if (q > 0.5) {
        w=Math.exp(q) - 1.0;
      }
 else {
        w=((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) * q + e1) * q;
      }
      if (c * u * sign_u <= w * Math.exp(e - 0.5 * t * t)) {
        final double x=s + 0.5 * t;
        return (x * x / theta);
      }
    }
  }
}

public static double nextRandom ( double k , double theta , Random random ) { final double q1 = 0.0416666664 , q2 = 0.0208333723 , q3 = 0.0079849875 ; final double q4 = 0.0015746717 , q5 = - 0.0003349403 , q6 = 0.0003340332 ; final double q7 = 0.0006053049 , q8 = - 0.0004701849 , q9 = 0.0001710320 ; final double a1 = 0.333333333 , a2 = - 0.249999949 , a3 = 0.199999867 ; final double a4 = - 0.166677482 , a5 = 0.142873973 , a6 = - 0.124385581 ; final double a7 = 0.110368310 , a8 = - 0.112750886 , a9 = 0.104089866 ; final double e1 = 1.000000000 , e2 = 0.499999994 , e3 = 0.166666848 ; final double e4 = 0.041664508 , e5 = 0.008345522 , e6 = 0.001353826 ; final double e7 = 0.000247453 ; if ( k < 1.0 ) { final double b = 1.0 + 0.36788794412 * k ; while ( true ) { final double p = b * random . nextDouble ( ) ; if ( p <= 1.0 ) { final double gds = Math . exp ( Math . log ( p ) / k ) ; if ( Math . log ( random . nextDouble ( ) ) <= - gds ) { return ( gds / theta ) ; } } else { final double gds = - Math . log ( ( b - p ) / k ) ; if ( Math . log ( random . nextDouble ( ) ) <= ( ( k - 1.0 ) * Math . log ( gds ) ) ) { return ( gds / theta ) ; } } } } else { final double ss , s , d ; if ( k != - 1.0 ) { ss = k - 0.5 ; s = Math . sqrt ( ss ) ; d = 5.656854249 - 12.0 * s ; } else { ss = 0.0 ; s = 0.0 ; d = 0.0 ; } final double v1 , v12 ; { double tv1 , tv2 , tv12 ; do { tv1 = 2.0 * random . nextDouble ( ) - 1.0 ; tv2 = 2.0 * random . nextDouble ( ) - 1.0 ; tv12 = tv1 * tv1 + tv2 * tv2 ; } while ( tv12 > 1.0 ) ; v1 = tv1 ; v12 = tv12 ; } final double b , c , si , q0 ; { final double t = v1 * Math . sqrt ( - 2.0 * Math . log ( v12 ) / v12 ) ; final double x = s + 0.5 * t ; final double gds = x * x ; if ( t >= 0.0 ) { return ( gds / theta ) ; } final double un = random . nextDouble ( ) ; if ( d * un <= t * t * t ) { return ( gds / theta ) ; } if ( k != - 1.0 ) { final double r = 1.0 / k ; q0 = ( ( ( ( ( ( ( ( q9 * r + q8 ) * r + q7 ) * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r + q2 ) * r + q1 ) * r ; if ( k > 3.686 ) { if ( k > 13.022 ) { b = 1.77 ; si = 0.75 ; c = 0.1515 / s ; } else { b = 1.654 + 0.0076 * ss ; si = 1.68 / s + 0.275 ; c = 0.062 / s + 0.024 ; } } else { b = 0.463 + s - 0.178 * ss ; si = 1.235 ; c = 0.195 / s - 0.079 + 0.016 * s ; } } else { b = 0.0 ; c = 0.0 ; si = 0.0 ; q0 = 0.0 ; } if ( x > 0.0 ) { final double v = t / ( s + s ) ; final double q ; if ( Math . abs ( v ) > 0.25 ) { q = q0 - s * t + 0.25 * t * t + ( ss + ss ) * Math . log ( 1.0 + v ) ; } else { q = q0 + 0.5 * t * t * ( ( ( ( ( ( ( ( a9 * v + a8 ) * v + a7 ) * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v + a2 ) * v + a1 ) * v ; } if ( Math . log ( 1.0 - un ) <= q ) { return ( gds / theta ) ; } } } while ( true ) { double e , u , sign_u , t ; do { e = - Math . log ( random . nextDouble ( ) ) ; u = random . nextDouble ( ) ; u = u + u - 1.0 ; sign_u = ( u > 0 ) ? 1.0 : - 1.0 ; t = b + ( e * si ) * sign_u ; } while ( t <= - 0.71874483771719 ) ; final double v = t / ( s + s ) ; final double q ; if ( Math . abs ( v ) > 0.25 ) { q = q0 - s * t + 0.25 * t * t + ( ss + ss ) * Math . log ( 1.0 + v ) ; } else { q = q0 + 0.5 * t * t * ( ( ( ( ( ( ( ( a9 * v + a8 ) * v + a7 ) * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v + a2 ) * v + a1 ) * v ; } if ( q <= 0.0 ) { continue ; } final double w ; if ( q > 0.5 ) { w = Math . exp ( q ) - 1.0 ; } else { w = ( ( ( ( ( ( e7 * q + e6 ) * q + e5 ) * q + e4 ) * q + e3 ) * q + e2 ) * q + e1 ) * q ; } if ( c * u * sign_u <= w * Math . exp ( e - 0.5 * t * t ) ) { final double x = s + 0.5 * t ; return ( x * x / theta ) ; } } } }
generate a random value with the generators parameters . along the lines of - j . h . ahrens , u . dieter ( 1974 ) : computer methods for sampling from gamma , beta , poisson and binomial distributions , computing 12 , 223 - 246 . - j . h . ahrens , u . dieter ( 1982 ) : generating gamma variates by a modified rejection technique , communications of the acm 25 , 47 - 54 .
