#include "ranlib.h"
#include <stdio.h>
#include <math.h>
#include <string.h>
#define	min(a,b) ((a) <= (b) ? (a) : (b))

extern void pr_ans( float av, float avtr, float	var, float vartr,
		    float xmin,	float xmax);
extern void statistics(float *x,long n,float *av,
		 float *var,float *xmin,float *xmax);
extern void trstatistics(char pr_type[4],float *parin,float *av,float *var);
extern void pr_ln_ar(long *array, long larray);

void statistics(float *x,long n,float *av,float	*var,
		 float *xmin,float *xmax)
/*
**********************************************************************
     SUBROUTINE	STATISTICS( X, N, AV, VAR)
	       compute STATistics
			      Function
     Computes AVerage and VARiance of array X(N).
**********************************************************************
*/
{
static long i;
static float sum;

    *xmin = *xmax = *x;
    sum	= 0.0;
    for(i=0; i<n; i++) {
	sum += *(x+i);
	if(*(x+i) < *xmin) *xmin = *(x+i);
	if(*(x+i) > *xmax) *xmax = *(x+i);
    }
    *av	= sum/(float)n;
    sum	= 0.0;
    for(i=0; i<n; i++) sum += pow(*(x+i)-*av,2.0);
    *var = sum/(float)(n-1);
}
int main(int argc,char*	argv[])
/*
     Interactive test for PHRTSD
*/
{
#define	mxwh 11L
static long K1 = 1;
static long i,is1,is2,itmp,iwhich,mxint,nperm,nrep,ntot,perm[500],ntry;
static float av,avtr,var,vartr,xmin,xmax,array[1000],param[3],pevt;
static char pr_type[4],phrase[100];
static long iarray[1000];

    puts(" Tests most generators of specific distributions.");
    puts(" Generates 1000 deviates: reports mean and variance.");
    puts(" Also	reports	theoretical mean and variance.");
    puts(" If theoretical mean or var doesn't exist prints -1.");
    puts(" For permutations, generates one permutation of 1..n");
    puts("     and prints it.");
    puts(" For uniform integers	asks for upper bound, number of");
    puts("     replicates per integer in 1..upper bound.");
    puts("     Prints table of num times each integer generated.");
/*
     Menu for choosing tests
*/
S10:
    puts(" Enter number	corresponding to choice:");
    puts("	(0) Exit this program");
    puts("	(1) Generate Chi-Square	deviates");
    puts("	(2) Generate noncentral	Chi-Square deviates");
    puts("	(3) Generate F deviates");
    puts("	(4) Generate noncentral	F  deviates");
    puts("	(5) Generate random permutation");
    puts("	(6) Generate uniform integers");
    puts("	(7) Generate uniform reals");
    puts("	(8) Generate beta deviates");
    puts("	(9) Generate binomial outcomes");
    puts("     (10) Generate Poisson outcomes");
    puts("     (11) Generate Exponential deviates");
    scanf("%d",&iwhich);
    if(!(iwhich	< 0 || iwhich >	mxwh)) goto S20;
    printf(" Choices are 1..%12ld - try	again.\n",mxwh);
    goto S10;
S20:
    if(iwhich == 0) {
	puts(" Normal termination rn tests");
	return 0;
      }
    puts(" Enter phrase	to initialize rn generator");
    scanf("%s",phrase);
    phrtsd(phrase,&is1,&is2);
    setall(is1,is2);
    if(1 != iwhich) goto S40;
/*
     Chi-square	deviates
*/
    strcpy(pr_type,"chis");
    puts(" Enter (real)	df for the chi-square generation");
    scanf("%f",param);
    for(i=0; i<1000; i++) *(array+i) = genchi(*param);
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av,	avtr, var, vartr, xmin,	xmax );
    goto S280;
S40:
    if(2 != iwhich) goto S60;
/*
     Noncentral	Chi-square deviates
*/
    strcpy( pr_type, "ncch");
    puts(" Enter (real)	df");
    puts("	 (real)	noncentrality parameter");
    scanf("%f%f",param,(param+1));
    for(i=0; i<1000; i++) *(array+i) = gennch(*param,*(param+1));
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av,	avtr, var, vartr, xmin,	xmax );
    goto S280;
S60:
    if(3 != iwhich) goto S80;
/*
     F deviates
*/
    strcpy( pr_type, "f	  ");
    puts(" Enter (real)	df of the numerator");
    puts("	 (real)	df of the denominator");
    scanf("%f%f",param,(param+1));
    for(i=0; i<1000; i++) *(array+i) = genf(*param,*(param+1));
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av,	avtr, var, vartr, xmin,	xmax );
    goto S280;
S80:
    if(4 != iwhich) goto S100;
/*
     Noncentral	F deviates
*/
    strcpy( pr_type, "ncf ");
    puts(" Enter (real)	df of the numerator");
    puts("	 (real)	df of the denominator");
    puts("	 (real)	noncentrality parameter");
    scanf("%f%f%f",param,(param+1),(param+2));
    for(i=0; i<1000; i++) *(array+i) = gennf(*param,*(param+1),*(param+2));
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av,	avtr, var, vartr, xmin,	xmax );
    goto S280;
S100:
    if(5 != iwhich) goto S140;
/*
     Random permutation
*/
S110:
    puts(" Enter size of permutation - range 1 to 500");
    scanf("%d",&nperm);
    if(!(nperm < 1 || nperm > 500)) goto S120;
    printf(" Permutation size must be between 1	and 500	- try again!\n");
    goto S110;
S120:
    printf("	   Random Permutation Generated	- Size%12ld\n",nperm);
    for(i=1; i==20; i++) *(perm+i-1) = 0;
    for(i=1; i<=nperm; i++) *(perm+i-1)	= i;
    genprm(perm,nperm);
    puts(" Perm	Generated");
    pr_ln_ar(perm, nperm);
    goto S280;
S140:
    if(6 != iwhich) goto S170;
/*
     Uniform integer
*/
    puts(" Enter maximum uniform integer - upper limit 1000");
    scanf("%d",&mxint);
    if (mxint <	1 || mxint > 1000){
	puts(" Maximum integer must be between 1 and 1000");
	goto S140;
      }
    puts(" Enter number	of replications	per integer");
    scanf("%d",&nrep);
    for(i=1; i<=1000; i++) *(iarray+i-1) = 0;
    ntot = mxint*nrep;
    for(i=1; i<=ntot; i++) {
	itmp = ignuin(K1,mxint);
	*(iarray+itmp-1) += 1;
    }
    puts("	   Counts of Integers Generated");
    pr_ln_ar(iarray,mxint);
    goto S280;
S170:
    if(7 != iwhich) goto S190;
/*
     Uniform real
*/
    strcpy( pr_type, "unif");
    puts(" Enter Low then High bound for uniforms");
    scanf("%f%f",param,(param+1));
    for(i=1; i<=1000; i++) *(array+i-1)	= genunf(*param,*(param+1));
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av,	avtr, var, vartr, xmin,	xmax );
    goto S280;
S190:
    if(8 != iwhich) goto S210;
/*
     Beta deviate
*/
    strcpy( pr_type, "beta");
    puts(" Enter A, B for Beta deviate");
    scanf("%f%f",param,(param+1));
    for(i=1; i<=1000; i++) *(array+i-1)	= genbet(*param,*(param+1));
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av,	avtr, var, vartr, xmin,	xmax );
    goto S280;
S210:
    if(9 != iwhich) goto S240;
/*
     Binomial outcomes
*/
    strcpy( pr_type, "bin ");
    printf(" Enter number of trials, Prob event	for binomial outcomes\n");
    scanf("%d%f",&ntry,&pevt);
    for(i=1; i<=1000; i++) *(iarray+i-1) = ignbin(ntry,pevt);
    for(i=1; i<=1000; i++) *(array+i-1)	= *(iarray+i-1);
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    *param = ntry;
    *(param+1) = pevt;
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av,	avtr, var, vartr, xmin,	xmax );
    goto S280;
S240:
    if(10 != iwhich) goto S250;
/*
     Poisson outcomes
*/
    strcpy( pr_type, "pois");
    puts(" Enter mean for Poisson generation");
    scanf("%f",param);
    for(i=1; i<=10000; i++) {
      printf("%d\n", ignpoi(*param));
    }
    for(i=1; i<1000; i++) {
      *(iarray+i-1) = ignpoi(*param);
    }
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av,	avtr, var, vartr, xmin,	xmax );
    goto S280;
S250:
    if(11 != iwhich) goto S270;
/*
     Exponential outcomes
*/
    strcpy( pr_type, "exp ");
    puts(" Enter mean for Exponential generation");
    scanf("%f",param);
    for(i=1; i<=1000; i++) *(array+i-1)	= genexp(*param);
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av,	avtr, var, vartr, xmin,	xmax );
    goto S280;
S280:
S270:
    goto S10;
#undef mxwh
}

void trstatistics(char pr_type[4],float	*parin,float *av,float *var)
/*
**********************************************************************
     SUBROUTINE	TRSTATISTICS( TYPE, PARIN, AV, VAR )
	       TRue STATistics
     Returns mean and variance for a number of statistical distribution
     as	a function of their parameters.
			      Arguments
     TYPE --> Character	string indicating type of distribution
	     'chis' chisquare
	     'ncch' noncentral chisquare
	     'f'    F (variance	ratio)
	     'ncf'  noncentral f
	     'unif' uniform
	     'beta' beta distribution
			 CHARACTER*(4) TYPE
     PARIN --> Array containing	parameters of distribution
	      chisquare
	       PARIN(1)	is df
	      noncentral chisquare
	       PARIN(1)	is df
	       PARIN(2)	is noncentrality parameter
	      F	(variance ratio)
	       PARIN(1)	is df numerator
	       PARIN(2)	is df denominator
	      noncentral F
	       PARIN(1)	is df numerator
	       PARIN(2)	is df denominator
	       PARIN(3)	is noncentrality parameter
	      uniform
	       PARIN(1)	is LOW bound
	       PARIN(2)	is HIGH	bound
	      beta
	       PARIN(1)	is A
	       PARIN(2)	is B
			 REAL PARIN(*)
	      binonial
	       PARIN(1)	is Number of trials
	       PARIN(2)	is Prob	Event at Each Trial
	      poisson
	       PARIN(1)	is Mean
	      exponential
	       PARIN(1)	is Reciprocal of Mean
     AV	<-- Mean of specified distribution with	specified parameters
			 REAL AV
     VAR <-- Variance of specified distribution	with specified paramete
			 REAL VAR
			      Note
     AV	and Var	will be	returned -1 if mean or variance	is infinite
*/
{
static float a,b,range;

    if(strcmp("chis",pr_type) !=0 ) goto S10;
    *av	= *parin;
    *var = 2.0**parin;
    goto S170;
S10:
    if(strcmp("ncch",pr_type) != 0) goto S20;
    a =	*parin+*(parin+1);
    b =	*(parin+1)/a;
    *av	= a;
    *var = 2.0*a*(1.0+b);
    goto S170;
S20:
    if(strcmp("f   ",pr_type) != 0) goto S70;
    if(!(*(parin+1) <= 2.0001))	goto S30;
    *av	= -1.0;
    goto S40;
S30:
    *av	= *(parin+1)/(*(parin+1)-2.0);
S40:
    if(!(*(parin+1) <= 4.0001))	goto S50;
    *var = -1.0;
    goto S60;
S50:
    *var = 2.0*pow(*(parin+1),2.0)*(*parin+*(parin+1)-2.0)/(*parin*pow(*(parin
      +1)-2.0,2.0)*(*(parin+1)-4.0));
S60:
    goto S170;
S70:
    if(strcmp("ncf ",pr_type) != 0) goto S120;
    if(!(*(parin+1) <= 2.0001))	goto S80;
    *av	= -1.0;
    goto S90;
S80:
    *av	= *(parin+1)*(*parin+*(parin+2))/((*(parin+1)-2.0)**parin);
S90:
    if(!(*(parin+1) <= 4.0001))	goto S100;
    *var = -1.0;
    goto S110;
S100:
    a =	pow(*parin+*(parin+2),2.0)+(*parin+2.0**(parin+2))*(*(parin+1)-2.0);
    b =	pow(*(parin+1)-2.0,2.0)*(*(parin+1)-4.0);
    *var = 2.0*pow(*(parin+1)/ *parin,2.0)*(a/b);
S110:
    goto S170;
S120:
    if(strcmp("unif",pr_type) != 0) goto S130;
    range = *(parin+1)-*parin;
    *av	= *parin+range/2.0;
    *var = pow(range,2.0)/12.0;
    goto S170;
S130:
    if(strcmp("beta",pr_type) != 0) goto S140;
    *av	= *parin/(*parin+*(parin+1));
    *var = *av**(parin+1)/((*parin+*(parin+1))*(*parin+*(parin+1)+1.0));
    goto S170;
S140:
    if(strcmp("bin ",pr_type) != 0) goto S150;
    *av	= *parin**(parin+1);
    *var = *av*(1.0-*(parin+1));
    goto S170;
S150:
    if(strcmp("pois",pr_type) != 0) goto S155;
    *av	= *var = *parin;
    goto S170;
S155:
    if(strcmp("exp ",pr_type) != 0) goto S160;
    *av	= *parin;
    *var = *av * *av;
    goto S170;
S160:

    puts("Unimplemented	pr_type	in TRSTATISTICS");
    exit(1);
S170:
    return;
}

void pr_ans( float av, float avtr, float var, float vartr,
	     float xmin, float xmax)
{

    printf(" Mean Generated	     %15.7g	True	      %15.7g\n",
	   (double)av, (double)avtr);
    printf(" Variance Gener	     %15.7g	True	      %15.7g\n",
	   (double)var,	(double)vartr);
    printf(" Minimum		     %15.7g	Maximum	      %15.7g\n",
	   (double)xmin, (double)xmax);
}

void pr_ln_ar(long *array, long	larray)
/*
  Print	array of longs;	length of array	is larray.  Print 10 to	a
  line using 7 characters for each field.
*/
{
    int	lo, hi,	i, up;
    up = larray-1;
    for	(lo=0; lo <=up;	lo +=10) {
	hi = min(lo+9,up);
	for (i=lo; i<=hi; i++)
	    printf(" %7ld",*(array+i));
    puts("");
    }
}
