/********************************************************************/
/**                                                                **/
/**   tdrg.c                                                       **/
/**                                                                **/
/**   A universal generator for multivariate T-concave             **/
/**   distributions                                                **/
/**                                                                **/
/**   generators for (truncated) gamma distributed variates        **/
/**                                                                **/
/**   authors: Josef Leydold, Guenter Tirler and Gerhard Derflinger**/
/**   last modified: Mon Mar  9 13:38:00 MET 1998                  **/
/**                                                                **/
/********************************************************************/              

/* the header file */
#include "tdrmv.h"

/*-----------------------------------------------------------------*/
#if GAMMA == 2
/*-----------------------------------------------------------------*/

/*******************************************************************

  tdrg()
  
  inversion-rejection method see

  @Article{Hoermann:95toms,
  author =       {W. H{\o}rmann},
  title =        {A rejection technique for sampling from {$T$}-concave distributions},
  journal =      {ACM Trans. Math. Software},
  year =         {1995},
  volume =       {21},
  number =       {2},
  pages =        {182--193}
  }

  ******************************************************************/


/** local definitions **/

#define GNTP (TDRG_NTP+TDRG_NTP_L)  /* maximum number of touching points for gamma distribution */
#define ST_SIZE 5                   /* relative size of search table for hat of gamma distribution */

/* structures */
static struct {                     /* data for hat of gamma distribution */
  double tp[GNTP];           /* coordinates of touching points */
  double logf[GNTP];         /* log(f(x)) at touching point of gamma density */
  double derlogf[GNTP];      /* derivative of log(f(x)) at touching point of gamma density */
  double bound[GNTP];        /* boundaries between touching points */
  double squeeze[GNTP-1];    /* slopes of squeeze */
  double vol[GNTP];          /* accumulated volume below hat */
  int st_size;               /* size of search table */
  int *st;                   /* search table for volumes */
} ghat;

/*-----------------------------------------------------------------*/

/* functions */
void tdrg_setup(C_ROOT *cr);           /* setup for generator of gamma distributed variates with shape parameter N */
void db_tdrg(C_ROOT *cr);              /* print parameters into log files */
double tdrg(C_ROOT *cr, CONE *c);   
/* generate variates with gamma distribution with shape parameters N and c->beta */

/* find optimal touching points */
static void optPoint(double (*density)(double x),double *bd,double cc,int nb,double x0,int mode,double step,
		     int in,double *xi);
static int  setas(double (*density)(double x),double *s,double *etas,double *alt,double *x,double *y,double *bd,
		  char *undef,double *dmax,int k,int j,int *mj,double cc, int in,char *lbd,double *brd);
static double opt_integrate(double *x, double *y, double *etas, int j, double cc);
static double area(double (*density)(double x),double *s,double *etas,double *x,double *y, double *bd,char *lbd,
		int *iou,int k,int m,int nb,double cc);
static double brf(double a,double b, double c, double p, double q, double r);
static double w(double x1,double x2,double xm,double y1,double y2,double ym);
static double dgamma(double x);    /* density function of gamma distribution ( * constant ) */
static void tdr_find_mode(double (*density)(double), char *lbd, double *bd, double *x, double *y, int mode, double step);
/* find mode of univariate density */
static void tdr_find_bracket(double (*density)(double), char *lbd, double *bd, double *x, double *y, double cc);
/* search for proper bracket or mode */

/*-----------------------------------------------------------------*/

extern void db_print_string(char *); /* print character string into log file */
extern double uniform();             /* uniform random number generator */
extern void warning(char *msg, char *filename, int line); /* write erro message */
extern void fatal(char *msg, char *filename, int line); /* fatal error, exit */

/*-----------------------------------------------------------------*/

/*******************************************************************/

/*-----------------------------------------------------------------*/

void tdrg_setup(C_ROOT *cr)
/*
  setup for generator of gamma distributed variates
*/
{
  int i,j;           /* aux counter */
  double bd[3];      /* bounds for domain of gamma distribution (see optPoint) */
  double mode;       /* mode of distribution */
  int ntp;           /* total number of construction points */

#if RECTANGLE == 1
  double max,tmp;    /* maximal height of pyramides */
  int *stx;          /* search table for finding interval which contains a given x */
  CONE *c;
  
  /* find maximal height of pyramides */
  max = 0.;
  for( c=cr->c; c!=NULL; c=c->next ) {   /* all cones */
    tmp = c->height * c->beta;   /* we have to consider this parameter (see rp_get_point()) */
    max = MAX(max,tmp);
    if( c == cr->last_c ) break;
  }

  /* initialize some parameters */
  bd[0] = 0.;       /* left boundary = 0. */
  bd[1] = -1.;       /* ignore */
  bd[2] = max;    /* right boundary = infinity */
  mode = ( (N-1.)<=max ) ? N-1. : max ;  /* the mode */
  /* find optimal touching points */
  if( TDRG_NTP_L < 0 ) FATAL("TDRG_NTP_L must be an integer >= 0");
  optPoint(dgamma,bd,0.,TDRG_NTP,mode,1,mode/2.,6,ghat.tp+TDRG_NTP_L);
  /* find additional points for the left part of the density */
  if( TDRG_NTP_L > 0 )
    ghat.tp[0] = ghat.tp[TDRG_NTP_L] * 0.5;
  if( TDRG_NTP_L > 1 ) {
    bd[2] = ghat.tp[TDRG_NTP_L] * TDRG_NTP_L / (TDRG_NTP_L+1.);
    mode = bd[2];  /* the mode */
    optPoint(dgamma,bd,0.,TDRG_NTP_L,mode,1,mode/2.,6,ghat.tp);
  }
  /* total number of construction points */
  ntp = TDRG_NTP+TDRG_NTP_L;
  /* store maximum height of pyramids */
  cr->max_gamma = max;    
#else
  /* initialize some parameters */
  bd[0] = 0.;       /* left boundary = 0. */
  bd[1] = 1.;       /* ignore */
  bd[2] = bd[1];    /* right boundary = infinity */
  mode = N-1.;  /* the mode */
  /* find optimal touching points */
  optPoint(dgamma,bd,0.,TDRG_NTP,mode,1,mode/2.,6,ghat.tp);
  /* total number of construction points */
  ntp = TDRG_NTP;
#endif

  /* store number of touching points */
  cr->tdrg_ntp = ntp;

  /* calculate (log(f(x)))' at touching point */
  for( i=0; i<ntp; i++ )
    while( 1 ) {
      ghat.derlogf[i] = (N-1.)/ghat.tp[i] - 1.;
      if( fabs(ghat.derlogf[i]) > 1.e-8 )
	/* o.k. */
	break;
      else
	/* we do not want to deal with the special case (log(f(x)))' == 0 */
	ghat.tp[i] -= 1.e-4*ghat.tp[i];
    }

  /* calculate log(f(x)) at touching point */
  for( i=0; i<ntp; i++ )
    ghat.logf[i] = (N-1.)*log(ghat.tp[i]) - ghat.tp[i];
  
  /* calculate boundaries between touching points */
  ghat.bound[0] = 0.;
  for( i=1; i<ntp; i++ )
    ghat.bound[i] = ( ghat.logf[i] - ghat.logf[i-1] - ghat.tp[i]*ghat.derlogf[i] + ghat.tp[i-1]*ghat.derlogf[i-1] ) /
      ( ghat.derlogf[i-1] - ghat.derlogf[i] );

  /* calculate slope of squeeze between touching points */
  for( i=0; i<ntp-1; i++ )
    ghat.squeeze[i] = ( ghat.logf[i+1] - ghat.logf[i] ) / ( ghat.tp[i+1] - ghat.tp[i] );

  /* calculate volume below hat */
  for( i=0; i<ntp-1; i++ )
    /* bounded regions */
    ghat.vol[i] = ( exp( ghat.logf[i] + ghat.derlogf[i] * ( ghat.bound[i+1] - ghat.tp[i] ) )
		    - exp( ghat.logf[i] + ghat.derlogf[i] * ( ghat.bound[i] - ghat.tp[i] ) ) ) / ghat.derlogf[i];

  /* last region */
#if RECTANGLE == 1
  /* bounded region */
  i = ntp-1;
  ghat.vol[i] = ( exp( ghat.logf[i] + ghat.derlogf[i] * ( max - ghat.tp[i] ) )
		  - exp( ghat.logf[i] + ghat.derlogf[i] * ( ghat.bound[i] - ghat.tp[i] ) ) ) / ghat.derlogf[i];
#else
  /* unbounded region */
  ghat.vol[ntp-1] = - exp( ghat.logf[ntp-1] + ghat.derlogf[ntp-1] * ( ghat.bound[ntp-1] - ghat.tp[ntp-1] ) )
    / ghat.derlogf[ntp-1];
#endif

  /* accumulated sum of volumes */
  for( i=1; i<ntp; i++ )
    ghat.vol[i] += ghat.vol[i-1];

  /* make a search table for volumes */
  /* size of search table */
  ghat.st_size = ST_SIZE*ntp;
  /* allocate memory */
  ghat.st = (int *)malloc((size_t) (ghat.st_size*sizeof(int)));
  if (!ghat.st) FATAL("allocation failure");
  /* create the table */
  for( j=0,i=0; j<ghat.st_size; j++ ) {
    while( i<ntp-1 && ( ghat.vol[i]*ghat.st_size < (j+1)*ghat.vol[ntp-1] ) )
      i++;
    ghat.st[j] = i;
  }

#if RECTANGLE == 1
  /* now calculate the volume of the hat above the domain of the truncated gamma distribution */

  /* make a search table for x values */
  /* allocate memory */
  stx = (int *)malloc((size_t) (ghat.st_size*sizeof(int)));
  if (!stx) FATAL("allocation failure");
  /* create the table */
  for( j=0,i=0; j<ghat.st_size; j++ ) {
    while( i<ntp-1 && ( ghat.bound[i+1]*ghat.st_size < (j+1)*max ) )
      i++;
    stx[j] = i;
  }


  for( c=cr->c; c!=NULL; c=c->next ) {
    /* upper bound of domain of truncated gamma density */
    max = c->height * c->beta;    /* we reuse the variable max */

    /* find greatest intervall that covers domain */
    i = ((int)(ghat.st_size*max/cr->max_gamma));    /* aux value (thus we do not need another variable */
    if( i >= ghat.st_size ) i = ghat.st_size-1;
    j = stx[i];
    for( ; j>=0 && ghat.bound[j]>max; j-- ) ;

    /* calculate volume above last interval */
    c->tdrg_vol = ( exp( ghat.logf[j] + ghat.derlogf[j] * ( max - ghat.tp[j] ) )
		    - exp( ghat.logf[j] + ghat.derlogf[j] * ( ghat.bound[j] - ghat.tp[j] ) ) ) / ghat.derlogf[j];
    /* compute volume above domain */
    c->tdrg_vol += (j) ? ghat.vol[j-1] : 0.;

    if( c == cr->last_c ) break;
  }

  /* we do not need the search table stx any more */
  free( stx );

#endif  

} /* end of tdrg_setup() */

/*-----------------------------------------------------------------*/

#if DEBUG > 0
void db_tdrg(C_ROOT *cr)
/*
  print parameters into log files.
  we have to compose the log message, since ghat is local.
*/
{
    char s[100];
    int i;

    db_print_string("\nTDRG Gamma random number generator:\t");
#if RECTANGLE == 1
    sprintf(s,"(%d+%d touching points for hat)\n\n",TDRG_NTP,TDRG_NTP_L);
#else
    sprintf(s,"(%d touching points for hat)\n\n",TDRG_NTP);
#endif
    db_print_string(s);
    db_print_string("  left bound         x_i       log(f(x_i))   (log(f(x_i)))'     volume        squeeze\n\n");

    for( i=0; i<cr->tdrg_ntp-1; i++ ) {
      sprintf(s," %10f     %10f     %10f     %10f     %10f     %10f\n",
	      ghat.bound[i], ghat.tp[i], ghat.logf[i], ghat.derlogf[i], ghat.vol[i], ghat.squeeze[i]);
      db_print_string(s);
    }
    sprintf(s," %10f     %10f     %10f     %10f     %10f\n",
	    ghat.bound[i], ghat.tp[i], ghat.logf[i], ghat.derlogf[i], ghat.vol[i]);
    db_print_string(s);

    db_print_string("\n==================================================\n\n");
   
} /* end of db_tdrg() */
#endif

/*-----------------------------------------------------------------*/

double tdrg(C_ROOT *cr, CONE *c)
/*
  generate variates with gamma distribution with shape parameters N and c->beta
  density: f(x) = x^(N-1) exp(- beta * x)    ( * const )
*/
{
  double u;       /* uniformly distributed random point */
  double x;       /* gamma distributed random point */
  int in;         /* number of intervall in domain of hat */
  double hy;      /* value of transformed hat at x: T(h(x)) */
  double hk, hd;  /* slope and constant in term for linear hat */
  double vol;     /* volume below hat function */

  /* set volume below hat function */
#if RECTANGLE == 1
  vol = c->tdrg_vol;
#else
  vol = ghat.vol[cr->tdrg_ntp-1];
#endif

  while( 1 ) {    /* repeat untill point is accepted */

#if RECTANGLE == 1
    /* uniform random number: 0 <= u < volume */
    u = uniform() * vol;
    /* find interval in search table */
    in = ghat.st[(int)(ghat.st_size*u/ghat.vol[cr->tdrg_ntp-1])];
    /* find intervall (linear search from top) */
    for( ; in>0 && ghat.vol[in-1]>u; in-- ) ;
#else
    /* uniform random number: 0 <= u < 1 */
    u = uniform();
    /* find interval in search table */
    in = ghat.st[(int)(ghat.st_size*u)];
    /* uniform random number: 0 <= u <= Volume(hat) */
    u *= vol;
    /* find intervall (linear search from top) */
    for( ; in>0 && ghat.vol[in-1]>u; in-- ) ;
#endif

    /* reuse random number */
    /* 0 <= u < Volume(hat over interval) */
    if( in > 0 )
      u = u-ghat.vol[in-1];
    
    /* get random point: x = H^{-1}(u) */
    hk = ghat.derlogf[in];
    hd = ghat.logf[in]-ghat.tp[in]*ghat.derlogf[in];
    x = ( log( hk * u  + exp( hk * ghat.bound[in] + hd ) ) - hd ) / hk;

    /** now accept or reject **/

    /* linear hat (of transformed density) */
    hy = hk * x + hd;
    
    /* random point in -00 < u <= log hy */
    u = hy + log( uniform() );
    
    /* point (x,exp(u)) below (transformed) squeeze ? */
    if( x < ghat.tp[in] ) {   /* left of touching point */
      if( in>0 && ghat.logf[in-1]+ghat.squeeze[in-1]*(x-ghat.tp[in-1])>=u ){
#if DEBUG > 0
	if( DEBUG & DB_GAMMA )
	  ++cr->db_n_tdrg_accept;
#endif
	return x/(c->beta);
      }
    }
    else {                    /* right of touching point */
      if( in<cr->tdrg_ntp && ghat.logf[in]+ghat.squeeze[in]*(x-ghat.tp[in])>=u ) {
#if DEBUG > 0
	if( DEBUG & DB_GAMMA )
	  ++cr->db_n_tdrg_accept;
#endif
	return x/(c->beta);
      }
    }
    
    /* is point below (transformed) density */
    if( u<= ((double)(N-1))*log(x)-x ) {
#if DEBUG > 0
      if( DEBUG & DB_GAMMA )
	++cr->db_n_tdrg_accept;
#endif
      return x/(c->beta);
    }

#if DEBUG > 0
    if( DEBUG & DB_GAMMA )
      ++cr->db_n_tdrg_reject;
#endif

  }

} /* end of tdrg() */

/*-----------------------------------------------------------------*/

/*******************************************************************/

#define   NNP          100                      /*  maximum number of points left and right of the mode */
#define   max_NNP      2*NNP+1                  /*  maximum number of points   */

/*-------------------------------------------------------------------------------------------------------*/

static void optPoint(double (*density)(double x),double *bd,double cc,int nb,double x0,int mode, double step,int in,double *xi)
     /* This functin is called by "setup"
        This program calculates the asymtotic optimal points for a hat function and squeezes 
	It implements
		"THE OPTIMAL SELECTION OF HAT FUNCTIONS FOR REJECTION ALGORITHMS" ,
		     Gerhard DERFLINGER and Wolfgang HOERMANN,
		 Institute f Statistics,University of Busness Administration Vienna,
		 Augasse 2-6, A-1090 Vienna, AUSTRIA.     
	Version 1.1
	date of last correction: 3.3.98
	written in C (Compiled with GNU C Compiler "gcc", the digital C Compiler  and the C Compiler of Borland)
	written by Guenter Tirler,
		     Institute f Statistics,
		     University of Busness Administration Vienna,
		     Augasse 2-6, A-1090 Vienna, AUSTRIA

        modified by Josef Leydold		     

	This function calls the functions
		    "void   setas(...)"
		    "double area(...)"
		    "double w(...)"
                    "double brf(...)"	

	 input data:
	  density            : density function 
          bd[0],bd[1],bd[2]  : 3 values which describe the support of the distribution:
                                  bd[0] is the left boundary of the support
                                  bd[2] is the right boundary of the support
                                  if bd[0]=bd[1] : left boundary is infinite
                                  if bd[1]=bd[2] : right boundary is infinite					 
	   double   cc        :  parameter of the transformation
	   int      nb        :  number of asymtotic optimal points
	   double   x0        :  starting point : (density(x0)>0)
	   int      mode      :  not "0" (e.g.:1)  if  x0 is the  mode.
				     "0"           if  x0 is not the mode
	   double   step      :  step to find mode  
	   int      in        :  number of intervals (important for the numerical integration)

	output data:
	     double *xi         :  array of opt.points: xi[0],...,xi[nb-1]


	this function uses following libraries: <stdio>   : printf()
                                                <stdlib.h>: exit()
			        		<math.h>  : exp(), log()                              */
{

#define   dritl      1.0/3.0                 

		       
  double x[max_NNP];              /* optimal values of hat functions and squeezes, x[points_lr] = mode */
  double y[max_NNP];              /* distribution at the value x   */
  double s[max_NNP];              /* compare to function 'F' in algorithm (2)       */
  double etas[max_NNP];           /* first derivation of the transformed density   */
  double wm[3];                   /*        */
  double brd[3];
  char   lbd[3];                  /* lbd[0] status of the left boundary of the support : "i"(infinite) or "f"(finite) */
		  		  /* lbd[2] status of the right boundary of the support: "i"(infinite) or "f"(finite) */
  char   undef[3]={'f','f','f'};  /*                       "t" for "true" and "f" for "false"  */
  char   trunc[3];
  double dmax[3];                 /* number of intervals for numerical integration*/
  double alt[3];
  int    mj[3]={-NNP,0,NNP};      /* number of integration steps left and right from mode*/
  int    iou[3];
  double hof;                     /* constant in equation (7) */
  int    k;                       /* step left ("-1") or right("1")  */
  int    j;
  int    jj;
  int    ih;
  int    io;
  double dlo;
  double wen;
  double urdlo;
  double clo;
  double elo;
  double deltaz;
  double z;

  const double bord=1000.0;       /* number of 'hoermann'-distances where the program stops */


  /** initialization and several definitions **/

  for( k=0; k<max_NNP; k++)
    y[k] = -1.0;
  for( k=0; k<NNP; k++) {
    s[k] = 1.0;
    s[NNP+k+1] = -1.0;
  }
  s[NNP]=0.0;

  if (bd[0] != bd[1]) lbd[0] ='f'; else lbd[0]='i';    /* status of the left boundary of the support (finite/infinite) */
  if (bd[1] != bd[2]) lbd[2] ='f'; else lbd[2]='i';    /* status of the right boundary of the support (finite/infinite)*/

  hof = ( cc != 0. ) ? 0.25 : exp(-1.);    /*  compare to equation (7). 0 < hof < 0,368... */

  /* check location of starting point */
  if ( ((lbd[0]=='f') && (x0<bd[0])) || ((lbd[2]=='f') && (x0>bd[2])) )
    FATAL("error 1:point outside");


  /** find mode of density and bracket **/
  x[NNP] = x0;      /* starting for searching / mode if known */
  tdr_find_mode(density,lbd,bd,x+NNP-1,y+NNP-1,mode,step);
  /* x[NNP] is the mode of density, x[NNP-1] <= x[NNP] <= x[NNP+1] with y[NNP] >= y[NNP-1] and y[NNP] >= y[NNP+1] */

  /* Next step: find bracket such that
    hu < y[NNP+k]=density(x[NNP+k]) < ho  for k=1,-1.
    hu=0.64*hof, ho =hof/0.64
    "0.64" is a heuristical choice */
  tdr_find_bracket(density,lbd,bd,x+NNP-1,y+NNP-1,cc);

  /* now make correction corrections for special case */
  for( k=-1; k<=1; k=k+2) {
    if( (x[NNP+k]==x[NNP]) || (y[NNP+k]>=y[NNP]) ) {
      /* mode = boundary --> point is changed */
      x[NNP] = (2.0*x[NNP+k]+x[NNP-k]) / 3.;
      y[NNP] = density(x[NNP]);
    }
  }

  /** integration **/
  /* set parameters for numerical integration */
  for( k=-1; k<=1; k=k+2) {

    /* step left and right                       */
    /* heuristical  correction factor sqrt(...)  */
    brd[k+1] = (x[NNP+k]-x[NNP])*sqrt(y[NNP+k]/y[NNP]/hof); 

    /* maximum length of intervals for numerical integration*/
    dmax[k+1] = brd[k+1]/((double) in);                      

    brd[k+1] *= bord;

    if( lbd[k+1]=='f' ) {                /* bounded in direction k */
      /* number of integration steps left and right from mode*/
      mj[k+1] = ((int)((bd[k+1]-x[NNP])/dmax[k+1]));  
      mj[k+1] = k*(mj[k+1]+2);
      /* integration interval is devided to mj[k+1] parts of the left/right support*/
      dmax[k+1] = k*(bd[k+1]-x[NNP])/mj[k+1];
    }
  }

  /* update bracket to length of integration intervals */
  for( k=-1; k<=1; k=k+2 ) {
    x[NNP+k]=x[NNP]+dmax[k+1];
    y[NNP+k]=density(x[NNP+k]);
  }

  /* calculate parameters for hat function */
  alt[0] = opt_integrate(x,y,etas,0,cc);
  alt[2] = alt[0];

  for( k=-1; k<=1; k=k+2) {
    for( j=1; j<=2; j++ )
      setas(density,s,etas,alt,x,y,bd,undef,dmax,k,k*j,mj,cc,in,lbd,brd);   /* j=-1,-2 (k=-1) und j=1,2 (k=1) */
    iou[k+1] = 2*k;
    trunc[k+1] = 'f';
  }

  dlo = area(density,s,etas,x,y,bd,lbd,iou,1,0,nb,cc);

  do {

    urdlo = dlo;
    for( k=-1; k<=1; k=k+2) {
      if( trunc[k+1] == 't' ) continue;

      iou[k+1] = iou[k+1] + k;
      if( k*iou[k+1] == NNP )
	FATAL("error 4: dimension NNP exceeded");

      if( brd[k+1]*(x[NNP+iou[k+1]]-brd[k+1]) > 0.0 ) {
	trunc[k+1] = 't';
	iou[k+1] = iou[k+1]-k;
	WARNING("search for optimal left/right point stopped");
	continue;         
      }
      if( iou[k+1]==(mj[k+1]+k) ) 
	wen = dlo;
      else {
	if( (k*s[iou[k+1]+NNP]) < 0.0 ) {
	  if( !setas(density,s,etas,alt,x,y,bd,undef,dmax,k,iou[k+1],mj,cc,in,lbd,brd) ) {
	    /* Error */
	    trunc[k+1] = 't';
	    iou[k+1] = iou[k+1]-k;
	    WARNING("search for optimal left/right point stopped");
	    continue;   
	  }
	}
	wen = area(density,s,etas,x,y,bd,lbd,iou,k,0,nb,cc);
      }
      if( wen < dlo )
	dlo = wen;
      else {
	iou[k+1] = iou[k+1]-2*k;
	wen = area(density,s,etas,x,y,bd,lbd,iou,k,0,nb,cc);
	if( wen < dlo )
	  dlo = wen;
	else
	  iou[k+1] = iou[k+1]+k;
      }
    }
  } while (dlo<urdlo);

  for (k=-1;k<=1;k=k+2) {
    j=iou[k+1];
    ih=(nb+1+k*(nb-1))/2;
    if (trunc[k+1]=='t') {
      j=j-1; 
      xi[ih-1]=x[j+NNP];       
      wm[k+1]=s[j+NNP];
    }
    else {
      if (iou[k+1]==mj[k+1]) {
	xi[ih-1]=bd[k+1];
	wm[k+1]=s[NNP+iou[k+1]];
      }
      else {
	clo=area(density,s,etas,x,y,bd,lbd,iou,k,-1,nb,cc);
	elo=area(density,s,etas,x,y,bd,lbd,iou,k,1,nb,cc);
	xi[ih-1]=brf(x[j-1+NNP],x[j+NNP],x[j+1+NNP],clo,dlo,elo);
	if ( (xi[ih-1]-x[j+NNP])<0.0 ) 
	  jj=-1;
	else
	  jj=1;
	wm[k+1]=s[j+NNP]+(xi[ih-1]-x[j+NNP])*(s[j+jj+NNP]-s[j+NNP])/(x[j+jj+NNP]-x[j+NNP]);
      }
    }
  }

  deltaz=(wm[2]-wm[0])/(nb-1);
  z=wm[0]+deltaz;
  io=iou[0];
  for (j=2;j<=nb-1;j=j+1) {
    do
      io++;
    while (s[io+NNP]<=z);
    xi[j-1]=x[io-1+NNP]+(z-s[io-1+NNP])/(s[io+NNP]-s[io-1+NNP])*(x[io+NNP]-x[io-1+NNP]);
    io=io-1;
    z=z+deltaz;
  }

} /* end of optPoint() */

/*--------------------------------------------------------------------------------------------------------*/

static int setas(double (*density)(double x),double *s,double *etas,double *alt,double *x,double *y,double *bd,
             char *undef,double *dmax,int k,int j,int *mj,double cc, int in,char *lbd,double *brd)
/* this function is called by "optPoint"
   s[j+NNP] (compare to function 'F' in algorithm (2) ) 
   and etas[j+NNP] (first derivation of the transformed density) are computed
   x[j-1], y[j-1] and  x[j+1], y[j+1] for second derivation are computed.                
*/
{
  double aneu;
  double eps;
  const int   ll=120;
  const double fd=0.01;
  const double ff=1.0+fd;


  if( j == mj[k+1] ) {
    /* outside of region for numerical region --> extrapolation */
    s[j+NNP] = 3.0*s[j-k+NNP] - 3.0*s[j-2*k+NNP] + s[j-3*k+NNP];    /* extrapolation */
    etas[j+NNP] = 2.0*etas[j-k+NNP] - etas[j-2*k+NNP];              /* extrapolation for derivation */
    return 1;  /* o.k. */
  }

  /* region of numerical integration */

  if( undef[k+1]=='t' )
    x[j+k+NNP] = (x[j+NNP]+bd[k+1])/2.0;

  else {
    if( k*(mj[k+1]-j) > 1 ) {
      x[j+k+NNP] = x[j+NNP]+dmax[k+1];

      if ( brd[k+1]*(x[j+k+NNP]-brd[k+1]) > 0.0 )
	return 0;    /* Error */

      if ( (k*j) >= (ll*in) ) {
	dmax[k+1] = ff*dmax[k+1];
	if ( ((k*j)==(ll*in)) && (lbd[k+1]=='f') ) {
	  mj[k+1] = k+j+k*log(1.0+fd*(bd[k+1]-x[j+k+NNP])/(ff*dmax[k+1]))/log(fd);
	  if ( (k*mj[k+1]) < NNP )
	    mj[k+1] = k*k*mj[k+1];
	  else
	    mj[k+1] = k*NNP;
	}
      }
    }
    else if ( (j+k)==mj[k+1] ) {
      if (lbd[k+1]=='i')
	FATAL("error 4:dimension nd exceeded");
      else {
	x[j+k+NNP] = bd[k+1];
	y[j+k+NNP] = density(x[j+k+NNP]);
	eps = 0.000001*dmax[k+1];
	if ( (y[j+k+NNP]==0.0) || ( (k*(density(x[j+k+NNP]-eps)-y[j+k+NNP])/eps) > (1000.0*y[NNP]) ) )
	  undef[k+1] = 't';
	else
	  undef[k+1] = 'f';
	if (undef[k+1]=='t') {
	  x[j+k+NNP] = (x[j+NNP]+bd[k+1])/2.0;
	  mj[k+1] = k*NNP;
	}
	else {
	  goto m1;
	}
      }
    }
  }
  
  y[j+k+NNP]=density(x[j+k+NNP]);                                           /* new x[NNP-3]...x[NNP+3] */
m1:
  aneu = opt_integrate(x,y,etas,j,cc);
  s[j+NNP] = s[j-k+NNP]+(alt[k+1]+aneu)*(x[j+NNP]-x[j-k+NNP])/2.0;                /* trapezoid formula */
  alt[k+1] = aneu;

  return 1;  /* o.k. */
} /* end of setas() */

/*--------------------------------------------------------------------------------------------------------*/

static double opt_integrate(double *x, double *y, double *etas, int j, double cc)
/*
  integrate
*/
{
  double hl,hr;       /* aux variable necessary to compute ys and yss */
  double ys;          /* numerical approximation for first derivative */
  double yss;         /* numerical approximation for second derivative */
  double temp;

  /* approximation of first and second derivative */
  hr = (y[j+1+NNP]-y[j+NNP]) / (x[j+1+NNP]-x[j+NNP]) / (x[j+1+NNP]-x[j-1+NNP]);
  hl = (y[j+NNP]-y[j-1+NNP]) / (x[j+NNP]-x[j-1+NNP]) / (x[j+1+NNP]-x[j-1+NNP]);
  ys = (x[j+NNP]-x[j-1+NNP])*hr + (x[j+1+NNP]-x[j+NNP])*hl;    /* approximation for the first derivative of density */  
  yss= 2.0*(hr-hl);                                            /* approximation for the second derivative of density */

  /* first derivation of the transformed density: T'(f(x)) */
  if (cc!=0.0)
    etas[j+NNP]=(-cc)*exp((cc-1.0)*log(y[j+NNP]))*ys; 
  else         
    etas[j+NNP]=ys/y[j+NNP];

  /* integrating function */
  temp = (1.0-cc)*ys*ys/y[j+NNP] - yss;                            /*  t(x) */
  if (temp<0.0) 
    FATAL(" error 3:  f(x) not T-concave");

  return( (temp==0.) ? 0 : pow(temp,dritl)/2.0 );

} /* end of opt_integrate() */

/*--------------------------------------------------------------------------------------------------------*/

static double area(double (*density)(double x),double *s,double *etas,double *x,double *y, double *bd,char *lbd,int *iou,int k,int m,int nb,double cc)
/*
  this function is called by "optPoint"               
*/
{
  int     kk;
  int     j;
  double  c1;
  double  eta;          /* transformed density */
  double  etabd;
  double  ar;           /* value of the function */

  iou[k+1]=iou[k+1]+m;
  ar=(s[NNP+iou[2]]-s[NNP+iou[0]])*(s[NNP+iou[2]]-s[NNP+iou[0]])*(s[NNP+iou[2]]-s[NNP+iou[0]])/(nb-1)/(nb-1);
  for (kk=-1;kk<=1;kk=kk+2) {
    j=iou[kk+1];
    if (cc!=0.0) {
      c1=(cc+1.0)/cc;
      eta=exp(cc*log(y[j+NNP]));
      if (lbd[kk+1]=='f') {
	etabd=eta-etas[j+NNP]*(bd[kk+1]-x[j+NNP]);
	if (cc==(-0.50)) {
	  ar=ar+kk*(bd[kk+1]-x[j+NNP])/(etabd*eta);
	}
	else {
	  ar=ar+kk*(exp(c1*log(eta))-exp(c1*log(etabd)))/(c1*etas[j+NNP]);
	}
      }
      else 
	ar=ar+kk*exp(c1*log(eta))/(c1*etas[j+NNP]);
    }
    else {
      if (lbd[kk+1]=='f') {
	if (etas[j+NNP]!=0.0) {
	  etabd=log(y[j+NNP])+etas[j+NNP]*(bd[kk+1]-x[j+NNP]);
	  ar=ar-kk*(y[j+NNP]-exp(etabd))/etas[j+NNP];
	}
	else {
	  ar=ar+kk*(bd[kk+1]-x[j+NNP])*y[j+NNP];
	}
      }
      else {
	ar=ar-kk*y[j+NNP]/etas[j+NNP];
      }
    }
  }
  iou[k+1]=iou[k+1]-m;
  return ar;
} /* end of area */

/*--------------------------------------------------------------------------------------------------------*/

static double w(double x1,double x2,double xm,double y1,double y2,double ym)
/*
  interpolation function
 */
{
  return xm+(x1-xm) * exp( 1.0/(1-(log(ym/y2-1.0)/log(ym/y1-1.0)) ) * log((x2-xm)/(x1-xm)) );
} /* end of w() */

/*--------------------------------------------------------------------------------------------------------*/

static double brf(double a,double b, double c, double p, double q, double r)          
/*
  interpolation function
 */
{
  return ((p*(b*b-c*c)+q*(c*c-a*a)+r*(a*a-b*b))/(p*(b-c)+q*(c-a)+r*(a-b))/(2.0));
} /* end of brf() */

/*--------------------------------------------------------------------------------------------------------*/

static double dgamma(double x)
/* 
   density function of gamma distribution ( * constant)
*/
{
  int i;
  double factor = 1;

  for( i=1; i<N; i++ )
    factor *= x;
  return( factor * exp(-x) );

} /* end of dgamma() */

/*--------------------------------------------------------------------------------------------------------*/

static void tdr_find_mode(double (*density)(double), char *lbd, double *bd, double *x, double *y, int mode, double step)
/* 
   search for mode of density
   
   density ... pointer to density function
   lbd[3]  ... status of left (lbd[0]) and right (lbd[2]) bound. 'i' inifinite, 'f' finite.
   bd[3]   ... boundary of interval. bd[0] left, bd[2] right.
   x[3]    ... touching points x[0] <= x[1] <= x[2]. x[1] starting point.  (bracket for mode)
   y[3]    ... density at x[i]
   mode    ... whether x[1] is mode or not.
   step    ... stepsize in searching for mode.
*/
{
  int k;               /* direction of search */
  double r;            /* ratio of values of density */
  double xt, yt;       /* temp values for x and f(x) */

  /* density at starting point */
  y[1] = density( x[1] );
  if( y[1] == 0. )
    FATAL("error 2: f(x)=0.0 at the starting point!");

  /* initialize y */
  y[0] = -1.;
  y[2] = -1.;

  /* check direction of searching */
  if( (lbd[2]=='f') && (x[1]==bd[2]) && step > 0. )
    /* right boundary is finite and is equal to starting point --> goto left */
    step *= -1.;
  if( (lbd[0]=='f') && (x[1]==bd[0]) && step < 0. )
    /* left boundary is finite and is equal to starting point --> goto right */
    step *= -1.;
  /* check input */
  if( step == 0. ) {
    WARNING(" stepsize must be != 0.");
    step = 1.;
  }

  /* find a bracket for maximum: x[0] <= x[1] <= x[2] with f(x[1]) > f(x[0]) and f(x[1]) > f(x[2]) */
  while( 1 ) {

    /* direction of search */
    /* k = -1 --> left; k = 1 --> right */
    k = ( step < 0.0 ) ? -1 : 1;

    while( 1 ) {
      /* find a point left/right of x[1] */

      /* make a proper step */
      while( 1 ) {
	/* jump */
	x[1+k] = x[1] + step;
	
	if( (lbd[1+k]=='f') &&  ( (k*(x[1+k]-bd[1+k])) > 0.0) ) {
	  /* value not within the support */
	  x[1+k] = bd[1+k];             /* jump to boundary instead and ... */
	  step = x[1+k]-x[1];     /* change step size */
	}
	
	/* density at next point */
	y[1+k] = density(x[1+k]);

	/* density o.k. */
	if( y[1+k] == 0.0 )
	  step=step/4.0;  /* reduce stepsize */
	else
	  break;          /* o.k. */
      }
	
      /* ratio of values of the density */
      r = y[1+k] / y[1];

      /* check the ratio */
      if( r <= 1.0 )
	/* f(x[1]) > f(x[1+k])  --> o.k. */
	break;

      /* else: x[1] is not a maximum */
      
      /* --> shift entries */
      x[1-k] = x[1];
      y[1-k] = y[1];
      x[1] = x[1+k];
      y[1] = y[1+k];

      if( r < 4.0 )        /* small change:(heuristical: r<4!)  */
	step=2.0*step;     /* far away from mode -> step is increased by a factor 2 */

      if( r < 1.0002 )
	step=4.0*step;     /* very small change  -> step is increased by a factor 4 */
    }

    if( y[1-k] < 0.0 )      /* we did not go into other direction */
      step=-step;          /* change direction */
    else
      break;
  }

  /* mode known ? */
  if( mode )  /* yes */
    return;

  /* mode unknown --> approximation of the mode */
  while ( (y[1]*y[1]) > (1.1*(y[0]*y[2])) ) {

    /* direction for searching in big intervall */
    k = ( ((x[2]-x[1])-(x[1]-x[0])) < 0.0 ) ? -1 : 1 ;

    /* new point between x[1] and x[1+k] */
    xt = 0.25*x[1+k] + 0.75*x[1];                                           /* new x-value   */
    yt = density(xt);

    if( yt > y[1]) {
      /* new point is nearer to mode */
      x[1-k] = x[1];
      y[1-k] = y[1];
      x[1] = xt;
      y[1] = yt;
    }
    else {
      /* shrink bracket */
      x[1+k] = xt;
      y[1+k] = yt;
    }
  }

} /* end of tdr_find_mode() */

/*-----------------------------------------------------------------*/

static void tdr_find_bracket(double (*density)(double), char *lbd, double *bd, double *x, double *y, double cc)
/* 
   search for proper bracket or mode 
   
   density ... pointer to density function
   lbd[3]  ... status of left (lbd[0]) and right (lbd[2]) bound. 'i' inifinite, 'f' finite.
   bd[3]   ... boundary of interval. bd[0] left, bd[2] right.
   x[3]    ... touching points x[0] <= x[1] <= x[2]. x[1] mode. hu < y[k] < ho  for k=2,2
   y[3]    ... density at x[i]
   cc      ... parameter of the transformation

   hof     ... constant in equation (7) 
*/
{
  int k;                /* direction of search */
  char is_boundary;     /* indicates if point is on boundary */
  double r;             /* ratio of values of density */
  double hof;           /* constant in equation (7) */
  double hl,hu;         /* lower/upper bound for tolerance for the solution of equation (7) */
  double xt, yt;        /* temp values for x and f(x) */

  /* set tolerance */
  hof = ( cc != 0. ) ? 0.25 : exp(-1.);    /*  compare to equation (7). 0 < hof < 0,368... */
  hl = 0.64*hof;        /* lower bound of tolerance, "0.64" is a heuristical value */
  hu = hof/0.64;        /* upper bound of tolerance, "0.64" is a heuristical value */


  for( k=-1; k<=1; k=k+2) {
    /* both directions */

    /* point on boundary ? */
    is_boundary = ( (lbd[k+1]=='f') && (x[1+k]==bd[k+1]) ) ? 't' : 'f' ;

    /* ratio of values of the density */
    r = y[1+k] / y[1];

    if( r < hl ) {     
      /*  y[1+k] < hl * y[1] */

      while( 1 ) {
	/* new point */
	xt = x[1]+(x[1+k]-x[1])/2.;
	yt = density(xt);

	/* check ratio */
	r = yt / y[1];
	if( r >= hl)
	  break;

	/* ratio not ok --> exchange point */
	x[1+k] = xt;
      }

      /* check upper bound */
      if( r <= hu )
	/* o.k. */
	x[1+k] = xt;
      else
	/* not o.k. --> use interpolation */
	x[1+k] = w(xt,x[1+k],x[1],yt,y[1+k],y[1]);

      /* at last we need the value of density at new point */
      y[1+k] = density(x[1+k]);
    }
    else {   /* r >= hl  <=>  y[1+k] > hl * y[1]  --> o.k. */

      if( (r>hu) && (is_boundary=='f') ) {
	/*  y[1+k] > hu * y[1]  and boundary not reached */
	
	while( 1 ) {

	  /* new point */
	  xt = x[1]+2.0*(x[1+k]-x[1]);

	  /* point on boundary ? */
	  is_boundary = ( (lbd[k+1]=='f') && ( (k*(xt-bd[k+1]))>=0. ) ) ? 't' : 'f' ;

	  if( is_boundary == 't' )
	    /* boundary reached -> new point = boundary point */
	    xt = bd[k+1];

	  /* density at new point */
	  yt = density(xt);

	  /*ratio of values */
	  r = yt / y[1];
	  if( !((r>hu) && (is_boundary=='f')) )
	    /* break loop if r <= hu or boundary is reached */
	    break;

	  /* exchange point */
	  x[1+k]=xt;
	}

	if( r > 0. ) {
	  /* f(x) != 0 at boundary */

	  /* check upper bound */
	  if( r >= hl)
	    /* o.k. */
	    x[1+k]=xt;
	  else
	    /* not o.k. --> use interpolation */
	    x[1+k]=w(x[1+k],xt,x[1],y[1+k],yt,y[1]);

	  /* at last we need the value of density at new point */
	  y[1+k]=density(x[1+k]);
	}
      }
    }
  }
} /* end of tdr_find_bracket() */

/*-----------------------------------------------------------------*/
#endif
/*-----------------------------------------------------------------*/