xref: /onnv-gate/usr/src/cmd/filebench/common/gamma_dist.c (revision 5184:da60d2b4a9e2)
1*5184Sek110237 /*
2*5184Sek110237  * CDDL HEADER START
3*5184Sek110237  *
4*5184Sek110237  * The contents of this file are subject to the terms of the
5*5184Sek110237  * Common Development and Distribution License (the "License").
6*5184Sek110237  * You may not use this file except in compliance with the License.
7*5184Sek110237  *
8*5184Sek110237  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9*5184Sek110237  * or http://www.opensolaris.org/os/licensing.
10*5184Sek110237  * See the License for the specific language governing permissions
11*5184Sek110237  * and limitations under the License.
12*5184Sek110237  *
13*5184Sek110237  * When distributing Covered Code, include this CDDL HEADER in each
14*5184Sek110237  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15*5184Sek110237  * If applicable, add the following below this CDDL HEADER, with the
16*5184Sek110237  * fields enclosed by brackets "[]" replaced with your own identifying
17*5184Sek110237  * information: Portions Copyright [yyyy] [name of copyright owner]
18*5184Sek110237  *
19*5184Sek110237  * CDDL HEADER END
20*5184Sek110237  */
21*5184Sek110237 /*
22*5184Sek110237  * Copyright 2007 Sun Microsystems, Inc.  All rights reserved.
23*5184Sek110237  * Use is subject to license terms.
24*5184Sek110237  */
25*5184Sek110237 
26*5184Sek110237 #pragma ident	"%Z%%M%	%I%	%E% SMI"
27*5184Sek110237 
28*5184Sek110237 #include <stdlib.h>
29*5184Sek110237 #include <math.h>
30*5184Sek110237 
31*5184Sek110237 /*
32*5184Sek110237  * This is valid for 0 < a <= 1
33*5184Sek110237  *
34*5184Sek110237  * From Knuth Volume 2, 3rd edition, pages 586 - 587.
35*5184Sek110237  */
36*5184Sek110237 static double
37*5184Sek110237 gamma_dist_knuth_algG(double a)
38*5184Sek110237 {
39*5184Sek110237 	double p, U, V, X, q;
40*5184Sek110237 
41*5184Sek110237 	p = M_E/(a + M_E);
42*5184Sek110237 G2:
43*5184Sek110237 	U = drand48();
44*5184Sek110237 	do {
45*5184Sek110237 		V = drand48();
46*5184Sek110237 	} while (V == 0);
47*5184Sek110237 
48*5184Sek110237 	if (U < p) {
49*5184Sek110237 		X = pow(V, 1/a);
50*5184Sek110237 		/* q = e^(-X) */
51*5184Sek110237 		q = exp(-X);
52*5184Sek110237 	} else {
53*5184Sek110237 		X = 1 - log(V);
54*5184Sek110237 		q = pow(X, a-1);
55*5184Sek110237 	}
56*5184Sek110237 
57*5184Sek110237 	/*
58*5184Sek110237 	 * X now has density g, and q = f(X)/cg(X)
59*5184Sek110237 	 */
60*5184Sek110237 	U = drand48();
61*5184Sek110237 	if (U >= q)
62*5184Sek110237 		goto G2;
63*5184Sek110237 	return (X);
64*5184Sek110237 }
65*5184Sek110237 
66*5184Sek110237 /*
67*5184Sek110237  * This is valid for a > 1
68*5184Sek110237  *
69*5184Sek110237  * From Knuth Volume 2, 3rd edition, page 134.
70*5184Sek110237  */
71*5184Sek110237 static double
72*5184Sek110237 gamma_dist_knuth_algA(double a)
73*5184Sek110237 {
74*5184Sek110237 	double U, Y, X, V;
75*5184Sek110237 
76*5184Sek110237 A1:
77*5184Sek110237 	U = drand48();
78*5184Sek110237 	Y = tan(M_PI*U);
79*5184Sek110237 	X = (sqrt((2*a) - 1) * Y) + a - 1;
80*5184Sek110237 
81*5184Sek110237 	if (X <= 0)
82*5184Sek110237 		goto A1;
83*5184Sek110237 
84*5184Sek110237 	V = drand48();
85*5184Sek110237 	/* V > (1 + Y^2) * exp((a - 1) * log(X / (a - 1)) - sqrt(2a - 1) * Y) */
86*5184Sek110237 	if (V > ((1 + (Y*Y)) * exp((a-1) * log(X/(a-1)) - sqrt(2*a -1) * Y)))
87*5184Sek110237 		goto A1;
88*5184Sek110237 
89*5184Sek110237 	return (X);
90*5184Sek110237 }
91*5184Sek110237 
92*5184Sek110237 double
93*5184Sek110237 gamma_dist_knuth(double a, double b)
94*5184Sek110237 {
95*5184Sek110237 	if (a <= 1.0)
96*5184Sek110237 		return (b * gamma_dist_knuth_algG(a));
97*5184Sek110237 	else
98*5184Sek110237 		return (b * gamma_dist_knuth_algA(a));
99*5184Sek110237 }
100