/*
    gamma.c
    $Id: gamma.c,v 1.3 2018/05/08 04:57:59 daichi Exp $

*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define RANDOM ((double)rand()/(double)RAND_MAX)  /* 実際にはMTなどを推奨 */

double
exprand (void)
{
	return (- log(RANDOM));
}

double
sgamrand (double a)
{
	double x, y, z;
	double u, v, w, b, c, e;
	int accept = 0;
	if (a <= 1)
	{
		/* a <= 1. Johnk's generator. Devroye (1986) p.418 */
		e = exprand();
		do {
			x = pow(RANDOM, 1 / a);
			y = pow(RANDOM, 1 / (1 - a));
		} while (x + y > 1);
		return (e * x / (x + y));
	} else {
		/* a > 1. Best's rejection algorithm. Devroye (1986) p.410 */
		b = a - 1;
		c = 3 * a - 0.75;
		while (accept != 1)
		{
			/* generate */
			u = RANDOM;
			v = RANDOM;
			w = u * (1 - u);
			y = sqrt(c / w) * (u - 0.5);
			x = b + y;
			if (x >= 0)
			{
				z = 64 * w * w * w * v * v;
				if (z <= 1 - (2 * y * y) / x)
				{
					accept = 1;
				} else {
					if (log(z) <= 2 * (b * log(x / b) - y))
						accept = 1;
				}
			}
		}
		return x;
	}
}

double
gamrand (double a, double b)
{
	return sgamrand(a) / b;
}

#if 1

int
main (int argc, char *argv[])
{
	if (argc < 4) {
		printf ("usage: %% gamma a b n\n");
		exit (0);
	}
	double a = atof(argv[1]), b = atof(argv[2]);
	int i, n = atoi(argv[3]);

	for (i = 0; i < n; i++)
		printf("%lf\n", gamrand (a, b));
	

}
#endif