#include <math.h>
#define PI 3.141592654

float poidev(float xm, int *idum)
{
	float gammln(float xx);
	float ran1(int *idum);
	static float sq,alxm,g,oldm=(-1.0);
	float em,t,y;

	if (xm < 12.0) {
		if (xm != oldm) {
			oldm=xm;
			g=exp(-xm);
		}
		em = -1;
		t=1.0;
		do {
			++em;
			t *= ran1(idum);
		} while (t > g);
	} else {
		if (xm != oldm) {
			oldm=xm;
			sq=sqrt(2.0*xm);
			alxm=log(xm);
			g=xm*alxm-gammln(xm+1.0);
		}
		do {
			do {
				y=tan(PI*ran1(idum));
				em=sq*y+xm;
			} while (em < 0.0);
			em=floor(em);
			t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
		} while (ran1(idum) > t);
	}
	return em;
}

float poidev_copy(float xm, int *idum)
{
	float gammln(float xx);
	float ran1_copy(int *idum);
	static float sq1,alxm1,g1,oldm1=(-1.0);
	float em,t,y;

	if (xm < 12.0) {
		if (xm != oldm1) {
			oldm1=xm;
			g1=exp(-xm);
		}
		em = -1;
		t=1.0;
		do {
			++em;
			t *= ran1_copy(idum);
		} while (t > g1);
	} else {
		if (xm != oldm1) {
			oldm1=xm;
			sq1=sqrt(2.0*xm);
			alxm1=log(xm);
			g1=xm*alxm1-gammln(xm+1.0);
		}
		do {
			do {
				y=tan(PI*ran1_copy(idum));
				em=sq1*y+xm;
			} while (em < 0.0);
			em=floor(em);
			t=0.9*(1.0+y*y)*exp(em*alxm1-gammln(em+1.0)-g1);
		} while (ran1_copy(idum) > t);
	}
	return em;
}

#undef PI
/* (C) Copr. 1986-92 Numerical Recipes Software )1!. */