Home | History | Annotate | Download | only in common
      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 }
    143