/*   Alea.c:  

These are the most common generators I use in my C programs. 
Warning: No support nor warranty associated with these!

*/

#include <math.h>
#include <stdio.h>
#include <malloc.h>
#include "Avoga.h"

double kiss();
void nor();
double normal_cdf();
double normal_quantile();
double Gam_1();
double Gam_2();
double Gam_3();
double Gam();
double Beta();
double Trun_Gam();
double Trunc_norm_left();
double Trunc_norm();
int n_optimal();
double inter_ri();
double gamma_right();
double integer();
double inter_le();
double gamma_left ();


/*__________________Kiss, The generator____________________________*/

double kiss(i,j,k)
	unsigned long int *i,*j,*k;

/* Generator proposed by Marsaglia and Zaman, 1993. See
   Robert and Casella (1999, pages 41-43) for details.  
   Watch out: the last line
	x = ((double) (*i+*j+*k)*exp(-32*log(2.0)));
   must be calibrated, depending on the precision 
   of the computer. */
{
  double x;

  *j = *j ^ (*j<<17);
  *k = (*k ^ (*k<<18)) & 0x7FFFFFFF;
  *i = 69069*(*i)+23606797;
  *j ^= (*j>>15);
  *k ^= (*k>>13);
  x = ((double) (*i+*j+*k)*exp(-32*log(2.0)));
  return(x);
}

/*__________________Normal (Box-Muller)____________________________*/

void nor (t1, t2, i0, j0, k0)
	unsigned long int *i0,*j0,*k0;
	double   *t1,*t2;

/* The most standard normal generator;
   see Robert and Casella (1999, page 46, Algorithm A.4)
   for details. */
{
  double z,u2;
  unsigned long int i,j,k;

  i=*i0; j=*j0; k=*k0;
  u2=kiss(&i,&j,&k);
  if (!((u2>0.)&&(u2<1.)))
	u2=kiss(&i,&j,&k);
  z = sqrt( (-2.0*( log( u2 ) ) ) );
  u2 = 2*PI*kiss(&i,&j,&k);
  *i0=i; *j0=j; *k0=k;
  *t1 = (double) ( z*cos(u2) );
  *t2 = (double) ( z*sin(u2) );
}

/*__________________Gamma(a,1), a>1____________________________*/

double Gam_1 (alf, i0, j0, k0)
	double alf;
	unsigned long int *i0,*j0,*k0;

/* Best (1978) gamma generator. Only applies when alf > 1;
   see Robert and Casella (1999, page 55, Algorithm A.7)
   for details.  */
{
  double b,d,u1,u2,v,w;
  unsigned long int i,j,k;
  boolean stoop;
  
  i=*i0; j=*j0; k=*k0;
  b = alf -1;
  d = 3*alf - 0.75;
  stoop = FALSE;

  while (! stoop)
  {  u1 = kiss(&i,&j,&k);
     w = u1*(1-u1);
     v = b + sqrt(d/w)*(u1-0.5);
     stoop = (v>0) && (kiss(&i,&j,&k)<(exp(b-v+b*log(v/b))/sqrt(64*w*w*w)));
  }
  *i0=i; *j0=j; *k0=k;
  return(v);
}

/*__________________Gamma(a,1) a<1____________________________*/

double Gam_2(alf,i0,j0,k0)
	double alf;
	unsigned long int *i0,*j0,*k0;

/* Gamma generator when   alf < 1 ;
   see Robert and Casella (1999, page 47, Example 2.2.5) */
{
  double x;
  unsigned long int i,j,k;

  i=*i0; j=*j0; k=*k0;
  x =exp(log(kiss(&i,&j,&k))/alf)*Gam_1(alf+1,&i,&j,&k);
  *i0=i; *j0=j; *k0=k;
  return (x);
}

/*__________________Gamma(a,1), a>>1____________________________*/

double Gam_3 (alf, i0, j0, k0)
	double alf;
	unsigned long int *i0, *j0, *k0;

/* The standard normal approximation for very
   large values of  alf . */
{
  double x,y;

  nor(&x,&y,i0,j0,k0);
  return(x*sqrt(alf)+alf);
}

/*__________________Gamma(a,1), summary____________________________*/

double Gam (alf, i0, j0, k0)
	double alf;
	unsigned long int *i0, *j0, *k0;
{
  double x;
  unsigned long int i,j,k;

  i=*i0; j=*j0; k=*k0;
  if (alf>1.1)
  {  if (alf<20.)
	x = Gam_1(alf,&i,&j,&k);
     else
	x = Gam_3(alf,&i,&j,&k);
  }
  else
  { if (alf==1.)
	x = -log(kiss(&i,&j,&k));
    else
	x = Gam_2(alf,&i,&j,&k);
  }

  *i0 =i; *j0=j; *k0=k;
  return (x);
}

double Beta(alf,bet,i0,j0,k0)
	double alf,bet;
	unsigned long int *i0, *j0, *k0;
{
  double x,y;
  unsigned long int i,j,k;

  i=*i0; j=*j0; k=*k0; 
  x = Gam (alf,&i,&j,&k);
  y = Gam (bet,&i,&j,&k);
  *i0=i; *j0=j; *k0=k;
  return ( x/(x+y) );
}


/*__________________Truncated Gamma(a,b), x<1____________________________*/

double Trun_Gam(a,b,i0,j0,k0)
   double a,b;
   unsigned long int *i0, *j0, *k0;

/* Gamma generator when  x < 1. 
   Watch out for approximations, like those
   induced by the stopping rule  infloop<1000  */
{
   double x0, y, M, ga, un;
   int x, infloop;
   unsigned long int u, v, w;

   u=*i0; v=*j0; w=*k0;
   x0=(a-0.5)/b; infloop=0;

   if (x0>1)
   {
     y=0;
     while ((y<1) && (infloop<1000))
     {
       y=Gam(a,&u,&v,&w)/b;
       infloop += 1;
     }
   }
   else
   {
     if (a>1)
     {
       ga= 0.5*(b-a+sqrt((b-a)*(b-a)+4*b));
       M=pow((a-1)/(b-ga), a-1)*exp(1-a);
       x=0;
       while ((x<1) && (infloop<1))
       {
	 y=-(log(kiss(&u,&v,&w))/ga)+1;
	 infloop += 1;
	 if (kiss(&u,&v,&w)*M < pow(y,a-1)*exp(-(b-ga)*y))
			x=2;
       }
      }
       else
       {
       x=0;
       while ((x<1) && (infloop<1000))
       {
	 y=-(log(kiss(&u,&v,&w))/b)+1;
	 infloop += 1;
	 if (kiss(&u,&v,&w) < pow(y,a-1));
	  x=2;
       }
     }
   }

   if (infloop>999)
     y=1/kiss(&u,&v,&w);

   *i0=u; *j0=v; *k0=w;
   return(y);


}

/*__________________Truncated normal N(0,1), x>a____________________________*/

double Trunc_norm_left(a,i0,j0,k0)
    double a;
    unsigned long int *i0, *j0, *k0;

{
  double x,a_star;
  boolean stop;

  if (a<0.0)
  {
    stop=0;
    while (stop==0)
    {
      nor(&x,&x,i0,j0,k0);
      stop=(x>a);
    }
  }
  else
  {
    a_star=0.5*(a+sqrt(a*a+4.0));
    stop=0;
    while (stop==0)
    {
      x=a-log(kiss(i0,j0,k0))/a_star;
      stop=(log(kiss(i0,j0,k0))<x*(a_star-0.5*x)-a_star*a_star*0.5);
    }
  }

  return(x);
}

/*__________________Optimal truncated normal N(0,1), a<x<b
		    (Robert, Stat&Computing, 1995)        ____________________________*/

double Trunc_norm(u,v,i0,j0,k0)
    double u,v;
    unsigned long int *i0, *j0, *k0;

{
  double x,a,b,boun,boo;
  boolean stop;

  if (v<0)
  {
    a=-v; b=-u;
  }
  else
  {
    a=u; b=v;
  }

  if (b>a+3*exp(-a*a-1/(a*a))/(a+sqrt(a*a+4.0)))
  {
    stop=0;
    while (stop==0)
    {
      x=Trunc_norm_left(a,i0,j0,k0);
      stop=(x<b);
    }
  }
  else
  {
    boo=0.0;
    if (b<0.0)
      boo=b*b;
    if (a>0.0)
      boo=a*a;

    stop=0;
    while (stop==0)
    {
      x=(b-a)*kiss(i0,j0,k0)+a;
      boun=boo-x*x;
      stop=(2.0*log(kiss(i0,j0,k0))<boun);
    }
  }

  if (v>0)
    return(x);
  else
    return(-x);
}

/*__________________Approximative Normal_cdf
		    (Abramowitz&Stegun, 1964)___________________________*/

double normal_cdf(x)
	double  x;

{
   double z,t;

   if (x<0.0)
      t=1.0/(1.0+nrmcd_p*(-x));
   else
      t=1.0/(1.0+nrmcd_p*x);

   z=nrmcd_b5*t;      z=(nrmcd_b4+z)*t;
   z=(nrmcd_b3+z)*t;  z=(nrmcd_b2+z)*t;
   z=(nrmcd_b1+z)*t;  z*=exp(-0.5*x*x)/sqrt_pi;

   if (x>0.0)
	return(1.0-z);
   else
	return(z);
}

/*__________________Approximative Normal_quantile function
		    (Abramowitz&Stegun, 1964)             __________________*/

double normal_quantile(x)
	double x;

{
   double y,z;

   if (x<0.5)
     y=sqrt(-2.0*log(x));
   else
     y=sqrt(-2.0*log(1.0-x));

   z=nrmqtl_a0+nrmqtl_a1*y;
   z/=1.0+(nrmqtl_b1+nrmqtl_b2*y)*y;

   if (x<0.5)
     return(z-y);
   else
          return(y-z);

}


/* _______________________Right truncated gamma on [0,t]________________*/ 

/* The following procedures follow from Anne Philippe's thesis
   and paper (Statistics and Computing, 1997) on truncated gammas.
   See also Robert and Casella (1999, pages 68-69, Exercises 2.44
   and 2.45 for details.		*/

int n_optimal(b)           /* Optimal n */
	double b;
{
  double nr,q;

  q=1.64;
  nr=floor(pow(q+sqrt(q*q+4.*b),2.)/4.-1.);
  return((int)nr);
}


double inter_ri(a,b,i0,j0,k0)
	double a,b;
	unsigned long int *i0,*j0,*k0;
{
  int n,i,j;
  double    *wl,*wlc;
  double     x,u,test=0,y,z,yy,zz;

  n=n_optimal(b);
  wl=(double *)malloc(n+2*sizeof(double));
  wlc=(double *)malloc(n+2*sizeof(double));
  wl[0]=1.;
  wlc[0]=1.;

  for (i=1; i<=n;i++)
  {
	 wl[i]=wl[i-1]*b/(a+i);
	 wlc[i]=wlc[i-1]+wl[i];
		 }

  for(i=0; i<=n; i++) 
	wlc[i]=wlc[i]/wlc[n];

  y=1.0; yy=1.0;

  for (i=1; i<=n; i++) 
  {
    yy=yy*b/i; y=y+yy;
		}

  while (test == 0.0)
  {
      u=kiss(i0,j0,k0);
      j=0;

      while(u>wlc[j])
	j++;

      x=Beta(a,(double)(j),i0,j0,k0);
      zz=1.0; z=1.0;

      for (i=1; i<=n; i++) {
           zz*=(1.-x)*b/i;
	   z+=zz;
	       }

       y*=exp(-b*x)/z;

       if (kiss(i0,j0,k0)<y)
	   test=1.0;
                }
		 
   free(wl);
   free(wlc);
   return(x);
}

double gamma_right(a,b,t,i0,j0,k0)
	double a,b,t;
        unsigned long int *i0,*j0,*k0;
{ 
	return(inter_ri(a,b*t,i0,j0,k0)*t);
};



/* _______________________left truncated gamma on [t,infinity[
			  (Anne Philippe, 1996)_____________________________*/


double integer(a,b,i0,j0,k0)
	double a,b;
        unsigned long int *i0,*j0,*k0;
{
  double *wl,*wlc,u,x;
  int i;

  wl=(double *)malloc((a+1)*sizeof(double));
  wlc=(double *)malloc((a+1)*sizeof(double));

  wl[1]=1.0;
  wlc[1]=1.0;

  for(i=2; i<=(int)(a) ;i++ )
  {
	wl[i]=wl[i-1]*(a-i+1)/b;
	wlc[i]=wlc[i-1]+wl[i];
		  }

  for(i=1; i<=(int)(a); i++)
       wlc[i]/=wlc[(int)(a)];

  u=kiss(i0,j0,k0);

  i=1;

  while(u>wlc[i])
	i++;

  x=Gam((double)(i),i0,j0,k0)/b+1.0;

  free(wl);
  free(wlc);
  return(x);
}

double inter_le(a,b,i0,j0,k0)
	double a,b;
	unsigned long int *i0,*j0,*k0;
{
  double test=0.0,u,x,y,M;

  if (a<1.0) { 
	M=1.0;

	while (test == 0.0){ 
	   x=1.-(1./b)*log(1.-kiss(i0,j0,k0));
	   y=1./pow(x,1.-a);

	   if (kiss(i0,j0,k0)< y/M) test=1.0;
			      }
  }
  else
  { 
 	if (a<b){ 
	  M=exp(floor(a)-a);

	  while (test == 0){ 
	     x=integer(floor(a), b*floor(a)/a,i0,j0,k0);
	     y=pow(x, a-floor(a))*exp(-x*b*(1.-floor(a)/a));

	     if (kiss()< y/M) test=1.0;
			}
		}
	else
	{
 	  M=exp(floor(a)-a)*pow(a/b,a-floor(a));

	  while (test == 0.0){ 
	     x=integer(floor(a), b+floor(a)-a,i0,j0,k0);
	     y=pow(x, a-floor(a))*exp(-x*(-floor(a)+a));

	     if (kiss()< y/M) 
			test=1.0;
			}
		}
	}

  return(x);
}

double gamma_left(a,b,t,i0,j0,k0)
	double a,b,t;
        unsigned long int *i0,*j0,*k0;
{
	return(inter_le(a, b*t,i0,j0,k0)*t);
}