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 2007 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) 38 { 39 double p, U, V, X, q; 40 41 p = M_E/(a + M_E); 42 G2: 43 U = drand48(); 44 do { 45 V = drand48(); 46 } while (V == 0); 47 48 if (U < p) { 49 X = pow(V, 1/a); 50 /* q = e^(-X) */ 51 q = exp(-X); 52 } else { 53 X = 1 - log(V); 54 q = pow(X, a-1); 55 } 56 57 /* 58 * X now has density g, and q = f(X)/cg(X) 59 */ 60 U = drand48(); 61 if (U >= q) 62 goto G2; 63 return (X); 64 } 65 66 /* 67 * This is valid for a > 1 68 * 69 * From Knuth Volume 2, 3rd edition, page 134. 70 */ 71 static double 72 gamma_dist_knuth_algA(double a) 73 { 74 double U, Y, X, V; 75 76 A1: 77 U = drand48(); 78 Y = tan(M_PI*U); 79 X = (sqrt((2*a) - 1) * Y) + a - 1; 80 81 if (X <= 0) 82 goto A1; 83 84 V = drand48(); 85 /* V > (1 + Y^2) * exp((a - 1) * log(X / (a - 1)) - sqrt(2a - 1) * Y) */ 86 if (V > ((1 + (Y*Y)) * exp((a-1) * log(X/(a-1)) - sqrt(2*a -1) * Y))) 87 goto A1; 88 89 return (X); 90 } 91 92 double 93 gamma_dist_knuth(double a, double b) 94 { 95 if (a <= 1.0) 96 return (b * gamma_dist_knuth_algG(a)); 97 else 98 return (b * gamma_dist_knuth_algA(a)); 99 } 100