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 }