1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 /*
  22  * Copyright 2008 Sun Microsystems, Inc.  All rights reserved.
  23  * Use is subject to license terms.
  24  */
  25 
  26 #pragma ident   "%Z%%M% %I%     %E% SMI"
  27 
  28 #include <stdlib.h>
  29 #include <math.h>
  30 
  31 /*
  32  * This is valid for 0 < a <= 1
  33  *
  34  * From Knuth Volume 2, 3rd edition, pages 586 - 587.
  35  */
  36 static double
  37 gamma_dist_knuth_algG(double a, double (*src)(unsigned short *),
  38     unsigned short *xi)
  39 {
  40         double p, U, V, X, q;
  41 
  42         p = M_E/(a + M_E);
  43 G2:
  44         /* get a random number U */
  45         U = (*src)(xi);
  46 
  47         do {
  48                 /* get a random number V */
  49                 V = (*src)(xi);
  50 
  51         } while (V == 0);
  52 
  53         if (U < p) {
  54                 X = pow(V, 1/a);
  55                 /* q = e^(-X) */
  56                 q = exp(-X);
  57         } else {
  58                 X = 1 - log(V);
  59                 q = pow(X, a-1);
  60         }
  61 
  62         /*
  63          * X now has density g, and q = f(X)/cg(X)
  64          */
  65 
  66         /* get a random number U */
  67         U = (*src)(xi);
  68 
  69         if (U >= q)
  70                 goto G2;
  71         return (X);
  72 }
  73 
  74 /*
  75  * This is valid for a > 1
  76  *
  77  * From Knuth Volume 2, 3rd edition, page 134.
  78  */
  79 static double
  80 gamma_dist_knuth_algA(double a, double (*src)(unsigned short *),
  81     unsigned short *xi)
  82 {
  83         double U, Y, X, V;
  84 
  85 A1:
  86         /* get a random number U */
  87         U = (*src)(xi);
  88 
  89         Y = tan(M_PI*U);
  90         X = (sqrt((2*a) - 1) * Y) + a - 1;
  91 
  92         if (X <= 0)
  93                 goto A1;
  94 
  95         /* get a random number V */
  96         V = (*src)(xi);
  97 
  98         if (V > ((1 + (Y*Y)) * exp((a-1) * log(X/(a-1)) - sqrt(2*a -1) * Y)))
  99                 goto A1;
 100 
 101         return (X);
 102 }
 103 
 104 /*
 105  * fetch a uniformly distributed random number using the drand48 generator
 106  */
 107 /* ARGSUSED */
 108 static double
 109 default_src(unsigned short *xi)
 110 {
 111         return (drand48());
 112 }
 113 
 114 /*
 115  * Sample the gamma distributed random variable with gamma 'a' and
 116  * result mulitplier 'b', which is usually mean/gamma. Uses the default
 117  * drand48 random number generator as input
 118  */
 119 double
 120 gamma_dist_knuth(double a, double b)
 121 {
 122         if (a <= 1.0)
 123                 return (b * gamma_dist_knuth_algG(a, default_src, NULL));
 124         else
 125                 return (b * gamma_dist_knuth_algA(a, default_src, NULL));
 126 }
 127 
 128 /*
 129  * Sample the gamma distributed random variable with gamma 'a' and
 130  * multiplier 'b', which is mean / gamma adjusted for the specified
 131  * minimum value. The suppled random number source function is
 132  * used to optain the uniformly distributed random numbers.
 133  */
 134 double
 135 gamma_dist_knuth_src(double a, double b,
 136     double (*src)(unsigned short *), unsigned short *xi)
 137 {
 138         if (a <= 1.0)
 139                 return (b * gamma_dist_knuth_algG(a, src, xi));
 140         else
 141                 return (b * gamma_dist_knuth_algA(a, src, xi));
 142 }