/* Perfect Gibbs simulation of auto-Gamma: ONLY COALESCENCE TIMES!! */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define eps 0

double vgamma(double alpha)
{    
    const double e = 2.71828182845905;
    double b,p,u0,u1,x,a,u,ustar,aa,c1,c2,c3,c4,c5,w,u2;
    
    if (alpha < 1.0){
	b = (alpha + e)/e;
    a10:
	p = b*drand48();
	if (p >= 1.0){
	    printf("Out of bounds!\n");
	    goto a20;
	}
	x = exp(log(p)/alpha);
	if (x >= -1.0*log(drand48()))
	    goto a10;
	return(x);
    a20:
	x = -1.0*log((b-p)/alpha);
	if (exp(log(x)*(alpha-1.0)) <= drand48())
	    goto a10;
	return(x);
    }
    
    if (alpha == 1.0)
	{
	    a = 0.0;
	b10:
	    u = drand48();
	    u0 = u;
	b20:
	    ustar = drand48();
	    if (u <= ustar)
		goto b30;
	    u = drand48();
	    if (u < ustar)
		goto b20;
	    a = a + 1.0;
	    goto b10;
	b30:
	    return(a+u0);
	}
    
    if (alpha > 1.0)
	{
	    c1 = alpha - 1.0;
	    aa = 1.0/c1;
	    c2 = aa*(alpha - 1.0/(6.0*alpha));
	    c3 = 2.0*aa;
	    c4 = c3 + 2.0;
	    if (alpha >  2.5)
		c5 = 1.0/sqrt(alpha);
	c10:
	    u1 = drand48();
	    u2 = drand48();
	    if (alpha <= 2.5)
		goto c20;
	    u1 = u2 + c5*(1.0 - 1.86*u1);
	    if (u1 <= 0.0 || u1 >= 1.0)
		goto c10;
	c20:
	    w = c2*u2/u1;
	    if (c3*u1 + w + 1/w <= c4)
		goto c30;
	    if (c3*log(u1) - log(w) + w >= 1.0)
		goto c10;
	c30:
	    x = c1*w;
	    return(x);
	} 
}

main()
{
    
    int     i,j,coal;
    long    t;
    double  G,U[11],L[11],sumL,sumU,newsumL,newsumU, mean,
    std,alpha[11],beta[11];
    
    /* Initialization */
    
    srand48(76252);
    
    mean = 0.000; 
    std = 0.000;
    
    alpha[0] = 18.030;
    alpha[1] =  6.802;
    alpha[2] =  2.802;
    alpha[3] =  6.802;
    alpha[4] = 15.802;
    alpha[5] =  4.802;
    alpha[6] = 20.802;
    alpha[7] =  2.802;
    alpha[8] =  2.802;
    alpha[9] =  5.802;
    alpha[10]= 23.802;
    
    beta[0] =   1.000;
    beta[1] =  94.320;
    beta[2] =  15.720;
    beta[3] =  62.880;
    beta[4] = 125.760;
    beta[5] =   5.240;
    beta[6] =  31.440;
    beta[7] =   1.048;
    beta[8] =   1.048;
    beta[9] =   2.096;
    beta[10] = 10.480; 
    
    
    /* Perfect simulation starts - forwards in time */
    
    for(j = 1; j <= 10000; j++)
	{ 
	    sumL = 0;
	    sumU = 0;
	    t = 0;
	    
	    for(i = 0; i<= 10; i++)
		{
		    L[i] = 0;
		    U[i] = vgamma(alpha[i]) / beta[i];
		    if (i > 0)
			sumU = sumU + U[i];
		}
	    do
		{
		    t = t + 1; 
		    
		    for(i = 0; i <= 10; i++)  
			      {
				  G = vgamma(alpha[i]);
				  if (i == 0)
				      {
					  U[0] = G / (beta[0] + sumL);
					  L[0] = G / (beta[0] + sumU);
				      }
				  else
				      {
					  U[i] = G / (beta[i] + L[0]);
					  L[i] = G / (beta[i] + U[0]);
					  
				      }
				  
				  if (i == 0)
				      {
					  newsumL = 0;
					  newsumU = 0;
				      }
				  else
				      {
					  newsumL = newsumL + L[i];
					  newsumU = newsumU + U[i];
				      }
			      }

		    sumL = newsumL;
		    sumU = newsumU;
		    
		    coal = 1;
		    i = 0;
		  
		    while (coal == 1 && i <= 10)
		      {
			  if (U[i]-L[i] > eps)
			      coal = -1;
			  i++;
		      }
	    
		}
	    while(coal == -1 && t <= 1000000000);
	    
	    if (t >= 1000000000)
		printf("Too small time for generating forwards in time!!!!");
	    mean = mean + t;
	    std = std + t * t;
	    
      }

    mean = mean/10000;
    std = sqrt((std-(10000*mean*mean))/9999) / 100;
    printf("Coalescence time: mean = %.5f, std = %.5f, epsilon = %.14f\n", mean, std, eps);
}