#include "k.h"
#define NUM_VRTS 2 /* no. of gamma variates */
/* 
*	Generate data from a compound multinomial (beta/binomial)
*   p is generated as Dirchlet with baseline prob 1/2,1/2.
*
*   Input m - number of alleles, theta 
* 	 *idum - pointer to random number work field.
*
*	Output rslt - vector of length 3: counts of each genotype 
*/
void datnull(int m, double theta, long *idum, int *rslt) {

int num=NUM_VRTS,i;														

double alphadot,alpha[2*NUM_VRTS],c_1[2*NUM_VRTS]={1.,1.,1.,1.},p,		
		z[2*NUM_VRTS]; 			/* storage for gamma variates   */

alphadot=(1.-theta)/theta;			/* calc alphadot		*/

for (i=0;i<2*num;i++) alpha[i]=1./2.*alphadot;	/* calc alpha			*/

rgamma(&num,alpha,c_1,z,idum);			/* get two gamma variates	*/
p=z[0]/(z[0]+z[1]);

datmult(m,p,idum,rslt); 
return;
}