/****************************************************************/
/*                                                              */
/*           FILE: stat_dist.c                                  */
/*           WRITTEN BY: Jonathan G. Fiscus                     */
/*           DATE: April 14 1989                                */
/*                 NATIONAL INSTITUTE OF STANDARDS              */
/*                         AND TECHNOLOGY                       */
/*                    SPEECH RECOGNITION GROUP                  */
/*                                                              */
/*           USAGE: This file initializes the distribution      */
/*                  structures for Z1tail, Z2tail, and X2       */
/*                  distributions                               */
/*                                                              */
/****************************************************************/
#include "sctk.h"

#define MAX_Z	10
#define MAX_X2	32

Z_STRUCT Z2tail[MAX_Z] = {
   {  2.575,      "99%",      "01%",      0.99   },
   {  2.329,      "98%",      "02%",      0.98   },
   {  2.170,      "97%",      "03%",      0.97   },
   {  2.055,      "96%",      "04%",      0.96   },
   {  1.960,      "95%",      "05%",      0.95   },
   {  1.881,      "94%",      "06%",      0.94   },
   {  1.811,      "93%",      "07%",      0.93   },                 
   {  1.751,      "92%",      "08%",      0.92   },
   {  1.695,      "91%",      "09%",      0.91   },
   {  1.645,      "90%",      "10%",      0.90   }
};


Z_STRUCT Z1tail[MAX_Z] = {
   {  2.330,      "99%",      "01%",      0.99   },
   {  2.055,      "98%",      "02%",      0.98   },
   {  1.881,      "97%",      "03%",      0.97   },
   {  1.751,      "96%",      "04%",      0.96   },
   {  1.645,      "95%",      "05%",      0.95   },
   {  1.555,      "94%",      "06%",      0.94   },
   {  1.474,      "93%",      "07%",      0.93   },
   {  1.405,      "92%",      "08%",      0.92   },
   {  1.340,      "91%",      "09%",      0.91   },
   {  1.381,      "90%",      "10%",      0.90   }
};


X2_STRUCT X2 = {
    {  99.0,  98.0,  95.0,  90.0,  80.0,  70.0,  50.0,  30.0,  20.0,  10.0,
	     5.0,   2.0,   1.0,   0.0999999  },
    { "99%", "98%", "95%", "90%", "80%", "70%", "50%", "30%", "20%", "10%",
	    "5%",  "2%",  "1%",   ".1%" },
    {  "1%",  "2%",  "5%", "10%", "20%", "30%", "50%", "70%", "80%", "90%", 
	   "95%", "98%", "99%", "99.9%" },
    {{ "df1",
       { 0.00016,
         0.00063,
         0.0039,
         0.016,
         0.064,
         0.15,
         0.46,
         1.07,
         1.64,
         2.71,
         3.84,
         5.41,
         6.64,
         10.83
       }
     },
     { "df2",
       { 0.02,
         0.04,
         0.10,
         0.21,
         0.45,
         0.71,
         1.39,
         2.41,
         3.22,
         4.60,
         5.99,
         7.82,
         9.21,
         13.82
       }
     },
     { "df3",
       { 0.12,
         0.18,
         0.35,
         0.58,
         1.00,
         1.42,
         2.37,
         3.66,
         4.64,
         6.25,
         7.82,
         9.84,
         11.34,
         16.27
       }
     },
     { "df4",
       { 0.30,
         0.43,
         0.71,
         1.06,
         1.65,
         2.20,
         3.36,
         4.88,
         5.99,
         7.78,
         9.49,
         11.67,
         13.28,
         18.46
       }
     },
     { "df5",
       { 0.55,
         0.75,
         1.14,
         1.61,
         2.34,
         3.00,
         4.35,
         6.06,
         7.29,
         9.24,
         11.07,
         13.39,
         15.09,
         20.52
       }
     },
     { "df6",
       { 0.87,
         1.13,
         1.64,
         2.20,
         3.07,
         3.83,
         5.35,
         7.23,
         8.56,
         10.64,
         12.59,
         15.03,
         16.81,
         22.46
       }
     },
     { "df7",
       { 1.24,
         1.56,
         2.17,
         2.83,
         3.82,
         4.67,
         6.35,
         8.38,
         9.80,
         12.02,
         14.07,
         16.62,
         18.48,
         24.32
       }
     },
     { "df8",
       { 1.65,
         2.03,
         2.73,
         3.49,
         4.59,
         5.53,
         7.34,
         9.52,
         11.03,
         13.36,
         15.51,
         18.17,
         20.09,
         26.12
       }
     },
     { "df9",
       { 2.09,
         2.53,
         3.32,
         4.17,
         5.38,
         6.39,
         8.34,
         10.66,
         12.24,
         14.68,
         16.92,
         19.68,
         21.67,
         27.88
       }
     },
     { "df10",
       { 2.56,
         3.06,
         3.94,
         4.86,
         6.18,
         7.27,
         9.34,
         11.78,
         13.44,
         15.99,
         18.31,
         21.16,
         23.21,
         29.59
       }
     },
     { "df11",
       { 3.05,
         3.61,
         4.58,
         5.58,
         6.99,
         8.15,
         10.34,
         12.90,
         14.63,
         17.28,
         19.68,
         22.62,
         24.72,
         31.26
       }
     },
     { "df12",
       { 3.57,
         4.18,
         5.23,
         6.30,
         7.81,
         9.03,
         11.34,
         14.01,
         15.81,
         18.55,
         21.03,
         24.05,
         26.22,
         32.91
       }
     },
     { "df13",
       { 4.11,
         4.76,
         5.89,
         7.04,
         8.63,
         9.93,
         12.34,
         15.12,
         16.98,
         19.81,
         22.36,
         25.47,
         27.69,
         34.53
       }
     },
     { "df14",
       { 4.66,
         5.37,
         6.57,
         7.79,
         9.47,
         10.82,
         13.34,
         16.22,
         18.15,
         21.06,
         23.68,
         26.87,
         29.14,
         36.12
       }
     },
     { "df15",
       { 5.23,
         5.98,
         7.26,
         8.55,
         10.31,
         11.72,
         14.34,
         17.32,
         19.31,
         22.31,
         25.00,
         28.26,
         30.58,
         37.70
       }
     },
     { "df16",
       { 5.81,
         6.61,
         7.96,
         9.31,
         11.15,
         12.62,
         15.34,
         18.42,
         20.46,
         23.54,
         26.30,
         29.63,
         32.00,
         39.29
       }
     },
     { "df17",
       { 6.41,
         7.26,
         8.67,
         10.08,
         12.00,
         13.53,
         16.34,
         19.51,
         21.62,
         24.77,
         27.59,
         31.00,
         33.41,
         40.75
       }
     },
     { "df18",
       { 7.02,
         7.91,
         9.39,
         10.86,
         12.86,
         14.44,
         17.34,
         20.60,
         22.76,
         25.99,
         28.87,
         32.35,
         34.80,
         42.31
       }
     },
     { "df19",
       { 7.63,
         8.57,
         10.12,
         11.65,
         13.72,
         15.35,
         18.34,
         21.69,
         23.90,
         27.20,
         30.14,
         33.69,
         36.19,
         43.82
       }
     },
     { "df20",
       { 8.26,
         9.24,
         10.85,
         12.44,
         14.58,
         16.27,
         19.34,
         22.78,
         25.04,
         28.41,
         31.41,
         35.02,
         37.57,
         45.32
       }
     },
     { "df21",
       { 8.90,
         9.92,
         11.59,
         13.24,
         15.44,
         17.18,
         20.34,
         23.86,
         26.17,
         29.62,
         32.67,
         36.34,
         38.93,
         46.80
       }
     },
     { "df22",
       { 9.54,
         10.60,
         12.34,
         14.04,
         16.31,
         18.10,
         21.24,
         24.94,
         27.30,
         30.81,
         33.92,
         37.66,
         40.29,
         48.27
       }
     },
     { "df23",
       { 10.20,
         11.29,
         13.09,
         14.85,
         17.19,
         19.02,
         22.34,
         26.02,
         28.43,
         32.01,
         35.17,
         38.97,
         41.64,
         49.73
       }
     },
     { "df24",
       { 10.86,
         11.99,
         13.85,
         15.66,
         18.06,
         19.94,
         23.34,
         27.10,
         29.55,
         33.20,
         36.42,
         40.27,
         42.98,
         51.18
       }
     },
     { "df25",
       { 11.52,
         12.70,
         14.61,
         16.47,
         18.94,
         20.87,
         24.34,
         28.17,
         30.68,
         34.38,
         37.65,
         41.57,
         44.31,
         52.62
       }
     },
     { "df26",
       { 12.20,
         13.41,
         15.38,
         17.29,
         19.82,
         21.79,
         25.34,
         29.25,
         31.80,
         35.56,
         38.88,
         42.86,
         45.64,
         54.05
       }
     },
     { "df27",
       { 12.88,
         14.12,
         16.15,
         18.11,
         20.70,
         22.72,
         26.34,
         30.32,
         32.91,
         36.74,
         40.11,
         44.14,
         46.96,
         55.48
       }
     },
     { "df28",
       { 13.56,
         14.85,
         16.93,
         18.94,
         21.59,
         23.65,
         27.34,
         31.39,
         34.03,
         37.92,
         41.34,
         45.42,
         48.28,
         56.89
       }
     },
     { "df29",
       { 14.26,
         15.57,
         17.71,
         19.77,
         22.48,
         24.58,
         28.34,
         32.46,
         35.14,
         39.09,
         42.56,
         46.69,
         49.59,
         58.30
       }
     },
     { "df30",
       { 14.95,
         16.31,
         18.49,
         20.60,
         23.36,
         25.51,
         29.34,
         33.53,
         36.25,
         40.26,
         43.77,
         47.96,
         50.89,
         59.70
       }
     }
    }

};

/**********************************************************/
/*    print to stdout the chi squared distribution table  */
/**********************************************************/
void dump_X2_table(void)
{
    int i,j;
    
    printf("\n\n\t\t\t\tX2 DISTRIBUTION TABLE\n\n");
    printf("\t");
    for (i=MIN_X2_PER;i<MAX_X2_PER+1; i++)
        printf("%s\t",X2.per_str[i]);
    printf("\n\n");
    for (i=MIN_DF; i<MAX_DF+1; i++){
        printf("%s\t",X2.df[i].str);
        for (j=MIN_X2_PER;j<MAX_X2_PER+1; j++)
            printf("%2.3f\t",X2.df[i].level[j]);
        printf("\n");
    }
    printf("\n\n");
}

/**********************************************************************/
/*  given a list of integers, calculate the mean, variance,           */
/*  standard deviation and Z_statistic                                */
/*  Sep 17 1992: changed the variance divisor to be n-1               */
void calc_mean_var_std_dev_Zstat(int *Z_list, int num_Z, double *mean, double *variance, double *std_dev, double *median, double *Z_stat)
{
    int i;
    double tmp = 0.0;

    for (i=0;i<num_Z;i++){
       tmp+=Z_list[i];
    }
    *mean = tmp/(double)num_Z;
    {   int *ptr; 
	alloc_singZ(ptr,num_Z,int,0);
	sort_int_arr(Z_list, num_Z, ptr, INCREASING);
	if ((num_Z % 2) == 0)
	    *median = (double)(Z_list[ptr[num_Z/2]] +Z_list[ptr[num_Z/2-1]]) /2;
	else
	    *median = (double)(Z_list[ptr[num_Z/2]]);

	free_singarr(ptr,int);
    }

    tmp = 0.0;
    for (i=0;i<num_Z;i++)
       tmp+= (Z_list[i] - *mean) *  (Z_list[i] - *mean);

    if (tmp != 0.0){
      *variance = tmp / (double)(num_Z - 1);
       *std_dev = sqrt(*variance);
       *Z_stat = *mean / (sqrt(*variance) / sqrt((double)num_Z));
    }
    else
       *variance = *std_dev = *Z_stat = 0.0;
}

/**********************************************************************/
/*  given a list of double, calculate the mean, variance,             */
/*  standard deviation and Z_statistic                                */
void calc_mean_var_std_dev_Zstat_double(double *Z_list, int num_Z, double *mean, double *variance, double *std_dev, double *median, double *Z_stat)
{
    int i;
    double tmp = 0.0;

    if (num_Z == 0){
	*mean = *median = *variance = *std_dev = *Z_stat = 0.0;
	return;
    }

    for (i=0;i<num_Z;i++){
       tmp+=Z_list[i];
    }
    *mean = tmp/(double)num_Z;
    {   if (num_Z == 1) 
	    *median = Z_list[0];
	else {
	    int *ptr; 
	
	    alloc_singZ(ptr,num_Z,int,0);
	    sort_double_arr(Z_list, num_Z, ptr, INCREASING);
	    if ((num_Z % 2) == 0)
		*median = (double)(Z_list[ptr[num_Z/2]] +
				   Z_list[ptr[num_Z/2 - 1]]) / 2;
	    else
		*median = Z_list[ptr[num_Z/2]];
	    free_singarr(ptr,int);
	}
    }

    tmp = 0.0;
    for (i=0;i<num_Z;i++)
       tmp+= (Z_list[i] - *mean) *  (Z_list[i] - *mean);

    if (tmp != 0.0){
       *variance = tmp / (double)(num_Z - 1);
       *std_dev = sqrt(*variance);
       *Z_stat = *mean / (sqrt(*variance) / sqrt((double)num_Z));
    }
    else
       *variance = *std_dev = *Z_stat = 0.0;
}

/**********************************************************************/
/*     a general analysis routine to test whether or not the the      */
/*     Z_statistic indicates a significant difference                 */
int print_Z_analysis(double Z_stat)
{
    int i;

    printf("%s                                  Reject if\n","");
    printf("%s                         Z > (+%1.3f) or Z < (-%1.3f)\n","",
                                                     Z2tail[GEN_Z_PER].z,
                                                     Z2tail[GEN_Z_PER].z);
    printf("\n");
    printf("%s                                  Z = %1.3f\n","",Z_stat);

    printf("\n\n%s\t\tSUMMARY:\n\t\t-------\n\n","");
    if (fabs(Z_stat) > Z2tail[GEN_Z_PER].z){
     printf("\tThere is a significant difference between the test results\n");
        printf("\tusing a %s confidence level.\n",
                                   Z2tail[GEN_Z_PER].str);
        printf("\n");
        printf("\tFurther, the probablity of there being a difference is\n");
        for (i=GEN_Z_PER;i>0;i--)
            if (fabs(Z_stat) < Z2tail[i-1].z)
                break;
        if (i==MAX_Z_PER)
            printf("\tgreater that %s.\n",Z2tail[0].str);
        else
            printf("\tbetween %s to %s.",Z2tail[i].str,Z2tail[i-1].str);
        return(TEST_DIFF);
    }
    else{
        printf("\tThere is no significant difference between the test\n");
        printf("\tresults using a %s confidence level.\n",
                                            Z2tail[GEN_Z_PER].str);
        printf("\n");
        printf("\tFurther, the probablity of there being a difference is\n");
        for (i=GEN_Z_PER;i<MIN_Z_PER;i++)
            if (fabs(Z_stat) > Z2tail[i+1].z)
                break;
        if (i==MIN_Z_PER)
            printf("\tless than %s.\n",Z2tail[i].str);
        else
            printf("\tbetween %s to %s.",Z2tail[i+1].str,Z2tail[i].str);
        return(NO_DIFF);
    }
}

int Z_pass(double Z_stat)
{
    if (fabs(Z_stat) > Z2tail[GEN_Z_PER].z)
        return(TEST_DIFF);
    else
        return(NO_DIFF);
}

/**********************************************************************/
/* calc_two_sample_z_test_double by Brett 5/10/93                      */ 
/*     a general routine to perform a two-sample z test on two lists  */
/*     of doubles                                                      */
/**********************************************************************/
void calc_two_sample_z_test_double(double *l1, double *l2, int num_l1, int num_l2, double *Z)
{
  double mean_l1,var_l1,sd_l1,med_l1,Z_stat_l1;
  double mean_l2,var_l2,sd_l2,med_l2,Z_stat_l2;
  calc_mean_var_std_dev_Zstat_double(l1,num_l1,&mean_l1,&var_l1,&sd_l1,
				    &med_l1,&Z_stat_l1);
  calc_mean_var_std_dev_Zstat_double(l2,num_l2,&mean_l2,&var_l2,&sd_l2,
				    &med_l2,&Z_stat_l2);
  *Z=((mean_l1-mean_l2)/sqrt(((sd_l1*sd_l1)/(double)num_l1)+((sd_l2*sd_l2)/(double)num_l2)));
  /* print_Z_analysis(*Z);*/
}

/******************************************************************/
/*   Compute the accumulated binomial distribution for the given  */
/*   parameters.  The Formula was taken from                      */
/*   Statistics: Probability, Inference and Decision              */
/*   By Robert L. Winkler and William L. Hays,  Page 206          */
/*                                                                */
/*                        / n \      r      n-r                   */
/*     P( R = r | n, p) = |   |  *  p   *  q                      */
/*                        \ r /                                   */
/*                                                                */
/*        n = number of trials                                    */
/*        R = number of successes in 'n' trials                   */
/*        p = Probability of success                              */
/*        q = 1 - p                                               */
/*                                                                */
/*   Unfortunately, when n is large (> 50?) this algorithm runs   */
/*   out of precision.  Instead, the binomial distribution can be */
/*   approximated by the normal distribution, if np(1-p) > 5.     */
/*   [ "Statistics: Basic Techniques for Solving Applied Problems"*/
/*     Stephen A. Book pps. 91-92]                                */
/******************************************************************/
double compute_acc_binomial(int R, int n, double p)
{
    int r, dbg=0;
    double sum=0.0;

    /* printf("Binomial 7,5,0.5 = %f\n",(double)compute_acc_binomial(1,1,0.1));*/

    if (n * p * (1.0 - p) <= 5.0){
      for (r=0; r <= R; r++) {
        if (dbg) printf("Binomial [nCr(%d,%d) = %f] * ",
                        n,r,n_CHOOSE_r(n,r));
        if (dbg) printf("[pow(%5.3f,%2d) = %f]",
                        p,r,pow(p,(double)r));
        if (dbg) printf(" * [pow(%5.3f,%2d) = %f",
                        p,n-r,pow(1.0-p,(double)(n-r)));
        sum += n_CHOOSE_r(n,r) * pow(p,(double)r) * pow(1.0-p,(double)(n-r));	
        if (dbg) printf("] = %f\n",sum);
      }
    } else {
      /* use the normal approximation */
      double mean, stdev, zstat, flip;
      mean = n*p;
      stdev = sqrt(n*p*(1.0-p));
      flip = (double)R * 0.5;
      zstat = (flip - mean) / stdev;
      sum = normprob((double)fabs(zstat));
      if (zstat < 0.0)  sum = 1.0 - sum;
      if (dbg) printf("Binomial Normal Approx. mean= %f stdev= %f flip= %f"
		      " zstat= %f accdist= %f\n",mean,stdev,flip,zstat,sum);
    }
    return(sum);
}


/*******************************************************************/
/*  Return the result of multipying all the numbers in the sequence*/
/*  from f to to                                                   */
/*******************************************************************/
double seq_mult(int f, int t)
{  
    double sum=1.0;
    int i;
    for (i=f; i<=t; i++)
        sum *= i;
    return(sum);
}

/*******************************************************************/
/* return the result for the N choose R counting Function          */
/*******************************************************************/
double n_CHOOSE_r(int n, int r)
{
    return(seq_mult(r+1,n) / seq_mult(1,n-r));
}

/*******************************************************************/
/*   The following was lifted from Bill Fisher */

/* File normpr1.h                             */
/* Fnorm[i] is the cumulative standard normal */
/*probability of znorm[i]. Probabilities are  */
/* from Winkler & Hays, pp. xiii.             */

/* Expanded by JGF using the dataplot commands*/
/*   let y = sequence 0 0.01 6                */
/*   let x = norcdf (y)                       */
/*   print y x                                */
/* to generate the extra points.              */


/* #ifndef NORMPROB1_HEADER
   #define NORMPROB1_HEADER */

  static int n_normp = 515;
static double znorm[515] =
    {0.0000000E+00, 0.1000000E-01, 0.2000000E-01, 0.3000000E-01, 0.4000000E-01,
     0.5000000E-01, 0.6000000E-01, 0.7000000E-01, 0.8000000E-01, 0.9000000E-01,
     0.9999999E-01, 0.1100000E+00, 0.1200000E+00, 0.1300000E+00, 0.1400000E+00,
     0.1500000E+00, 0.1600000E+00, 0.1700000E+00, 0.1800000E+00, 0.1900000E+00,
     0.2000000E+00, 0.2100000E+00, 0.2200000E+00, 0.2300000E+00, 0.2400000E+00,
     0.2500000E+00, 0.2600000E+00, 0.2700000E+00, 0.2800000E+00, 0.2900000E+00,
     0.3000000E+00, 0.3100000E+00, 0.3200000E+00, 0.3300000E+00, 0.3400000E+00,
     0.3500000E+00, 0.3600000E+00, 0.3700000E+00, 0.3800000E+00, 0.3900000E+00,
     0.4000000E+00, 0.4100000E+00, 0.4200000E+00, 0.4300000E+00, 0.4400000E+00,
     0.4500000E+00, 0.4600000E+00, 0.4700000E+00, 0.4800000E+00, 0.4900000E+00,
     0.5000000E+00, 0.5100000E+00, 0.5200000E+00, 0.5300000E+00, 0.5400000E+00,
     0.5500000E+00, 0.5600000E+00, 0.5700000E+00, 0.5800000E+00, 0.5900000E+00,
     0.6000000E+00, 0.6100000E+00, 0.6200000E+00, 0.6300000E+00, 0.6400000E+00,
     0.6500000E+00, 0.6600000E+00, 0.6700000E+00, 0.6800000E+00, 0.6900000E+00,
     0.7000000E+00, 0.7100000E+00, 0.7200000E+00, 0.7300000E+00, 0.7400000E+00,
     0.7500000E+00, 0.7600000E+00, 0.7700000E+00, 0.7800000E+00, 0.7900000E+00,
     0.8000000E+00, 0.8100000E+00, 0.8200000E+00, 0.8300000E+00, 0.8400000E+00,
     0.8500000E+00, 0.8600000E+00, 0.8700000E+00, 0.8800000E+00, 0.8900000E+00,
     0.9000000E+00, 0.9100000E+00, 0.9200000E+00, 0.9300000E+00, 0.9400000E+00,
     0.9500000E+00, 0.9600000E+00, 0.9700000E+00, 0.9800000E+00, 0.9899999E+00,
     0.1000000E+01, 0.1010000E+01, 0.1020000E+01, 0.1030000E+01, 0.1040000E+01,
     0.1050000E+01, 0.1060000E+01, 0.1070000E+01, 0.1080000E+01, 0.1090000E+01,
     0.1100000E+01, 0.1110000E+01, 0.1120000E+01, 0.1130000E+01, 0.1140000E+01,
     0.1150000E+01, 0.1160000E+01, 0.1170000E+01, 0.1180000E+01, 0.1190000E+01,
     0.1200000E+01, 0.1210000E+01, 0.1220000E+01, 0.1230000E+01, 0.1240000E+01,
     0.1250000E+01, 0.1260000E+01, 0.1270000E+01, 0.1280000E+01, 0.1290000E+01,
     0.1300000E+01, 0.1310000E+01, 0.1320000E+01, 0.1330000E+01, 0.1340000E+01,
     0.1350000E+01, 0.1360000E+01, 0.1370000E+01, 0.1380000E+01, 0.1390000E+01,
     0.1400000E+01, 0.1410000E+01, 0.1420000E+01, 0.1430000E+01, 0.1440000E+01,
     0.1450000E+01, 0.1460000E+01, 0.1470000E+01, 0.1480000E+01, 0.1490000E+01,
     0.1500000E+01, 0.1510000E+01, 0.1520000E+01, 0.1530000E+01, 0.1540000E+01,
     0.1550000E+01, 0.1560000E+01, 0.1570000E+01, 0.1580000E+01, 0.1590000E+01,
     0.1600000E+01, 0.1610000E+01, 0.1620000E+01, 0.1630000E+01, 0.1640000E+01,
     0.1650000E+01, 0.1660000E+01, 0.1670000E+01, 0.1680000E+01, 0.1690000E+01,
     0.1700000E+01, 0.1710000E+01, 0.1720000E+01, 0.1730000E+01, 0.1740000E+01,
     0.1750000E+01, 0.1760000E+01, 0.1770000E+01, 0.1780000E+01, 0.1790000E+01,
     0.1800000E+01, 0.1810000E+01, 0.1820000E+01, 0.1830000E+01, 0.1840000E+01,
     0.1850000E+01, 0.1860000E+01, 0.1870000E+01, 0.1880000E+01, 0.1890000E+01,
     0.1900000E+01, 0.1910000E+01, 0.1920000E+01, 0.1930000E+01, 0.1940000E+01,
     0.1950000E+01, 0.1960000E+01, 0.1970000E+01, 0.1980000E+01, 0.1990000E+01,
     0.2000000E+01, 0.2010000E+01, 0.2020000E+01, 0.2030000E+01, 0.2040000E+01,
     0.2050000E+01, 0.2060000E+01, 0.2070000E+01, 0.2080000E+01, 0.2090000E+01,
     0.2100000E+01, 0.2110000E+01, 0.2120000E+01, 0.2130000E+01, 0.2140000E+01,
     0.2150000E+01, 0.2160000E+01, 0.2170000E+01, 0.2180000E+01, 0.2190000E+01,
     0.2200000E+01, 0.2210000E+01, 0.2220000E+01, 0.2230000E+01, 0.2240000E+01,
     0.2250000E+01, 0.2260000E+01, 0.2270000E+01, 0.2280000E+01, 0.2290000E+01,
     0.2300000E+01, 0.2310000E+01, 0.2320000E+01, 0.2330000E+01, 0.2340000E+01,
     0.2350000E+01, 0.2360000E+01, 0.2370000E+01, 0.2380000E+01, 0.2390000E+01,
     0.2400000E+01, 0.2410000E+01, 0.2420000E+01, 0.2430000E+01, 0.2440000E+01,
     0.2450000E+01, 0.2460000E+01, 0.2470000E+01, 0.2480000E+01, 0.2490000E+01,
     0.2500000E+01, 0.2510000E+01, 0.2520000E+01, 0.2530000E+01, 0.2540000E+01,
     0.2550000E+01, 0.2560000E+01, 0.2570000E+01, 0.2580000E+01, 0.2590000E+01,
     0.2600000E+01, 0.2610000E+01, 0.2620000E+01, 0.2630000E+01, 0.2640000E+01,
     0.2650000E+01, 0.2660000E+01, 0.2670000E+01, 0.2680000E+01, 0.2690000E+01,
     0.2700000E+01, 0.2710000E+01, 0.2720000E+01, 0.2730000E+01, 0.2740000E+01,
     0.2750000E+01, 0.2760000E+01, 0.2770000E+01, 0.2780000E+01, 0.2790000E+01,
     0.2800000E+01, 0.2810000E+01, 0.2820000E+01, 0.2830000E+01, 0.2840000E+01,
     0.2850000E+01, 0.2860000E+01, 0.2870000E+01, 0.2880000E+01, 0.2890000E+01,
     0.2900000E+01, 0.2910000E+01, 0.2920000E+01, 0.2930000E+01, 0.2940000E+01,
     0.2950000E+01, 0.2960000E+01, 0.2970000E+01, 0.2980000E+01, 0.2990000E+01,
     0.3000000E+01, 0.3010000E+01, 0.3020000E+01, 0.3030000E+01, 0.3040000E+01,
     0.3050000E+01, 0.3060000E+01, 0.3070000E+01, 0.3080000E+01, 0.3090000E+01,
     0.3100000E+01, 0.3110000E+01, 0.3120000E+01, 0.3130000E+01, 0.3140000E+01,
     0.3150000E+01, 0.3160000E+01, 0.3170000E+01, 0.3180000E+01, 0.3190000E+01,
     0.3200000E+01, 0.3210000E+01, 0.3220000E+01, 0.3230000E+01, 0.3240000E+01,
     0.3250000E+01, 0.3260000E+01, 0.3270000E+01, 0.3280000E+01, 0.3290000E+01,
     0.3300000E+01, 0.3310000E+01, 0.3320000E+01, 0.3330000E+01, 0.3340000E+01,
     0.3350000E+01, 0.3360000E+01, 0.3370000E+01, 0.3380000E+01, 0.3390000E+01,
     0.3400000E+01, 0.3410000E+01, 0.3420000E+01, 0.3430000E+01, 0.3440000E+01,
     0.3450000E+01, 0.3460000E+01, 0.3470000E+01, 0.3480000E+01, 0.3490000E+01,
     0.3500000E+01, 0.3510000E+01, 0.3520000E+01, 0.3530000E+01, 0.3540000E+01,
     0.3550000E+01, 0.3560000E+01, 0.3570000E+01, 0.3580000E+01, 0.3590000E+01,
     0.3600000E+01, 0.3610000E+01, 0.3620000E+01, 0.3630000E+01, 0.3640000E+01,
     0.3650000E+01, 0.3660000E+01, 0.3670000E+01, 0.3680000E+01, 0.3690000E+01,
     0.3700000E+01, 0.3710000E+01, 0.3720000E+01, 0.3730000E+01, 0.3740000E+01,
     0.3750000E+01, 0.3760000E+01, 0.3770000E+01, 0.3780000E+01, 0.3790000E+01,
     0.3800000E+01, 0.3810000E+01, 0.3820000E+01, 0.3830000E+01, 0.3840000E+01,
     0.3850000E+01, 0.3860000E+01, 0.3870000E+01, 0.3880000E+01, 0.3890000E+01,
     0.3900000E+01, 0.3910000E+01, 0.3920000E+01, 0.3930000E+01, 0.3940000E+01,
     0.3950000E+01, 0.3960000E+01, 0.3970000E+01, 0.3980000E+01, 0.3990000E+01,
     0.4000000E+01, 0.4010000E+01, 0.4020000E+01, 0.4030000E+01, 0.4040000E+01,
     0.4050000E+01, 0.4060000E+01, 0.4070000E+01, 0.4080000E+01, 0.4090000E+01,
     0.4100000E+01, 0.4110000E+01, 0.4120000E+01, 0.4130000E+01, 0.4140000E+01,
     0.4150000E+01, 0.4160000E+01, 0.4170000E+01, 0.4180000E+01, 0.4190000E+01,
     0.4200000E+01, 0.4210000E+01, 0.4220000E+01, 0.4230000E+01, 0.4240000E+01,
     0.4250000E+01, 0.4260000E+01, 0.4270000E+01, 0.4280000E+01, 0.4290000E+01,
     0.4300000E+01, 0.4310000E+01, 0.4320000E+01, 0.4330000E+01, 0.4340000E+01,
     0.4350000E+01, 0.4360000E+01, 0.4370000E+01, 0.4380000E+01, 0.4390000E+01,
     0.4400000E+01, 0.4410000E+01, 0.4420000E+01, 0.4430000E+01, 0.4440000E+01,
     0.4450000E+01, 0.4460000E+01, 0.4470000E+01, 0.4480000E+01, 0.4490000E+01,
     0.4500000E+01, 0.4510000E+01, 0.4520000E+01, 0.4530000E+01, 0.4540000E+01,
     0.4550000E+01, 0.4560000E+01, 0.4570000E+01, 0.4580000E+01, 0.4590000E+01,
     0.4600000E+01, 0.4610000E+01, 0.4620000E+01, 0.4630000E+01, 0.4640000E+01,
     0.4650000E+01, 0.4660000E+01, 0.4670000E+01, 0.4680000E+01, 0.4690000E+01,
     0.4700000E+01, 0.4710000E+01, 0.4720000E+01, 0.4730000E+01, 0.4740000E+01,
     0.4750000E+01, 0.4760000E+01, 0.4770000E+01, 0.4780000E+01, 0.4790000E+01,
     0.4800000E+01, 0.4810000E+01, 0.4820000E+01, 0.4830000E+01, 0.4840000E+01,
     0.4850000E+01, 0.4860000E+01, 0.4870000E+01, 0.4880000E+01, 0.4890000E+01,
     0.4900000E+01, 0.4910000E+01, 0.4920000E+01, 0.4930000E+01, 0.4940000E+01,
     0.4950000E+01, 0.4960000E+01, 0.4970000E+01, 0.4980000E+01, 0.4990000E+01,
     0.5000000E+01, 0.5010000E+01, 0.5020000E+01, 0.5030000E+01, 0.5040000E+01,
     0.5050000E+01, 0.5060000E+01, 0.5070000E+01, 0.5080000E+01, 0.5090000E+01,
     0.5100000E+01, 0.5110000E+01, 0.5120000E+01, 0.5130000E+01, 0.5140000E+01
   };
     

static double Fnorm[515] =
    {0.5000000E+00, 0.5039895E+00, 0.5079783E+00, 0.5119665E+00, 0.5159534E+00,
     0.5199389E+00, 0.5239222E+00, 0.5279033E+00, 0.5318813E+00, 0.5358564E+00,
     0.5398278E+00, 0.5437953E+00, 0.5477583E+00, 0.5517168E+00, 0.5556699E+00,
     0.5596176E+00, 0.5635594E+00, 0.5674949E+00, 0.5714238E+00, 0.5753454E+00,
     0.5792598E+00, 0.5831660E+00, 0.5870644E+00, 0.5909541E+00, 0.5948349E+00,
     0.5987062E+00, 0.6025680E+00, 0.6064197E+00, 0.6102611E+00, 0.6140916E+00,
     0.6179113E+00, 0.6217195E+00, 0.6255157E+00, 0.6292999E+00, 0.6330717E+00,
     0.6368306E+00, 0.6405763E+00, 0.6443087E+00, 0.6480271E+00, 0.6517317E+00,
     0.6554217E+00, 0.6590970E+00, 0.6627572E+00, 0.6664021E+00, 0.6700314E+00,
     0.6736447E+00, 0.6772418E+00, 0.6808225E+00, 0.6843864E+00, 0.6879330E+00,
     0.6914625E+00, 0.6949742E+00, 0.6984682E+00, 0.7019441E+00, 0.7054015E+00,
     0.7088404E+00, 0.7122603E+00, 0.7156612E+00, 0.7190428E+00, 0.7224047E+00,
     0.7257469E+00, 0.7290691E+00, 0.7323712E+00, 0.7356527E+00, 0.7389138E+00,
     0.7421539E+00, 0.7453731E+00, 0.7485712E+00, 0.7517478E+00, 0.7549030E+00,
     0.7580364E+00, 0.7611480E+00, 0.7642375E+00, 0.7673050E+00, 0.7703501E+00,
     0.7733727E+00, 0.7763727E+00, 0.7793501E+00, 0.7823045E+00, 0.7852361E+00,
     0.7881446E+00, 0.7910300E+00, 0.7938920E+00, 0.7967306E+00, 0.7995459E+00,
     0.8023375E+00, 0.8051056E+00, 0.8078499E+00, 0.8105704E+00, 0.8132671E+00,
     0.8159399E+00, 0.8185887E+00, 0.8212136E+00, 0.8238145E+00, 0.8263912E+00,
     0.8289438E+00, 0.8314723E+00, 0.8339767E+00, 0.8364569E+00, 0.8389129E+00,
     0.8413447E+00, 0.8437524E+00, 0.8461357E+00, 0.8484949E+00, 0.8508300E+00,
     0.8531409E+00, 0.8554276E+00, 0.8576903E+00, 0.8599288E+00, 0.8621434E+00,
     0.8643339E+00, 0.8665004E+00, 0.8686431E+00, 0.8707618E+00, 0.8728567E+00,
     0.8749280E+00, 0.8769755E+00, 0.8789995E+00, 0.8809998E+00, 0.8829768E+00,
     0.8849303E+00, 0.8868605E+00, 0.8887675E+00, 0.8906513E+00, 0.8925123E+00,
     0.8943502E+00, 0.8961653E+00, 0.8979576E+00, 0.8997273E+00, 0.9014746E+00,
     0.9031994E+00, 0.9049020E+00, 0.9065824E+00, 0.9082408E+00, 0.9098772E+00,
     0.9114919E+00, 0.9130850E+00, 0.9146565E+00, 0.9162066E+00, 0.9177355E+00,
     0.9192433E+00, 0.9207301E+00, 0.9221961E+00, 0.9236414E+00, 0.9250662E+00,
     0.9264707E+00, 0.9278549E+00, 0.9292191E+00, 0.9305633E+00, 0.9318879E+00,
     0.9331928E+00, 0.9344782E+00, 0.9357445E+00, 0.9369916E+00, 0.9382198E+00,
     0.9394292E+00, 0.9406200E+00, 0.9417924E+00, 0.9429466E+00, 0.9440826E+00,
     0.9452007E+00, 0.9463010E+00, 0.9473839E+00, 0.9484493E+00, 0.9494974E+00,
     0.9505285E+00, 0.9515428E+00, 0.9525403E+00, 0.9535214E+00, 0.9544861E+00,
     0.9554346E+00, 0.9563671E+00, 0.9572838E+00, 0.9581849E+00, 0.9590705E+00,
     0.9599409E+00, 0.9607961E+00, 0.9616365E+00, 0.9624621E+00, 0.9632731E+00,
     0.9640697E+00, 0.9648522E+00, 0.9656205E+00, 0.9663751E+00, 0.9671159E+00,
     0.9678433E+00, 0.9685573E+00, 0.9692581E+00, 0.9699460E+00, 0.9706211E+00,
     0.9712835E+00, 0.9719334E+00, 0.9725711E+00, 0.9731966E+00, 0.9738102E+00,
     0.9744120E+00, 0.9750022E+00, 0.9755809E+00, 0.9761483E+00, 0.9767046E+00,
     0.9772499E+00, 0.9777845E+00, 0.9783084E+00, 0.9788218E+00, 0.9793249E+00,
     0.9798179E+00, 0.9803008E+00, 0.9807739E+00, 0.9812373E+00, 0.9816912E+00,
     0.9821357E+00, 0.9825709E+00, 0.9829971E+00, 0.9834142E+00, 0.9838226E+00,
     0.9842225E+00, 0.9846137E+00, 0.9849966E+00, 0.9853714E+00, 0.9857379E+00,
     0.9860966E+00, 0.9864475E+00, 0.9867907E+00, 0.9871263E+00, 0.9874546E+00,
     0.9877756E+00, 0.9880894E+00, 0.9883962E+00, 0.9886962E+00, 0.9889894E+00,
     0.9892759E+00, 0.9895560E+00, 0.9898296E+00, 0.9900969E+00, 0.9903582E+00,
     0.9906133E+00, 0.9908625E+00, 0.9911060E+00, 0.9913437E+00, 0.9915758E+00,
     0.9918025E+00, 0.9920238E+00, 0.9922398E+00, 0.9924506E+00, 0.9926564E+00,
     0.9928572E+00, 0.9930531E+00, 0.9932443E+00, 0.9934309E+00, 0.9936128E+00,
     0.9937903E+00, 0.9939634E+00, 0.9941322E+00, 0.9942968E+00, 0.9944574E+00,
     0.9946138E+00, 0.9947664E+00, 0.9949151E+00, 0.9950600E+00, 0.9952012E+00,
     0.9953388E+00, 0.9954728E+00, 0.9956035E+00, 0.9957307E+00, 0.9958547E+00,
     0.9959754E+00, 0.9960929E+00, 0.9962074E+00, 0.9963188E+00, 0.9964274E+00,
     0.9965330E+00, 0.9966358E+00, 0.9967359E+00, 0.9968332E+00, 0.9969280E+00,
     0.9970202E+00, 0.9971099E+00, 0.9971972E+00, 0.9972820E+00, 0.9973645E+00,
     0.9974448E+00, 0.9975229E+00, 0.9975988E+00, 0.9976726E+00, 0.9977443E+00,
     0.9978140E+00, 0.9978817E+00, 0.9979476E+00, 0.9980115E+00, 0.9980737E+00,
     0.9981341E+00, 0.9981928E+00, 0.9982498E+00, 0.9983051E+00, 0.9983588E+00,
     0.9984111E+00, 0.9984617E+00, 0.9985110E+00, 0.9985587E+00, 0.9986051E+00,
     0.9986500E+00, 0.9986937E+00, 0.9987361E+00, 0.9987772E+00, 0.9988170E+00,
     0.9988557E+00, 0.9988933E+00, 0.9989296E+00, 0.9989649E+00, 0.9989992E+00,
     0.9990323E+00, 0.9990645E+00, 0.9990957E+00, 0.9991259E+00, 0.9991552E+00,
     0.9991836E+00, 0.9992111E+00, 0.9992377E+00, 0.9992636E+00, 0.9992886E+00,
     0.9993128E+00, 0.9993362E+00, 0.9993590E+00, 0.9993810E+00, 0.9994023E+00,
     0.9994229E+00, 0.9994429E+00, 0.9994622E+00, 0.9994809E+00, 0.9994990E+00,
     0.9995165E+00, 0.9995335E+00, 0.9995499E+00, 0.9995657E+00, 0.9995810E+00,
     0.9995959E+00, 0.9996102E+00, 0.9996241E+00, 0.9996375E+00, 0.9996505E+00,
     0.9996630E+00, 0.9996752E+00, 0.9996868E+00, 0.9996982E+00, 0.9997091E+00,
     0.9997197E+00, 0.9997299E+00, 0.9997397E+00, 0.9997492E+00, 0.9997584E+00,
     0.9997673E+00, 0.9997759E+00, 0.9997842E+00, 0.9997922E+00, 0.9997999E+00,
     0.9998074E+00, 0.9998145E+00, 0.9998215E+00, 0.9998282E+00, 0.9998346E+00,
     0.9998409E+00, 0.9998469E+00, 0.9998527E+00, 0.9998583E+00, 0.9998636E+00,
     0.9998689E+00, 0.9998739E+00, 0.9998787E+00, 0.9998834E+00, 0.9998878E+00,
     0.9998922E+00, 0.9998963E+00, 0.9999003E+00, 0.9999042E+00, 0.9999080E+00,
     0.9999115E+00, 0.9999150E+00, 0.9999183E+00, 0.9999216E+00, 0.9999247E+00,
     0.9999276E+00, 0.9999305E+00, 0.9999332E+00, 0.9999359E+00, 0.9999385E+00,
     0.9999409E+00, 0.9999433E+00, 0.9999456E+00, 0.9999477E+00, 0.9999499E+00,
     0.9999519E+00, 0.9999538E+00, 0.9999557E+00, 0.9999575E+00, 0.9999592E+00,
     0.9999609E+00, 0.9999625E+00, 0.9999641E+00, 0.9999655E+00, 0.9999669E+00,
     0.9999683E+00, 0.9999696E+00, 0.9999709E+00, 0.9999721E+00, 0.9999732E+00,
     0.9999744E+00, 0.9999754E+00, 0.9999765E+00, 0.9999775E+00, 0.9999784E+00,
     0.9999793E+00, 0.9999802E+00, 0.9999810E+00, 0.9999819E+00, 0.9999826E+00,
     0.9999834E+00, 0.9999841E+00, 0.9999847E+00, 0.9999854E+00, 0.9999861E+00,
     0.9999866E+00, 0.9999872E+00, 0.9999878E+00, 0.9999883E+00, 0.9999888E+00,
     0.9999893E+00, 0.9999897E+00, 0.9999902E+00, 0.9999906E+00, 0.9999911E+00,
     0.9999915E+00, 0.9999918E+00, 0.9999922E+00, 0.9999925E+00, 0.9999928E+00,
     0.9999932E+00, 0.9999935E+00, 0.9999938E+00, 0.9999940E+00, 0.9999943E+00,
     0.9999946E+00, 0.9999948E+00, 0.9999951E+00, 0.9999953E+00, 0.9999955E+00,
     0.9999957E+00, 0.9999959E+00, 0.9999961E+00, 0.9999962E+00, 0.9999964E+00,
     0.9999966E+00, 0.9999968E+00, 0.9999969E+00, 0.9999970E+00, 0.9999972E+00,
     0.9999973E+00, 0.9999974E+00, 0.9999976E+00, 0.9999977E+00, 0.9999978E+00,
     0.9999979E+00, 0.9999980E+00, 0.9999981E+00, 0.9999982E+00, 0.9999983E+00,
     0.9999983E+00, 0.9999984E+00, 0.9999985E+00, 0.9999986E+00, 0.9999986E+00,
     0.9999987E+00, 0.9999987E+00, 0.9999988E+00, 0.9999989E+00, 0.9999989E+00,
     0.9999990E+00, 0.9999990E+00, 0.9999991E+00, 0.9999991E+00, 0.9999992E+00,
     0.9999992E+00, 0.9999992E+00, 0.9999993E+00, 0.9999993E+00, 0.9999993E+00,
     0.9999994E+00, 0.9999994E+00, 0.9999995E+00, 0.9999995E+00, 0.9999995E+00,
     0.9999995E+00, 0.9999995E+00, 0.9999996E+00, 0.9999996E+00, 0.9999996E+00,
     0.9999996E+00, 0.9999996E+00, 0.9999996E+00, 0.9999997E+00, 0.9999997E+00,
     0.9999997E+00, 0.9999997E+00, 0.9999998E+00, 0.9999998E+00, 0.9999998E+00,
     0.9999998E+00, 0.9999998E+00, 0.9999998E+00, 0.9999998E+00, 0.9999998E+00,
     0.9999998E+00, 0.9999998E+00, 0.9999998E+00, 0.9999999E+00, 0.9999999E+00
    };

/* #endif */

/* file normpr1.c */

 /************************************************/
 /*  double normprob(z)                          */
 /*  Returns an approximation to the cumulative  */
 /* standard normal probability of z. The value  */
 /* is found by table look-up;  when z lies      */
 /* between two table values, the lower one is   */
 /* returned, as a conservative approximation.   */
 /* If z is higher than any value in the table,  */
 /* the highest tabulated value is returned.     */
 /*                                              */
 /************************************************/
  double normprob(double z)
{int i, ix=0; double zx, p; boolean found;

    if (z < 0.0) zx = -z; else zx = z;
    found = F;
    for (i=1; (i < n_normp)&&(!found); i++)
      {if (zx < znorm[i]) {ix = i - 1; found = T;}
       if (zx == znorm[i]){ix = i;     found = T;}
      }
    if (!found) p = Fnorm[n_normp-1];
    else p = Fnorm[ix];
    return p;
   }