#include "randlib.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);extern void pr_fl_ar(float *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);}void main(int argc,char* argv[])/*     Interactive test for PHRTSD*/{#define mxwh 15L#define MXNCAT 100Lstatic long K1 = 1;static long i,is1,is2,itmp,iwhich,mxint,nperm,nrep,ntot,perm[500],ntry,ncat;static float av,avtr,var,vartr,xmin,xmax,array[1000],param[3],pevt,psum,rtry,             prob[MXNCAT];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.");    puts(" For multinomial asks for number of events to be");    puts("     classified, number of categories in which they");    puts("     are to be classified, and the probabilities that");    puts("     an event will be classified in the categories,");    puts("     for all but the last category.  Prints table of");    puts("     number of events by category, true probability");    puts("     associated with each category, and observed");    puts("     proportion of events in each category.");/*     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");    puts("     (12) Generate gamma deviates");    puts("     (13) Generate multinomial outcomes");    puts("     (14) Generate normal deviates");    puts("     (15) Generate negative binomial outcomes");    scanf("%ld",&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");        exit(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 S320;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 S320;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 S320;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 S320;S100:    if(5 != iwhich) goto S140;/*     Random permutation*/S110:    puts(" Enter size of permutation - range 1 to 500");    scanf("%ld",&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<=nperm; i++) *(perm+i-1) = i;    genprm(perm,nperm);    puts(" Perm Generated");    pr_ln_ar(perm, nperm);    goto S320;S140:    if(6 != iwhich) goto S170;/*     Uniform integer*/    puts(" Enter maximum uniform integer - upper limit 1000");    scanf("%ld",&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("%ld",&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 S320;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 S320;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 S320;S210:    if(9 != iwhich) goto S240;/*     Binomial outcomes*/    strcpy( pr_type, "bin ");    printf(" Enter number of trials, Prob event for binomial outcomes\n");    scanf("%ld%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 S320;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<=1000; i++) *(iarray+i-1) = ignpoi(*param);    for(i=1; i<=1000; i++) *(array+i-1) = *(iarray+i-1);    statistics(array,1000L,&av,&var,&xmin,&xmax);    trstatistics(pr_type,param,&avtr,&vartr);    pr_ans( av, avtr, var, vartr, xmin, xmax );    goto S320;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 S320;S270:    if(12 != iwhich) goto S280;/*     Gamma deviates*/    strcpy( pr_type, "gamm");    puts(" Enter (real) A, (real) R for Gamma deviate");    scanf("%f%f",param,(param+1));    for(i=1; i<=1000; i++) *(array+i-1) = gengam(*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 S320;S280: