Home | History | Annotate | Download | only in common
      1  5184  ek110237 /*
      2  5184  ek110237  * CDDL HEADER START
      3  5184  ek110237  *
      4  5184  ek110237  * The contents of this file are subject to the terms of the
      5  5184  ek110237  * Common Development and Distribution License (the "License").
      6  5184  ek110237  * You may not use this file except in compliance with the License.
      7  5184  ek110237  *
      8  5184  ek110237  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
      9  5184  ek110237  * or http://www.opensolaris.org/os/licensing.
     10  5184  ek110237  * See the License for the specific language governing permissions
     11  5184  ek110237  * and limitations under the License.
     12  5184  ek110237  *
     13  5184  ek110237  * When distributing Covered Code, include this CDDL HEADER in each
     14  5184  ek110237  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
     15  5184  ek110237  * If applicable, add the following below this CDDL HEADER, with the
     16  5184  ek110237  * fields enclosed by brackets "[]" replaced with your own identifying
     17  5184  ek110237  * information: Portions Copyright [yyyy] [name of copyright owner]
     18  5184  ek110237  *
     19  5184  ek110237  * CDDL HEADER END
     20  5184  ek110237  */
     21  5184  ek110237 /*
     22  6212  aw148015  * Copyright 2008 Sun Microsystems, Inc.  All rights reserved.
     23  5184  ek110237  * Use is subject to license terms.
     24  5184  ek110237  */
     25  5184  ek110237 
     26  5184  ek110237 #pragma ident	"%Z%%M%	%I%	%E% SMI"
     27  5184  ek110237 
     28  5184  ek110237 #include <stdlib.h>
     29  5184  ek110237 #include <math.h>
     30  5184  ek110237 
     31  5184  ek110237 /*
     32  5184  ek110237  * This is valid for 0 < a <= 1
     33  5184  ek110237  *
     34  5184  ek110237  * From Knuth Volume 2, 3rd edition, pages 586 - 587.
     35  5184  ek110237  */
     36  5184  ek110237 static double
     37  6212  aw148015 gamma_dist_knuth_algG(double a, double (*src)(unsigned short *),
     38  6212  aw148015     unsigned short *xi)
     39  5184  ek110237 {
     40  5184  ek110237 	double p, U, V, X, q;
     41  5184  ek110237 
     42  5184  ek110237 	p = M_E/(a + M_E);
     43  5184  ek110237 G2:
     44  6212  aw148015 	/* get a random number U */
     45  6212  aw148015 	U = (*src)(xi);
     46  6212  aw148015 
     47  5184  ek110237 	do {
     48  6212  aw148015 		/* get a random number V */
     49  6212  aw148015 		V = (*src)(xi);
     50  6212  aw148015 
     51  5184  ek110237 	} while (V == 0);
     52  5184  ek110237 
     53  5184  ek110237 	if (U < p) {
     54  5184  ek110237 		X = pow(V, 1/a);
     55  5184  ek110237 		/* q = e^(-X) */
     56  5184  ek110237 		q = exp(-X);
     57  5184  ek110237 	} else {
     58  5184  ek110237 		X = 1 - log(V);
     59  5184  ek110237 		q = pow(X, a-1);
     60  5184  ek110237 	}
     61  5184  ek110237 
     62  5184  ek110237 	/*
     63  5184  ek110237 	 * X now has density g, and q = f(X)/cg(X)
     64  5184  ek110237 	 */
     65  6212  aw148015 
     66  6212  aw148015 	/* get a random number U */
     67  6212  aw148015 	U = (*src)(xi);
     68  6212  aw148015 
     69  5184  ek110237 	if (U >= q)
     70  5184  ek110237 		goto G2;
     71  5184  ek110237 	return (X);
     72  5184  ek110237 }
     73  5184  ek110237 
     74  5184  ek110237 /*
     75  5184  ek110237  * This is valid for a > 1
     76  5184  ek110237  *
     77  5184  ek110237  * From Knuth Volume 2, 3rd edition, page 134.
     78  5184  ek110237  */
     79  5184  ek110237 static double
     80  6212  aw148015 gamma_dist_knuth_algA(double a, double (*src)(unsigned short *),
     81  6212  aw148015     unsigned short *xi)
     82  5184  ek110237 {
     83  5184  ek110237 	double U, Y, X, V;
     84  5184  ek110237 
     85  5184  ek110237 A1:
     86  6212  aw148015 	/* get a random number U */
     87  6212  aw148015 	U = (*src)(xi);
     88  6212  aw148015 
     89  5184  ek110237 	Y = tan(M_PI*U);
     90  5184  ek110237 	X = (sqrt((2*a) - 1) * Y) + a - 1;
     91  5184  ek110237 
     92  5184  ek110237 	if (X <= 0)
     93  5184  ek110237 		goto A1;
     94  5184  ek110237 
     95  6212  aw148015 	/* get a random number V */
     96  6212  aw148015 	V = (*src)(xi);
     97  6212  aw148015 
     98  5184  ek110237 	if (V > ((1 + (Y*Y)) * exp((a-1) * log(X/(a-1)) - sqrt(2*a -1) * Y)))
     99  5184  ek110237 		goto A1;
    100  5184  ek110237 
    101  5184  ek110237 	return (X);
    102  5184  ek110237 }
    103  5184  ek110237 
    104  6212  aw148015 /*
    105  6212  aw148015  * fetch a uniformly distributed random number using the drand48 generator
    106  6212  aw148015  */
    107  6212  aw148015 /* ARGSUSED */
    108  6212  aw148015 static double
    109  6212  aw148015 default_src(unsigned short *xi)
    110  6212  aw148015 {
    111  6212  aw148015 	return (drand48());
    112  6212  aw148015 }
    113  6212  aw148015 
    114  6212  aw148015 /*
    115  6212  aw148015  * Sample the gamma distributed random variable with gamma 'a' and
    116  6212  aw148015  * result mulitplier 'b', which is usually mean/gamma. Uses the default
    117  6212  aw148015  * drand48 random number generator as input
    118  6212  aw148015  */
    119  5184  ek110237 double
    120  5184  ek110237 gamma_dist_knuth(double a, double b)
    121  5184  ek110237 {
    122  5184  ek110237 	if (a <= 1.0)
    123  6212  aw148015 		return (b * gamma_dist_knuth_algG(a, default_src, NULL));
    124  5184  ek110237 	else
    125  6212  aw148015 		return (b * gamma_dist_knuth_algA(a, default_src, NULL));
    126  5184  ek110237 }
    127  6212  aw148015 
    128  6212  aw148015 /*
    129  6212  aw148015  * Sample the gamma distributed random variable with gamma 'a' and
    130  6212  aw148015  * multiplier 'b', which is mean / gamma adjusted for the specified
    131  6212  aw148015  * minimum value. The suppled random number source function is
    132  6212  aw148015  * used to optain the uniformly distributed random numbers.
    133  6212  aw148015  */
    134  6212  aw148015 double
    135  6212  aw148015 gamma_dist_knuth_src(double a, double b,
    136  6212  aw148015     double (*src)(unsigned short *), unsigned short *xi)
    137  6212  aw148015 {
    138  6212  aw148015 	if (a <= 1.0)
    139  6212  aw148015 		return (b * gamma_dist_knuth_algG(a, src, xi));
    140  6212  aw148015 	else
    141  6212  aw148015 		return (b * gamma_dist_knuth_algA(a, src, xi));
    142  6212  aw148015 }
    143