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
gamma_dist_knuth_algG(double a,double (* src)(unsigned short *),unsigned short * xi)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
gamma_dist_knuth_algA(double a,double (* src)(unsigned short *),unsigned short * xi)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
default_src(unsigned short * xi)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
gamma_dist_knuth(double a,double b)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
gamma_dist_knuth_src(double a,double b,double (* src)(unsigned short *),unsigned short * xi)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 }
143