/* PURPOSE:
 * This program generates a synthetic workload (i.e. a stream of jobs)
 *    that satisfies a Hyper-Erlang distribution(s) of Common Order.
 *
 * The input parameters used in this program were derived from the observed
 *   workload on a 322-node partition of the Cornell Theory Center (CTC) SP2,
 *   over the period from June 25, 1996 to September 12, 1996.
 *
 * Detailed description of the models can be found in the paper titled:
 *            "Modeling of Workload in MPPs"
 * by Jann J. et al, published in the book ISBN 3-540-63574-2 by Springer Verlag,
 * 1997, Lecture Notes in Computer Science vol. 1291, pp. 95-116;
 * Job Scheduling Stategies for Parallel Processing.
 *
 * INPUT DATA:
 * The model is actually composed of 10 sub-models, one for each range between
 * successive powers of 2 in the total range of 1 to 322.
 * For each range, there is a separate model for CPU time and interarrival times.
 * The parameters of each such range are:
 * 1. Fraction of jobs desired (fjobs)
 * 2. Minimum & maximum number of processes per job (pmin & pmax)
 * 3. Inter-arrival time: lambda1, lambda2, order of the HErlang distrib, p1(=prob)
 * 4. Service/CPU time:   lambda1, lambda2, order of the HErlang distrib, p1(=prob)
 * In addition there is a parameter that specifies the total number of jobs to
 * generate (njobs).
 *
 * OUTPUT:
 * For each job in the workload:
 * 1. Number of processors for the job
 * 2. Total CPU time used by the job (=service time)
 * 3. Arrival time of the job
 *
 * For the whole workload of njobs :
 * 1. 1st, 2nd & 3rd moments of the inter-arrival-time        of the jobs
 * 2. 1st, 2nd & 3rd moments of the service-time (= CPU-time) of the jobs
 *
 * NOTES:
 * To keep this program simple, we have eliminated all checkings of the input.
 * If you use the program often,  we strongly recommend you to add code to
 *     catch error conditions.
 *
 * AUTHOR: Joefon Jann joefon@watson.ibm.com
 *         IBM T. J. Watson Research Center,
 *         P. O. Box 218, Yorktown Heights, NY 10598, USA
 * DATE:   11/1997
 *
 * Extended by Dror Feitelson, Hebrew University, 3/1999.
 * Extensions include parameters for the full model (with 10 submodels),
 * selection of sub-models according to the distribution of job sizes,
 * integration of the job streams from the different sub-models,
 * and output using the standard workload format.
 * Extended again to include ASCI parameters, 7/2001.
 * Corrected shameful typos in ASCI parameters, 9/2004.
 * Also discovered that ASCI params seem to be erroneous anyway.
 */

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

/*******************************************************
 * These are things you may want to change
 *
 * The following dictates which model to use.
 * the default is to use CTC.
 * defining MODEL_ASCI replaces this with the ASCI data.
 *
 * NOTE: it seems that the ASCI data has an error.
 * It leads to extremely short runtimes and low utilization.
 *
#define MODEL_ASCI
 */
/*
 * The following causes output in the standard format.
 * The alternative is a line describing each job.
 */
#define STD_FORMAT
/*
 * The number of jobs created is defined here.
 */
#define NUM_JOBS 100000
/*
 *******************************************************/

#define	MAXMODELS 10

#define min(a,b) (((a) < (b)) ? (a) : (b))
#define RANSTREAM 4
#define RANPOOLSIZE 10000

typedef struct
{  double p1, lambda1, lambda2, ur;
   int    order;
} MODEL_DATA;

typedef struct
{  double fjobs;
   int    pmin, pmax;
   int    count_jobs;
} SIZE_DATA;

/* Global Variables: */
double ran_pool[RANSTREAM][RANPOOLSIZE];
int    ran_remain[RANSTREAM];
int    nmodels;

/* - - - - - Function prototypes: - - - - - - - - - */
int find_next_arrival(double *next_arrival);

double ur_to_HER_r(double ur, MODEL_DATA model);

double H_Er_cdf(double x, void *model);

void solve_bisec( double xleft_0,
                  double xright_0,
                  double tol,
                  double *x,
                  double *y,
                  void *model_data,
                  double (*f)(double, void *),
                  int *rc );

double fn_to_solve(double x, void *);

int    fac(int k);

double str_rand(int str);

/* =================================================== */
main(int argc, char **argv)
{  MODEL_DATA inter_arrival_model[MAXMODELS], cpu_model[MAXMODELS];
   double next_arrival[MAXMODELS];
   double inter_arrival_mu1, inter_arrival_mu2,
          inter_arrival_mu3;
   double cpu_mu1, cpu_mu2, cpu_mu3;
   double simtime, inter_arrival_time, cpu_time_used;
   double dtemp;

   SIZE_DATA job_size_model[MAXMODELS];
   int pmin, pmax, p_needed;
   int imodel;
   int njobs, ijob;
   int i,j;

   int maxnodes;
   char sysname[10];

   /* Beginning of the INPUT section - - - - - - - - - - - - - - - - - - */

   /* The number of jobs generated is given by njobs
    */
   njobs = NUM_JOBS;

#ifndef MODEL_ASCI
   /* this is for CTC */

   maxnodes = 322;
   sprintf(sysname, "CTC");

   /* The number of sub-models (nmodels) reflects the number of
    * lines from Tables 1 and 3 in the Results
    * section of the paper concurrently used for
    * generating the job stream.
    */
   nmodels = 10;

   /* The following section initializes the parameters of the sub-models.
    * The values are taken from tables 1 and 3 in the paper.
    */
   inter_arrival_model[0].lambda1 = 1.43e-04;
   inter_arrival_model[0].lambda2 = 2.05e-03;
   inter_arrival_model[0].order   = 1;
   inter_arrival_model[0].p1      = 8.56e-02;

   inter_arrival_model[1].lambda1 = 3.31e-05;
   inter_arrival_model[1].lambda2 = 9.94e-04;
   inter_arrival_model[1].order   = 1;
   inter_arrival_model[1].p1      = 1.67e-01;

   inter_arrival_model[2].lambda1 = 5.58e-05;
   inter_arrival_model[2].lambda2 = 1.13e-03;
   inter_arrival_model[2].order   = 1;
   inter_arrival_model[2].p1      = 1.55e-01;

   inter_arrival_model[3].lambda1 = 7.37e-05;
   inter_arrival_model[3].lambda2 = 1.88e-03;
   inter_arrival_model[3].order   = 1;
   inter_arrival_model[3].p1      = 2.38e-01;

   inter_arrival_model[4].lambda1 = 6.87e-05;
   inter_arrival_model[4].lambda2 = 7.16e-04;
   inter_arrival_model[4].order   = 1;
   inter_arrival_model[4].p1      = 1.30e-01;

   inter_arrival_model[5].lambda1 = 4.08e-05;
   inter_arrival_model[5].lambda2 = 4.98e-04;
   inter_arrival_model[5].order   = 1;
   inter_arrival_model[5].p1      = 1.49e-01;

   inter_arrival_model[6].lambda1 = 5.67e-05;
   inter_arrival_model[6].lambda2 = 4.91e-03;
   inter_arrival_model[6].order   = 1;
   inter_arrival_model[6].p1      = 4.79e-01;

   inter_arrival_model[7].lambda1 = 3.33e-05;
   inter_arrival_model[7].lambda2 = 6.07e-04;
   inter_arrival_model[7].order   = 2;
   inter_arrival_model[7].p1      = 5.03e-01;

   inter_arrival_model[8].lambda1 = 4.48e-06;
   inter_arrival_model[8].lambda2 = 2.78e-04;
   inter_arrival_model[8].order   = 1;
   inter_arrival_model[8].p1      = 2.78e-01;

   inter_arrival_model[9].lambda1 = 1.69e-06;
   inter_arrival_model[9].lambda2 = 1.99e-05;
   inter_arrival_model[9].order   = 1;
   inter_arrival_model[9].p1      = 1.18e-01;

   cpu_model[0].lambda1 = 1.23e-04;
   cpu_model[0].lambda2 = 4.72e-03;
   cpu_model[0].order   = 4;
   cpu_model[0].p1      = 4.60e-01;

   cpu_model[1].lambda1 = 8.04e-05;
   cpu_model[1].lambda2 = 1.02e-01;
   cpu_model[1].order   = 7;
   cpu_model[1].p1      = 3.94e-01;

   cpu_model[2].lambda1 = 3.23e-05;
   cpu_model[2].lambda2 = 1.61e-02;
   cpu_model[2].order   = 4;
   cpu_model[2].p1      = 2.51e-01;

   cpu_model[3].lambda1 = 1.17e-05;
   cpu_model[3].lambda2 = 9.18e-03;
   cpu_model[3].order   = 1;
   cpu_model[3].p1      = 3.76e-01;

   cpu_model[4].lambda1 = 7.65e-06;
   cpu_model[4].lambda2 = 6.40e-04;
   cpu_model[4].order   = 3;
   cpu_model[4].p1      = 2.87e-01;

   cpu_model[5].lambda1 = 2.80e-06;
   cpu_model[5].lambda2 = 2.51e-04;
   cpu_model[5].order   = 1;
   cpu_model[5].p1      = 3.75e-01;

   cpu_model[6].lambda1 = 1.59e-06;
   cpu_model[6].lambda2 = 7.00e-05;
   cpu_model[6].order   = 3;
   cpu_model[6].p1      = 2.78e-01;

   cpu_model[7].lambda1 = 1.19e-06;
   cpu_model[7].lambda2 = 3.99e-04;
   cpu_model[7].order   = 6;
   cpu_model[7].p1      = 2.07e-01;

   cpu_model[8].lambda1 = 1.62e-05;
   cpu_model[8].lambda2 = 3.15e-03;
   cpu_model[8].order   = 3;
   cpu_model[8].p1      = 1.48e-01;

   cpu_model[9].lambda1 = 5.18e-07;
   cpu_model[9].lambda2 = 4.51e-04;
   cpu_model[9].order   = 6;
   cpu_model[9].p1      = 1.65e-01;

   job_size_model[0].fjobs = 0.4042;
   job_size_model[0].pmin  = 1;
   job_size_model[0].pmax  = 1;
   job_size_model[0].count_jobs = 0;
   
   job_size_model[1].fjobs = 0.0717;
   job_size_model[1].pmin  = 2;
   job_size_model[1].pmax  = 2;
   job_size_model[1].count_jobs = 0;
      
   job_size_model[2].fjobs = 0.1197;
   job_size_model[2].pmin  = 3;
   job_size_model[2].pmax  = 4;
   job_size_model[2].count_jobs = 0;

   job_size_model[3].fjobs = 0.1164;
   job_size_model[3].pmin  = 5;
   job_size_model[3].pmax  = 8;
   job_size_model[3].count_jobs = 0;

   job_size_model[4].fjobs = 0.1361;
   job_size_model[4].pmin  = 9;
   job_size_model[4].pmax  = 16;
   job_size_model[4].count_jobs = 0;

   job_size_model[5].fjobs = 0.0789;
   job_size_model[5].pmin  = 17;
   job_size_model[5].pmax  = 32;
   job_size_model[5].count_jobs = 0;

   job_size_model[6].fjobs = 0.0491;
   job_size_model[6].pmin  = 33;
   job_size_model[6].pmax  = 64;
   job_size_model[6].count_jobs = 0;

   job_size_model[7].fjobs = 0.0132;
   job_size_model[7].pmin  = 65;
   job_size_model[7].pmax  = 128;
   job_size_model[7].count_jobs = 0;

   job_size_model[8].fjobs = 0.0064;
   job_size_model[8].pmin  = 129;
   job_size_model[8].pmax  = 256;
   job_size_model[8].count_jobs = 0;

   job_size_model[9].fjobs = 0.0037;
   job_size_model[9].pmin  = 257;
   job_size_model[9].pmax  = 322;
   job_size_model[9].count_jobs = 0;
#else
   /* this is for ASCI */

   maxnodes = 256;
   sprintf(sysname, "ASCI");

   nmodels = 6;

   inter_arrival_model[0].lambda1 = 1.02e-04;
   inter_arrival_model[0].lambda2 = 2.06e-03;
   inter_arrival_model[0].order   = 1;
   inter_arrival_model[0].p1      = 1.07e-01;

   inter_arrival_model[1].lambda1 = 1.69e-04;
   inter_arrival_model[1].lambda2 = 1.40e-03;
   inter_arrival_model[1].order   = 1;
   inter_arrival_model[1].p1      = 4.10e-01;

   inter_arrival_model[2].lambda1 = 1.94e-04;
   inter_arrival_model[2].lambda2 = 3.02e-03;
   inter_arrival_model[2].order   = 2;
   inter_arrival_model[2].p1      = 3.09e-01;

   inter_arrival_model[3].lambda1 = 3.17e-04;
   inter_arrival_model[3].lambda2 = 2.15e-03;
   inter_arrival_model[3].order   = 1;
   inter_arrival_model[3].p1      = 7.71e-01;

   inter_arrival_model[4].lambda1 = 8.94e-05;
   inter_arrival_model[4].lambda2 = 4.56e-03;
   inter_arrival_model[4].order   = 3;
   inter_arrival_model[4].p1      = 3.07e-01;

   inter_arrival_model[5].lambda1 = 7.94e-05;
   inter_arrival_model[5].lambda2 = 1.71e-03;
   inter_arrival_model[5].order   = 2;
   inter_arrival_model[5].p1      = 4.26e-01;

   cpu_model[0].lambda1 = 5.20e-04;
   cpu_model[0].lambda2 = 4.65e-03;
   cpu_model[0].order   = 1;
   cpu_model[0].p1      = 2.85e-01;

   cpu_model[1].lambda1 = 1.15e-03;
   cpu_model[1].lambda2 = 6.08e-02;
   cpu_model[1].order   = 5;
   cpu_model[1].p1      = 5.83e-01;

   cpu_model[2].lambda1 = 1.32e-04;
   cpu_model[2].lambda2 = 1.50e-01;
   cpu_model[2].order   = 1;
   cpu_model[2].p1      = 3.03e-01;

   cpu_model[3].lambda1 = 4.85e-04;
   cpu_model[3].lambda2 = 1.09e-02;
   cpu_model[3].order   = 2;
   cpu_model[3].p1      = 9.00e-01;

   cpu_model[4].lambda1 = 8.25e-04;
   cpu_model[4].lambda2 = 4.71e-02;
   cpu_model[4].order   = 3;
   cpu_model[4].p1      = 5.38e-01;

   cpu_model[5].lambda1 = 4.52e-04;
   cpu_model[5].lambda2 = 8.61e-03;
   cpu_model[5].order   = 3;
   cpu_model[5].p1      = 4.37e-01;

   job_size_model[0].fjobs = 0.3718;
   job_size_model[0].pmin  = 1;
   job_size_model[0].pmax  = 8;
   job_size_model[0].count_jobs = 0;
   
   job_size_model[1].fjobs = 0.1963;
   job_size_model[1].pmin  = 9;
   job_size_model[1].pmax  = 16;
   job_size_model[1].count_jobs = 0;
      
   job_size_model[2].fjobs = 0.1386;
   job_size_model[2].pmin  = 17;
   job_size_model[2].pmax  = 32;
   job_size_model[2].count_jobs = 0;

   job_size_model[3].fjobs = 0.2148;
   job_size_model[3].pmin  = 33;
   job_size_model[3].pmax  = 64;
   job_size_model[3].count_jobs = 0;

   job_size_model[4].fjobs = 0.0393;
   job_size_model[4].pmin  = 65;
   job_size_model[4].pmax  = 128;
   job_size_model[4].count_jobs = 0;

   job_size_model[5].fjobs = 0.0393;
   job_size_model[5].pmin  = 129;
   job_size_model[5].pmax  = 256;
   job_size_model[5].count_jobs = 0;
#endif
   /* End of the INPUT section - - - - - - - - - - - - - - - - - - */
   
#ifdef STD_FORMAT
   /* print header comments */
   printf("; Version: 2\n");
   printf("; Information: Jann'97 Workload Model for %s\n",sysname);
   printf(";              From the Parallel Workloads Archive\n");
   printf(";              http://www.cs.huji.ac.il/labs/parallel/workload\n");
   printf("; Acknowledge: Joefon Jann, joefon@watson.ibm.com\n");
   printf("; Conversion: Dror Feitelson, feit@cs.huji.ac.il\n");
   printf("; MaxNodes: %d\n;\n",maxnodes);
#endif

   inter_arrival_mu1=0;
   inter_arrival_mu2=0;
   inter_arrival_mu3=0;

   cpu_mu1=0;
   cpu_mu2=0;
   cpu_mu3=0;

   /* Each sub-model includes interarrival times from the previous job in
    * the same class. We need to interleave the 10 streams according to
    * arrival times, so we need to start with a first sample of the arrival
    * time from each sub model.
    */
   for (i=0 ; i<nmodels ; i++)
   {  dtemp=str_rand(0);
      next_arrival[i] = ur_to_HER_r(dtemp,inter_arrival_model[i]);
   }
   /* Now find the minimal one, which is the first arrival, and make it
    * time 0.
    */
   simtime = next_arrival[ find_next_arrival( next_arrival ) ];
   for (i=0 ; i<nmodels ; i++)
     next_arrival[i] -= simtime;
   simtime = 0;

   /* Here we are using the number of jobs (njobs) as the termination criterion.
    * One may choose to use simtime (i.e. total length of time represented
    *                              in this workload) as a termination criterion.
    * In that case, use simtime in the following WHILE statement.
    */
   for( ijob=0 ; ijob<njobs ; ijob++ )
   {  /* First find the sub-model with the next arrival
       */
     inter_arrival_time = simtime;
     simtime = next_arrival[ imodel=find_next_arrival(next_arrival) ];
     inter_arrival_time = simtime - inter_arrival_time;

     /* Tabulate and prepare next arrival time for this submodel
      */
     job_size_model[imodel].count_jobs++;
     dtemp=str_rand(0);
     next_arrival[imodel] = simtime + 
			    ur_to_HER_r(dtemp,inter_arrival_model[imodel]);

      /* The number of processors is uniform in [pmin,pmax] for
       * the selected sub-model
       */
      pmin = job_size_model[imodel].pmin;
      pmax = job_size_model[imodel].pmax;
      dtemp=str_rand(2);
      dtemp *= (pmax-pmin+1);
      j=dtemp;
      if( ((double) j) < dtemp ) { p_needed = j + pmin; }
      else                       { p_needed = j + pmin - 1; }

      /* Get CPU time according to the sub-model
       */
      dtemp=str_rand(1);
      cpu_time_used = ur_to_HER_r(dtemp,cpu_model[imodel]);

      if (cpu_time_used == 0.0)
      {  ijob--;
	 continue;
      }

      inter_arrival_mu1 +=  inter_arrival_time;
      inter_arrival_mu2 += (inter_arrival_time *
			    inter_arrival_time);
      inter_arrival_mu3 += (inter_arrival_time *
			    inter_arrival_time *
			    inter_arrival_time);

      cpu_mu1 +=  cpu_time_used;
      cpu_mu2 += (cpu_time_used *
		  cpu_time_used);
      cpu_mu3 += (cpu_time_used *
		  cpu_time_used *
		  cpu_time_used);

#ifdef STD_FORMAT
      /* Now output the data for this job, using the standard format.
       * the fields are:
       * job number = ijob+1 [start from 1, not from 0]
       * submit time = simtime
       * wait time = -1 [queueing is a scheduler's problem, not
       *		       part of the workload model]
       * run time = (cpu_time_used / p_needed)
       *	[cpu_time_used is cummulative on all processors, so the
       *	 running time is this divided by the number of processors]
       * used processors = p_needed
       * average CPU time = (cpu_time_used / p_needed) [same as runtime]
       * used memory = -1 [not modeled]
       * requested processors = -1 [not modeled]
       * requested time = -1 [not modeled]
       * requested memory = -1 [not modeled]
       * completed = 1 [ all jobs complete OK in the model ]
       * user ID = -1 [not modeled]
       * group ID = -1 [not modeled]
       * application = -1  [not modeled]
       * queue = -1  [not modeled]
       * partition = -1  [not modeled]
       * preceding job = -1 [dependencies not modeled]
       * think time = -1 [not modeled]
       */
      printf("%5d %12.4f -1 %12.4f %3d %12.4f -1 -1 -1 -1 1 -1 -1 -1 -1 -1 -1 -1\n",
	     ijob+1,
	     simtime,
	     (cpu_time_used / p_needed),
	     p_needed,
	     (cpu_time_used / p_needed)
	    );
#else
      printf( "Job %d used %d CPUs and %e cummulative CPUtime",
	      ijob, p_needed, cpu_time_used);
      printf( " arrived at time %e \n", simtime);
#endif
   }  /* end of for loop */

   printf(";\n;\n;\n; Model checks:\n");
   inter_arrival_mu1 /= (double)ijob;
   inter_arrival_mu2 /= (double)ijob;
   inter_arrival_mu3 /= (double)ijob;
   printf("; Inter-arrival Time Moments of the jobs are:\n");
   printf(";    mu1 = %e %s\n", inter_arrival_mu1,
	  (sysname[0]=='C')? "(should be 4.234e02)" : "");
   printf(";    mu2 = %e %s\n", inter_arrival_mu2,
	  (sysname[0]=='C')? "(should be 1.794e06)" : "");
   printf(";    mu3 = %e %s\n", inter_arrival_mu3,
	  (sysname[0]=='C')? "(should be 2.449e10)" : "");

   cpu_mu1 /= (double)ijob;
   cpu_mu2 /= (double)ijob;
   cpu_mu3 /= (double)ijob;
   printf("; CPU Time Moments of the jobs are:\n");
   printf(";    mu1 = %e %s\n", cpu_mu1,
	  (sysname[0]=='C')? "(should be 9.126e04)" : "");
   printf(";    mu2 = %e %s\n", cpu_mu2,
	  (sysname[0]=='C')? "(should be 2.601e11)" : "");
   printf(";    mu3 = %e %s\n", cpu_mu3,
	  (sysname[0]=='C')? "(should be 2.255e18)" : "");

   for (imodel=0 ; imodel<nmodels ; imodel++)
      printf("; Got %d jobs in range %d-%d (should be %d)\n",
	     job_size_model[imodel].count_jobs,
	     job_size_model[imodel].pmin, job_size_model[imodel].pmax,
	     (int)(job_size_model[imodel].fjobs*(double)ijob) );
}   /* end of main program */

/* ================================================== */
int find_next_arrival(double *next_arrival)
{  int    i, model;
   double min;

   min = next_arrival[0];
   model = 0;
   for (i=1 ; i<nmodels ; i++)
     if (next_arrival[i] < min)
     {  min = next_arrival[i];
        model = i;
     }
   return model;
}

/* ================================================== */
double ur_to_HER_r(double ur, MODEL_DATA m)
{  double xleft, xright;
   double x, y;
   int    rc;

   m.ur = ur;
   xleft = 1.e-20;        /* a small +ve number */
   xright = -log(1.0-ur) / min(m.lambda1,m.lambda2);
   xright *= 100;

   solve_bisec(xleft, xright,
               (double) 1.e-12,
               &x, &y, &m,
               fn_to_solve, &rc);
   if(rc !=0 )
   {  printf(";Unable to find the Random Number for ur = %e \n",m.ur);
      return(0.0);
   }
   return(x);
}

/* =================================================== */
double fn_to_solve(double x, void *model_data)
{  MODEL_DATA *m;

   m = (MODEL_DATA *) model_data;
   return( H_Er_cdf(x, model_data)  -  m->ur );
}

/* =================================================== */
double H_Er_cdf(double x, void *model_data)
{  MODEL_DATA *p;
   double cdf, t1, t2, dtemp;
   int k;
   t1=0;
   t2=0;
   p = (MODEL_DATA *) model_data;

   for(k=0;  k< p->order;  k++)
   {  dtemp=fac(k);
      t1 += (pow(x*p->lambda1,k)/dtemp);
      t2 += (pow(x*p->lambda2,k)/dtemp);
   }
   cdf = 1.0 -    p->p1 *exp(-(p->lambda1*x))*t1
             - (1-p->p1)*exp(-(p->lambda2*x))*t2;
   return(cdf);
 }

/* =================================================== */
void solve_bisec(double xleft_0,
                 double xright_0,
                 double tol,
                 double *x,
                 double *y,
                 void *model_data,
                 double (*f)(double, void *),
                 int *rc)

{  double xleft, yleft,
          xright, yright,
          xnext,ynext;
   int    ic;

   xleft=xleft_0; xright=xright_0;
   yleft=(*f)(xleft, model_data);
   yright=(*f)(xright,model_data);
   *rc=1; ic=0;

   if( yleft*yright  >  0.0 ) { return; }
   /* No solution in this interval */

   while(ic < 100000 )
   {  ic++;
      xnext=(xleft+xright)/2.0;
      ynext=(*f)(xnext,model_data);
      if( ynext*yright  >  0.0 )
      {   xright=xnext;
          yright=ynext;
      }
      else
      {   xleft=xnext;
          yleft=ynext;
      }
      if( fabs(ynext)  <  tol )
      {   *x=xnext; *y=ynext; *rc=0; return;}
   }
}

/* =================================================== */
int fac(int k)
{  int i,j;
   j=1;

   if(k==0) return (j);
   for(i=1;i<=k;i++) j*=i;
   return(j);
}

/* =================================================== */
double str_rand(int strm)
{  int i;
   if(  ran_remain[strm] <=  0 )
   {  for(i=0;i<RANPOOLSIZE;i++)
          ran_pool[strm][i]=drand48();
      ran_remain[strm]=RANPOOLSIZE;
   }
   ran_remain[strm]--;

   return( ran_pool[strm] [RANPOOLSIZE-ran_remain[strm]]) ;
}