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