/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */
/* From input file "matmod.p" */


#include <p2c/p2c.h>


/* mathematics modules
   thomas schneider
*/

/* end of the program */

/* begin module version */

#define version         "matmod 2.07 2007mar21 tds/gds"
/*   2007 Mar 21, 2.07: add simpson routine
     1994 Sep  5, 2.06: functional */
/* origin around 1980 */
/* end module version */

/* begin module describe.matmod */
/*
name
      matmod: mathematics modules

synopsis
      matmod(encseq: in, output: out)

files
      encseq: empty or the output of the encode program for testing parameter
         reading routines.
      output: the version of matmod is printed.  successful compilation
         and running of the program indicates that the modules are correct.

description
      self contained modules for mathematical manipulation.
      included is a procedure for linear regression analysis of data pairs,
      a random number generator, newton's method to find roots of functions,
      and routines for reading the parameters from the encode program.

see also
      delmod.p, module.p, encode.p

author
      thomas d. schneider and gary d. stormo

bugs
      none known

technical notes
      the constant n in procedure randomtest determines how many times
      the random number generator will be in a series of tests.  if n
      is small, the the test will be poor, if it is large then the test may
      take a long time.

*/
/* end module describe.matmod */

/* begin module matmod.const */

#define debugging       false   /* for debugging purposes */
/* end module matmod.const */

/* begin module vector.const */

#define maxvecpart      64
    /* the number of elements in a 'part' of a vector */
#define vpagewidth      64


/* the width (in characters) of the output vector
   file from procedure writevector */
/* end module vector.const */


/* begin module package.encode */
/* *********************************************************************** */
/* begin module encode.readencpar */
/* end module encode.readencpar */

/* begin module encode.functions */
/* end module encode.functions */

/* begin module encode.evprint */
/* end module encode.evprint */

/* *********************************************************************** */
/* end module package.encode */

/* begin module encode.info */
/* information about the encode modules

definitions

      vector.  a set of integers (elements) that represent the encoding
of one sequence.  a vector is made up of one or more parameter vectors.
the form of the vector is defined by the parameter list (the entire
linked list pointed to by a paramptr).

      parameter vector (pv).  a set of integers (elements) that represent the
encoding due to one parameter record.  each parameter vector is made up of
one or more window vectors.  the integers (elements) represent the
encoding of one region of the sequence.

      window vector (wv).  a set of integers (elements) that represent the
encoding of a sequence over a particular window.  the length of the
window vector is always four raised to the coding level.  (eg. for
mononucleotides the length is 4, for dinucleotides the length is 16, etc).
a window vector is made up of several elements.

element.  one of the integers in a vector.  this is the number of a
particular n-tide (eg a, gxc, etc) found in the window.

*/
/* end module encode.info */

/* begin module encode.type.param */
/* the following types allow the user to specify parameters which will
   be used to encode the sequences in the book. */
/* spaces are allowed between the encoded bases, and the number of
   bases to be skipped between each encoded pair of bases are kept
   in this linked list of integers */

typedef struct spacelist {
  long skips;   /* bases skipped to next coded base */
  struct spacelist *next;   /* points to next spacing number */
} spacelist;

/* the encoding parameters for each region are stored in these records,
   and the records for each region are connected into a linked list of
   all the encoding parameters for the entire sequences */
typedef enum {
  start, stop
} endpoints;   /* end points of a coding region */

typedef struct parameter {
  /* these are the instructions for doing the coding */
  long range[2];
  /* the bases to be coded by these
                          instructions, relative to the alignedbase */
  long window;   /* the number of bases included in a coding vector */
  long wshift;   /* movement of the window to the next site */
  long coding;
  /* the 'level' at which the coding is being done,
                      i.e., 1: mononucleotides; 2: dinucleotides; ... */
  spacelist *spaces;   /* the spacing between the encoded bases */
  long cshift;   /* shift to the next coding unit in the window */

  /* values calculated at read-in time */
  long wvlength;
  /* length of a window vector.  this is coding
     raised to the 4th power */
  long pvlength;
  /* parameter vector length.  the vector made up
     of a series of window vectors. */
  /* note: pvlength/wvlength is the number of windows in the
     parameter range */

  struct parameter *next;
      /* set of instructions for coding another region */
} parameter;

/* end module encode.type.param */

/* begin module vector.type */
/* vectors are useful for many applications, including the encoding of
   sequences.  this vector type is designed to be flexible enough to be
   used whenever one needs a vector.  it is designed as a linked list so
   it can contain as many elements as are ever needed.  the actual 'vector'
   is a record containing the length and a pointer to the first 'vectorpart'.
   that vectorpart is a record containing an array of the first 'maxvecpart'
   elements  (maxvecpart is a constant, from module vector.const, which must
   be included) and a pointer to the next 'vectorpart'.  the elements are
   of type real so that either integers or reals may be used.  */

typedef struct vectorpart {
  double numbers[maxvecpart];
  struct vectorpart *next;
} vectorpart;

typedef struct vector {
  long length;
  vectorpart *part;
} vector;

/* end module vector.type version = 'matmod 1.55  83 jul 7 tds/gds'; */


/* begin module matmod.var */
Static _TEXT encseq;   /* the output of the encode program */


Static jmp_buf _JL1;


/* end module matmod.var */

/* begin module halt */
Static Void halt()
{
  /* stop the program.  the procedure performs a goto to the end of the
     program.  you must have a label:
        label 1;
     declared, and also the end of the program must have this label:
        1: end.
     examples are in the module libraries.
     this is the only goto in the delila system. */
  printf(" program halt.\n");
  longjmp(_JL1, 1);
}


/* end module halt version = 'delmod 6.44 82 aug 29 tds/gds'; */

/* begin module linear.regression */
/* linear regression procedure
origin: 1980 august 5
      by thomas schneider
      copyright 1986 */

Static Void regress(control, x, y, sumx, sumy, sumxsqd, sumysqd, sumxy, ex,
		    ey, varx, vary, covxy, r, m, b, n)
Char control;
double x, y, *sumx, *sumy, *sumxsqd, *sumysqd, *sumxy, *ex, *ey, *varx, *vary,
       *covxy, *r, *m, *b;
long *n;
{
  /* described below */
  /* the data pairs */
  /* internal records */
  /* mean, variance and covariance of x and y: */
  /* correlation coefficient: */
  /* m = slope, b = y intercept as in y = mx + b */
  /* number of data pairs entered: */
  /* linear regression for variables x and y.
     the control variable has only three acceptable states:
        c clear all variables to zero
          (note that r is set to 2.0, so that a program can check whether
           or not the results have been calculated yet. the recommended test
           is not for equality, since round-off could cause problems.  use
                if r > 1.5 then (no calculations yet)
                           else (calculations done)    )
        e enter the x and y values into the internal sums
        r results are calculated from the sums
  */
  if (control != 'r' && control != 'e' && control != 'c') {
    printf(" linear regression control variable value, \"%c\" is not acceptable.\n",
	   control);
    printf(" it must be in [\"c\",\"e\",\"r\"]\n");
    halt();
    return;
  }

  switch (control) {

  case 'c':  /* clear */
    x = 0.0;
    y = 0.0;
    *sumx = 0.0;
    *sumy = 0.0;
    *sumxsqd = 0.0;
    *sumysqd = 0.0;
    *sumxy = 0.0;
    *ex = 0.0;
    *ey = 0.0;
    *varx = 0.0;
    *vary = 0.0;
    *covxy = 0.0;
    *r = 2.0;
    /* this (impossible) value can act as a flag that
               shows that the calculations have not been made. */

    *m = 0.0;
    *b = 0.0;
    *n = 0;
    break;

  case 'e':  /* enter data */
    *sumx += x;
    *sumy += y;
    *sumxsqd += x * x;
    *sumysqd += y * y;
    *sumxy += x * y;
    (*n)++;
    break;

  case 'r':  /* calculate results */
    /* check for conditions that would bomb the program */
    if (*n == 0) {
      printf("regress: n is 0; no samples recorded\n");
      halt();
    }
    if (*n == 1) {
      printf("regress: n is 1; regression impossible\n");
      halt();
    }

    *ex = *sumx / *n;
    *ey = *sumy / *n;

    *varx = *sumxsqd / *n - *ex * *ex;
    *vary = *sumysqd / *n - *ey * *ey;
    if (*varx == 0.0) {
      printf("regress: variance of x is zero; regression impossible\n");
      halt();
    }
    if (*vary == 0.0) {
      printf("regress: variance of y is zero; regression impossible\n");
      halt();
    }

    *covxy = *sumxy / *n - *ex * *ey;

    *r = *covxy / sqrt(*varx * *vary);

    *m = *covxy / *varx;
    *b = *ey - *m * *ex;
    break;
  }
}


#define pi              3.141592653589
    /* a value used to help mix up the result */


/* end module linear.regression */

/* begin module random */
Static Void random(seed)
double *seed;
{
  /* the seed value is used to generate a random number in the range 0 to 1.
the seed is then altered, to provide the seed of the next random number.
in other words, the seed is the random number.  to obtain the next random
number in the pseudo-random series, simply call random again with the same
seed.  note that the series will not be truely random, since the previous
value determines the next value.  this is "pseudo-random".

an initial seed of 4.0 is suggested.  note that this is not the first
random number.  you can use different seeds to get different pseudo-random
number series.

the expected mean and standard deviation of a flat distribution is
calculated using simple integrals:

         1
mean = / x dx = 0.5
         0
             1           2
variance = / (x - mean)  dx = 1/12;  stdev = sqrt(variance) = .2886751
             0

(the integral for the variance is solved easily by substituting y = x - mean.)

     the algorithm is derived from:
pascal programs for scientists and engineers, alan r. miller
sybex, berkeley ca (1981).
which in turn was an adaptation from hp-35 applications programs */
  double intermediate;   /* an intermediate result */

  intermediate = *seed + pi;
  intermediate = exp(5.0 * log(intermediate));
  *seed = intermediate - (long)intermediate;
}  /* random */

#undef pi


#define n               100   /* the number of calls to random per test */
#define tests           20   /* the number of tests to apply */

#define seed            4.0   /* the seed for generation */


/* end module random */

/* begin module random.test */
Static Void randomtest(table)
_TEXT *table;
{
  /* this procedure tests the random procedure by calling it many
times and printing a table representing the results,  the mean and the
standard deviation for a series of tests.  this random testing procedure
should be run every time random is used on a new computer system. */
  long test;   /* the number of the current test */
  long call;   /* the number of the current call */
  double x = seed;   /* the random variable */
  double sumx;   /* the sum of x during one test */
  double sumxsqd;   /* the sum of x squared */
  double mean;   /* the mean of x in one test */
  double std;   /* the standard deviation of x in one test */

  putc('\n', table->f);
  /* initialize things */
  /* start x.  note that this is not the first random number */
  fprintf(table->f, " test of the random number generator:\n\n");
  fprintf(table->f, "     mean     std dev\n");
  fprintf(table->f, "    (0.5000)  (0.2887)  (expected values)\n");
  fprintf(table->f, " ===================== %ld calls per test\n", (long)n);

  for (test = 1; test <= tests; test++) {
    /* initialize values for one test */
    sumx = 0.0;
    sumxsqd = 0.0;

    for (call = 1; call <= n; call++) {
      random(&x);   /* generate the next random number */
      sumx += x;
      sumxsqd += x * x;
    }

    mean = sumx / n;
    std = sqrt(n / (n - 1.0) * (sumxsqd / n - mean * mean));
    fprintf(table->f, " %10.4f%10.4f\n", mean, std);
  }
}  /* randomtest */

#undef n
#undef tests
#undef seed


/* end module random.test */

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

/* begin module newton.test.functions */
Static double f(x)
double x;
{
  /* a function for testing the newton's method function */
  return log(x);
}  /* f */


Static double fp(x)
double x;
{
  /* the derivative of the function f, for testing the newton's method function */
  return (1 / x);
}  /* fp */


/* end module newton.test.functions */

/* begin module newton */
Static double newton(initialx, tolerance, tries, passes)
double initialx, tolerance;
long tries, *passes;
{
  /* this function uses newton's method to find a root of the function
f(x).  f(x) is coded externally from this function, as is fp(x), the
derivative of f(x) with respect to x.  an initial guess for x, initialx,
must be given along with a tolerance.  the returned value will be within
the region [root-tolerance, root+tolerance].  if the root can not be
found within 'tries' tries, 'passes' will equal 'tries' and the current
value of x will be returned.  the user must determine if newton found
the root by comparing 'passes' to 'tries'. */
  double deltax;   /* the change in x at each pass */
  boolean searching = true;   /* test to continue the search for the root */
  double x = initialx;   /* the independent variable */

  *passes = 0;
  while (searching && *passes < tries) {
    /* writeln(output, passes:5, x:10:5);*/
    deltax = f(x) / fp(x);
    if (fabs(deltax) > tolerance) {
      x -= deltax;
      (*passes)++;
    } else
      searching = false;
  }
  return x;
}  /* newton */


#define initialx        0.00000001   /* makes for a nice demo */

#define maxtries        100   /* do not get into an infinite loop... */

#define tolerance       0.00000001   /* you may want to try zero also... */


/* end module newton */

/* begin module newton.test */
Static Void newtontest(fout)
_TEXT *fout;
{
  /* test the newton method function. */
  double answer;   /* the value returned by newton */
  long passes = 0;   /* number of times newton tried */
  long tries = 0;   /* how many times we ask newton to try */

  fprintf(fout->f, "\n test of newton's method.\n");
  fprintf(fout->f, " the function is the natural logarithm, for which the\n");
  fprintf(fout->f, " root should be 1.00.\n");
  fprintf(fout->f, " the tolerance in the answer is:  %10.8f\n", tolerance);
  fprintf(fout->f, " the initial guess is:%12c%10.8f\n", ' ', initialx);

  while (tries <= passes && tries < maxtries) {
    tries++;
    answer = newton(initialx, tolerance, tries, &passes);
    fprintf(fout->f, " tries: %4ld passes: %4ld answer: %10.8f\n",
	    tries, passes, answer);
  }
}  /* newtontest */

#undef initialx
#undef maxtries
#undef tolerance


/* end module newton.test */

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

/* begin module vector.io */
/* *********************************************************************** */
/* begin module vector.primitives */
/* end module vector.primitives */

/* begin module vector.readvector */
/* end module vector.readvector */

/* begin module vector.writevector */
/* end module vector.writevector */
/* *********************************************************************** */
/* end module vector.io */

/* begin module vector.functions */
/* *********************************************************************** */
/* begin module vector.io */
/* end module vector.io */

/* begin module vector.dotproduct */
/* end module vector.dotproduct */

/* begin module vector.magnitude */
/* end module vector.magnitude */

/* begin module vector.normalize */
/* end module vector.normalize */
/* *********************************************************************** */
/* end module vector.functions */

/* begin module vector.primitives */
/* these functions and procedures do manipulations on vectors. */

Static double vget(v, pos)
vector v;
long pos;
{
  /* this returns from vector 'v' the value of the element at position 'pos' */
  long i;

  if (pos > v.length || pos < 1) {
    printf(
      " error in call to function vget: position %ld is beyond the end of the vector\n",
      pos);
    halt();
  }

  /* move to the correct 'vectorpart' */
  for (i = 1; i <= (pos - 1) / maxvecpart; i++)
    v.part = v.part->next;

  /* get the proper array element from this part */
  return (v.part->numbers[(pos - 1) & (maxvecpart - 1)]);
}


Static Void vput(v, pos, number)
vector *v;
long pos;
double number;
{
  /* this puts into vector 'v' the value of 'number' at position 'pos' */
  long i;
  vectorpart *firstpart;   /* of the vector */

  if (pos > v->length || pos < 1) {
    printf(
      " error in call to function vput: position %ld is beyond the end of the vector\n",
      pos);
    halt();
  }

  /* move to the correct 'vectorpart' */
  firstpart = v->part;
  for (i = 1; i <= (pos - 1) / maxvecpart; i++)
    v->part = v->part->next;

  /* put the 'number' into the proper array element for this part */
  v->part->numbers[(pos - 1) & (maxvecpart - 1)] = number;
  v->part = firstpart;
}


Static Void makevector(v, l)
vector *v;
long l;
{
  /* create the linked list of vector-parts which are needed for a vector
     of length 'l' */
  long numparts;   /* number of parts needed for the vector */
  long i;   /* index to the parts of the vector */
  vectorpart *firstpart;   /* of the vector */
  vectorpart *newpart;   /* of the vector */

  if (l < 1) {
    printf(" makevector: positive length required\n");
    halt();
  }

  v->length = l;
  v->part = (vectorpart *)Malloc(sizeof(vectorpart));
  firstpart = v->part;
  numparts = (v->length - 1) / maxvecpart + 1;
  for (i = 1; i < numparts; i++) {
    newpart = (vectorpart *)Malloc(sizeof(vectorpart));
    v->part->next = newpart;
    v->part = newpart;
  }
  v->part->next = NULL;
  v->part = firstpart;
}


/* end module vector.primitives */

/* begin module vector.readvector */
Static Void readvector(thefile, v)
_TEXT *thefile;
vector *v;
{
  /* read the elements of v into v from 'thefile'.  v must already be set
     up (all the parts created and linked together) before calling.
     'thefile' must contain, from the cursor point until the end of the
      vector, only numbers, either integers or reals; otherwise it will
      bomb.  the vector ends with an end of line so that
      end of file can be tested for. */
  long i, j;   /* indecies */
  long numparts;   /* number of parts in the vector */
  long lastpart;   /* the number of elements in the last vector part */
  vectorpart *firstpart;   /* pointer to the first vector part */

  firstpart = v->part;
  numparts = (v->length - 1) / maxvecpart + 1;
  lastpart = ((v->length - 1) & (maxvecpart - 1)) + 1;
  for (i = 1; i < numparts; i++) {
    for (j = 0; j < maxvecpart; j++)
      fscanf(thefile->f, "%lg", &v->part->numbers[j]);
    v->part = v->part->next;
  }
  for (j = 0; j < lastpart; j++)
    fscanf(thefile->f, "%lg", &v->part->numbers[j]);
  fscanf(thefile->f, "%*[^\n]");
  getc(thefile->f);
  v->part = firstpart;
}


/* end module vector.readvector */

/* begin module vector.writevector */
Static Void writevector(thefile, v, y, z)
_TEXT *thefile;
vector v;
long y, z;
{
  /* writes the elements of 'v' to 'thefile'.  the integers y and z are:
      y: the size of the field for printing the number;
      z: the size of the decimal part of the field,
      if z = 0 then integers will be printed. */
  long pos = 0;   /* posititon in the vector */
  long i, j;   /* indecies */
  long numparts;   /* number of parts in the vector */
  long lastpart;   /* the number of elements in the last vector part */
  vectorpart *firstpart;   /* pointer to the first vector part */
  long x;   /* the number of elements to write on a line */

  x = (long)(vpagewidth / (y + 1.0));
  firstpart = v.part;
  numparts = (v.length - 1) / maxvecpart + 1;
  lastpart = ((v.length - 1) & (maxvecpart - 1)) + 1;
  if (z == 0) {
    for (i = 1; i < numparts; i++) {
      for (j = 0; j < maxvecpart; j++) {
	fprintf(thefile->f, " %*ld",
		(int)y, (long)floor(v.part->numbers[j] + 0.5));
	pos++;
	if (pos % x == 0)
	  putc('\n', thefile->f);
/* p2c: matmod.p, line 607:
 * Note: Using % for possibly-negative arguments [317] */
      }
      v.part = v.part->next;
    }
    for (j = 0; j < lastpart; j++) {
      fprintf(thefile->f, " %*ld",
	      (int)y, (long)floor(v.part->numbers[j] + 0.5));
      pos++;
      if (pos < v.length && pos % x == 0)
	putc('\n', thefile->f);
/* p2c: matmod.p, line 614:
 * Note: Using % for possibly-negative arguments [317] */
    }
  } else {
    for (i = 1; i < numparts; i++) {
      for (j = 0; j < maxvecpart; j++) {
	fprintf(thefile->f, " %*.*f", (int)y, (int)z, v.part->numbers[j]);
	pos++;
	if (pos % x == 0)
	  putc('\n', thefile->f);
/* p2c: matmod.p, line 622:
 * Note: Using % for possibly-negative arguments [317] */
      }
      v.part = v.part->next;
    }
    for (j = 0; j < lastpart; j++) {
      fprintf(thefile->f, " %*.*f", (int)y, (int)z, v.part->numbers[j]);
      pos++;
      if (pos < v.length && pos % x == 0)
	putc('\n', thefile->f);
/* p2c: matmod.p, line 629:
 * Note: Using % for possibly-negative arguments [317] */
    }
  }
  putc('\n', thefile->f);
  v.part = firstpart;
}


/* end module vector.writevector */

/* begin module vector.dotproduct */
Static double dotproduct(vectora, vectorb)
vector vectora, vectorb;
{
  /* this returns the dotproduct of 'vectora' and 'vectorb' */
  long i;   /* an index */
  double j = 0.0;   /* partial products */

  if (vectora.length != vectorb.length) {
    printf(" function dotproduct: vector lengths must be equal\n");
    halt();
  }
  for (i = 1; i <= vectora.length; i++)
    j += vget(vectora, i) * vget(vectorb, i);
  return j;
}


/* end module vector.dotproduct */

/* begin module vector.magnitude */
Static double magnitude(v)
vector *v;
{
  /* find the magnitude (length) of the vector v */
  return sqrt(dotproduct(*v, *v));
}


/* end module vector.magnitude */

/* begin module vector.normalize */
Static Void normalize(v)
vector *v;
{
  /* this replaces the vector v with the congruent vector of unit length,
     i.e.,  the resulting vector v has the property that v.v = 1   */
  long i;   /* index */
  double length;   /* of the unnormalized vector */
  long FORLIM;

  length = magnitude(v);
  FORLIM = v->length;
  for (i = 1; i <= FORLIM; i++)
    vput(v, i, vget(*v, i) / length);
}


/* end module vector.normalize */

/* begin module vector.vset */
Static Void vset(thevalue, v)
double thevalue;
vector *v;
{
  /* set the value thevalue into v */
  long i;   /* an index */
  long FORLIM;

  FORLIM = v->length;
  for (i = 1; i <= FORLIM; i++)
    vput(v, i, thevalue);
}


/* end module vector.vset */

/* begin module vector.vadd */
Static Void vadd(a, b)
vector a, *b;
{
  /* add the vector a into the vector b */
  long i;   /* an index */

  if (a.length != b->length) {
    printf(" function vadd: vector lengths must be equal\n");
    halt();
  }

  for (i = 1; i <= a.length; i++)
    vput(b, i, vget(*b, i) + vget(a, i));
}


/* end module vector.vadd */

/* begin module vector.vsub */
Static Void vsub(a, b)
vector a, *b;
{
  /* subtract the vector a from the vector b */
  long i;   /* an index */

  if (a.length != b->length) {
    printf(" function vsub: vector lengths must be equal\n");
    halt();
  }

  for (i = 1; i <= a.length; i++)
    vput(b, i, vget(*b, i) - vget(a, i));
}


/* end module vector.vsub */

/* begin module vector.vmult */
Static Void vmult(k, v)
double k;
vector *v;
{
  /* multiply the vector v by the scaler k */
  long i;   /* an index */
  long FORLIM;

  FORLIM = v->length;
  for (i = 1; i <= FORLIM; i++)
    vput(v, i, k * vget(*v, i));
}


/* end module vector.vmult */

/* begin module vector.extract */
Static Void extract(a, astart, astop, b, bstart)
vector a;
long astart, astop;
vector *b;
long bstart;
{
  /* extract the segment defined by [astart,astop] of vector a
and insert it into vector b starting at bstart, going in the
positive direction.  any segment of a may be used, so one may
reverse the order of the values.  all checking is done by vget
and vput. */
  long i;   /* index */
  long precalc;   /* a precalculated value to increase speed */

  if (astart <= astop) {
    precalc = bstart - astart;
    for (i = astart; i <= astop; i++)
      vput(b, i + precalc, vget(a, i));
    return;
  }
  precalc = bstart + astart;
  for (i = astart; i >= astop; i--)
    vput(b, precalc - i, vget(a, i));
}


/* end module vector.extract */

/* begin module arccos */
Static double arccos_(x)
double x;
{
  /* Calculate the arccosine.  Pascal only supplies arctan(angle),
so we have to calculate it ourselves.  Begin from

        2           2
     sin angle + cos angle = 1

               2
divide by cos angle to obtain a relationship between tan(angle) and
cos(angle).  Rearrange.
     The special case of x=0 must be dealt with, since this form
would give a division by zero.  The value that must be returned is
pi/2, since cos(pi/2)=0.  A nice way to find this is to note that at a
45 degree angle (pi/4) the tangent is 1, so that pi=4arctan(1),
and pi/2 = 2 * arctan(1).
     Angles larger than pi/2 can be had if x is negative. */
  double pi;   /* 3.1415... */

  if (x < -1.0 || x > 1.0) {
    printf(" function arccos(x): x was not -1<=x<=+1\n");
    halt();
  }

  if (x == 0.0)
    return (2.0 * atan(1.0));
  else if (x > 0.0)
    return atan(sqrt(1.0 / (x * x) - 1.0));
  else {
    pi = 4.0 * atan(1.0);
    return (pi - atan(sqrt(1.0 / (x * x) - 1.0)));
  }
}


/* end module arccos */

/* begin module vector.angle */
Static double angle(a, b)
vector *a, *b;
{
  /* find the angle between the vectors a and b.
to calculate the angle between the vectors we use
     dotproduct(a,b)=magnitude(a)*magnitude(b)*cos(angle)
which requires arccos(angle), supplied by module arccos. */
  return (arccos_(dotproduct(*a, *b) / (magnitude(a) * magnitude(b))));
}


/* end module vector.angle */

/* begin module vector.test */
Static Void vectortest(fout)
_TEXT *fout;
{
  /* test of some of the vector properties */
  vector a, b, c, d, e, f;   /* vectors for trying out the functions */
  long i;   /* index for vector */
  double pi;   /* 3.1415... but as accurate as the machine can calculate */

  makevector(&a, 2L);
  makevector(&b, 2L);
  makevector(&c, 2L);
  makevector(&d, 5L);
  makevector(&e, 3L);
  makevector(&f, 2L);

  vput(&a, 1L, 1.0);
  vput(&a, 2L, 0.0);
  vset(1.0, &b);
  for (i = 1; i <= 5; i++)
    vput(&d, i, (double)i);
  vput(&f, 1L, -1.0);
  vput(&f, 2L, 1.0);

  fprintf(fout->f, "\n test of vector math\n");
  fprintf(fout->f, " vector a: ");
  writevector(fout, a, 1L, 0L);
  fprintf(fout->f, " vector b: ");
  writevector(fout, b, 1L, 0L);
  fprintf(fout->f, " vector d: ");
  writevector(fout, d, 2L, 1L);
  fprintf(fout->f, " vector f: ");
  writevector(fout, f, 2L, 1L);

  fprintf(fout->f, " extract parts of d: \n");
  extract(d, 2L, 3L, &c, 1L);
  fprintf(fout->f, " from 2 to 3: ");
  writevector(fout, c, 1L, 0L);
  extract(d, 5L, 3L, &e, 1L);
  fprintf(fout->f, " from 5 to 3: ");
  writevector(fout, e, 1L, 0L);
  extract(d, 3L, 5L, &e, 1L);
  fprintf(fout->f, " from 3 to 5: ");
  writevector(fout, e, 1L, 0L);
  extract(d, 4L, 2L, &e, 1L);
  fprintf(fout->f, " from 4 to 2: ");
  writevector(fout, e, 1L, 0L);

  fprintf(fout->f, " a + b = ");
  vset(0.0, &c);
  vadd(a, &c);
  vadd(b, &c);
  writevector(fout, c, 1L, 0L);

  fprintf(fout->f, " a - b = ");
  vset(0.0, &c);
  vadd(a, &c);
  vsub(b, &c);
  writevector(fout, c, 1L, 0L);

  fprintf(fout->f, " a is multiplied by 5, a now is ");
  vmult(5.0, &a);
  writevector(fout, a, 1L, 0L);

  fprintf(fout->f, " magnitude(a) = %5.3f\n", magnitude(&a));
  fprintf(fout->f, " magnitude(b) = %5.3f\n\n", magnitude(&b));

  fprintf(fout->f, " test arccos function\n");
  pi = 2 * arccos_(0.0);
  fprintf(fout->f, " pi = %6.4f\n", pi);
  fprintf(fout->f, "cos(arccos(0.00 * pi/4.0)) * 4/pi = %6.4f\n",
	  cos(arccos_(0.00 * pi / 4.0)) * 4 / pi);
  fprintf(fout->f, "cos(arccos(0.12 * pi/4.0)) * 4/pi = %6.4f\n",
	  cos(arccos_(0.12 * pi / 4.0)) * 4 / pi);
  fprintf(fout->f, "cos(arccos(0.25 * pi/4.0)) * 4/pi = %6.4f\n",
	  cos(arccos_(0.25 * pi / 4.0)) * 4 / pi);
  fprintf(fout->f, "cos(arccos(0.37 * pi/4.0)) * 4/pi = %6.4f\n",
	  cos(arccos_(0.37 * pi / 4.0)) * 4 / pi);
  fprintf(fout->f, "cos(arccos(0.50 * pi/4.0)) * 4/pi = %6.4f\n",
	  cos(arccos_(0.50 * pi / 4.0)) * 4 / pi);
  fprintf(fout->f, "cos(arccos(0.75 * pi/4.0)) * 4/pi = %6.4f\n",
	  cos(arccos_(0.75 * pi / 4.0)) * 4 / pi);
  fprintf(fout->f, "cos(arccos(1.00 * pi/4.0)) * 4/pi = %6.4f\n",
	  cos(arccos_(1.00 * pi / 4.0)) * 4 / pi);

  fprintf(fout->f, " angle(a,b) = pi/%6.4f\n", pi / angle(&a, &b));
  fprintf(fout->f, " angle(b,f) = pi/%6.4f\n", pi / angle(&b, &f));
  fprintf(fout->f, " angle(a,f) = pi*%6.4f\n", angle(&a, &f) / pi);
}


/* end module vector.test */

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

/* begin module encode.readencpar */
/* read encode parameters */

Static long vectorsize(param)
parameter *param;
{
  /* determine the size of the vector generated by a string of parameter
     records, begin with 'param' */
  long size = 0;   /* of the vector, so far */

  while (param != NULL) {
    size += param->pvlength;
    param = param->next;
  }
  return size;
}


Static long paramsize(param)
parameter *param;
{
  /* determine the size of the vector generated by one parameter record */
  long rangesize;   /* the number of bases in the range */
  long numwindows;

  /* the number of windows which start in the
              range, taking into account the shifting */
  rangesize = param->range[(long)stop] - param->range[(long)start] + 1;
  numwindows = (long)((rangesize - 1.0) / param->wshift + 1);
  return (numwindows * param->wvlength);
}  /* paramsize */


Static Void readencpar(thefile, param, regions, vectorlength)
_TEXT *thefile;
parameter **param;
long *regions, *vectorlength;
{
  /* read the parameter list from 'thefile' and put the information into the
     linked parameter list, beginning with 'param'; the number of parameter
     records is returned in 'regions'.  the length of each vector is
     vectorlength */
  parameter *aparam;   /* from the parameter list */
  parameter *newparam;   /* the new parameter record from 'thefile' */
  spacelist *firstspaces;   /* in the spaces list */
  spacelist *newspaces;   /* for the spaces list */
  long i, j;   /* indices */
  Char ch;   /* to check for proper positioning */
  long FORLIM, FORLIM1;

  if (*thefile->name != '\0') {
    if (thefile->f != NULL)
      thefile->f = freopen(thefile->name, "r", thefile->f);
    else
      thefile->f = fopen(thefile->name, "r");
  } else
    rewind(thefile->f);
  if (thefile->f == NULL)
    _EscIO2(FileNotFound, thefile->name);
  RESETBUF(thefile->f, Char);
  if (BUFEOF(thefile->f)) {
    printf(" encoded sequence file is empty\n");
    halt();
  }
  aparam = (parameter *)Malloc(sizeof(parameter));
  if (*param == NULL)
    *param = (parameter *)Malloc(sizeof(parameter));
  aparam = *param;
  /* skip the header */
  for (i = 1; i <= 2; i++) {
    fscanf(thefile->f, "%*[^\n]");
    getc(thefile->f);
  }

  /* skip a blank line in the header */
  fscanf(thefile->f, "%*[^\n]");
  getc(thefile->f);

  /* read the number of regions */
  fscanf(thefile->f, "%ld%*[^\n]", regions);
  getc(thefile->f);

  FORLIM = *regions;
  for (i = 1; i <= FORLIM; i++) {
    /* read the range */

    fscanf(thefile->f, "%ld", aparam->range);
    do {
      ch = getc(thefile->f);
      if (ch == '\n')
	ch = ' ';
    } while (ch != 'o');
    fscanf(thefile->f, "%ld%*[^\n]", &aparam->range[(long)stop]);
    getc(thefile->f);
    /* read the window size and shift */
    fscanf(thefile->f, "%ld%*[^\n]", &aparam->window);
    getc(thefile->f);
    fscanf(thefile->f, "%ld%*[^\n]", &aparam->wshift);
    getc(thefile->f);
    /* read the coding level and arrangement */
    fscanf(thefile->f, "%ld", &aparam->coding);
    aparam->spaces = (spacelist *)Malloc(sizeof(spacelist));
    if (aparam->coding > 1) {
      do {
	ch = getc(thefile->f);
	if (ch == '\n')
	  ch = ' ';
      } while (ch != ':');
      firstspaces = (spacelist *)Malloc(sizeof(spacelist));
      aparam->spaces = firstspaces;
      fscanf(thefile->f, "%ld", &aparam->spaces->skips);
      FORLIM1 = aparam->coding;
      for (j = 2; j < FORLIM1; j++) {
	newspaces = (spacelist *)Malloc(sizeof(spacelist));
	aparam->spaces->next = newspaces;
	aparam->spaces = newspaces;
	fscanf(thefile->f, "%ld", &aparam->spaces->skips);
      }
      aparam->spaces->next = NULL;
      aparam->spaces = firstspaces;
    } else
      aparam->spaces = NULL;
    fscanf(thefile->f, "%*[^\n]");
    getc(thefile->f);
    fscanf(thefile->f, "%ld%*[^\n]", &aparam->cshift);
    getc(thefile->f);

    /* set up wvlength */
    aparam->wvlength = (long)floor(exp(aparam->coding * log(4.0)) + 0.5);
    /* set up pvlength */
    aparam->pvlength = paramsize(aparam);

    if (i < *regions) {
      newparam = (parameter *)Malloc(sizeof(parameter));
      aparam->next = newparam;
      aparam = newparam;
    }
  }
  aparam->next = NULL;

  /* read the vector size */
  fscanf(thefile->f, "%ld%*[^\n]", vectorlength);
  getc(thefile->f);
  if (*vectorlength == vectorsize(*param))
    return;
  printf(" vector length in the encoded file\n");
  printf(" does not correspond to the parameters\n");
  halt();
}


/* end module encode.readencpar version = 'matmod 1.55  83 jul 7 tds/gds'; */

/* begin module copy.dummy */
/* these routines are needed to keep enchead happy; they
are available from other module libraries */

Static Void copyaline(a, b)
_TEXT *a, *b;
{   /* dummy */
}


Static long copylines(a, b, i)
_TEXT *a, *b;
long i;
{
  return 0;   /* dummy */
}


/* end module copy.dummy */

/* begin module encode.enchead */
Static Void enchead(encseq, fout, length)
_TEXT *encseq, *fout;
long *length;
{
  /* copy the information from the start of the encoded vector
file encseq to fout.  to make fout have the same form as encseq,
it is necessary for the program which calls this procedure to
identify itself to the fout file on one line.
     the vector lengths is returned.  the encseq file is positioned
for reading of the vectors immediately after enchead.
     in contrast, the fout file is left at the end of the
vector length line so the user may add more information to the line.
after completing the line with writeln, a second writeln is
required before the vectors begin. */
  long regions;   /* number of independently encoded regions */
  long i;   /* the number of lines of parameters */

  if (*encseq->name != '\0') {
    if (encseq->f != NULL)
      encseq->f = freopen(encseq->name, "r", encseq->f);
    else
      encseq->f = fopen(encseq->name, "r");
  } else
    rewind(encseq->f);
  if (encseq->f == NULL)
    _EscIO2(FileNotFound, encseq->name);
  RESETBUF(encseq->f, Char);
  fscanf(encseq->f, "%*[^\n]");
  getc(encseq->f);   /* skip the previous program identification */
  copyaline(encseq, fout);   /* copy book information */
  copyaline(encseq, fout);   /* copy blank line */
  /* copy the encoding parameters to the fout */
  fscanf(encseq->f, "%ld", &regions);
  fprintf(fout->f, " %ld", regions);
  i = regions * 6 + 2;
  /* 6 lines for each parameter and
     2 lines for the coded regions */

  if (copylines(encseq, fout, i) != i) {
    printf(" encseq is missing parameters\n");
    halt();
  }

  /* now get the length */
  fscanf(encseq->f, "%ld%*[^\n]", length);
  getc(encseq->f);
  fprintf(fout->f, " %ld is the vector length", *length);
}


/* end module encode.enchead */

/* begin module encode.functions */
/* these functions are quite handy to help one know where one is
in an encoded vector.  they all take a parameter list as read by
readencpar, and the number of the element of the vector one is interested
in.  this element value runs from 0 to one less than the length of the
parameter vector.  examples are in evprint. */

Static boolean beginpv(aparam, element)
parameter *aparam;
long element;
{
  /* are we at the beginning of a parameter vector? */
  if (aparam != NULL)   /* is there a vector? */
    return (element == 0);
  else
    return true;   /* note circular use of aparam in evprint... */
}  /* beginpv */


Static boolean endpv(aparam, element)
parameter *aparam;
long element;
{
  /* are we at the end of a parameter vector?
the zero element is the first element of the parameter vector. */
  if (aparam != NULL)
    return (element >= aparam->pvlength);
  else
    return true;   /* past the end of all vectors */
}  /* endpv */


Static boolean beginwv(aparam, element)
parameter *aparam;
long element;
{
  /* are we at the beginning of a window vector in the current parameter
vector?  the zero element is the first element of the parameter vector. */
  if (aparam != NULL) {
    return (element % aparam->wvlength == 0);
/* p2c: matmod.p, line 1077:
 * Note: Using % for possibly-negative arguments [317] */
  } else
    return false;   /* no further parameters */
}  /* beginwv */


Static boolean endwv(aparam, element)
parameter *aparam;
long element;
{
  /* are we at the end of a window vector in the current parameter vector?
the zero element is the first element of the parameter vector. */
  if (aparam != NULL) {
    return ((element + 1) % aparam->wvlength == 0);
/* p2c: matmod.p, line 1086:
 * Note: Using % for possibly-negative arguments [317] */
  } else {
    return true;   /* no further parameters */
    /*;writeln(output,'element=',element:1,'   wvlength=',aparam^.wvlength:1); */
  }
}  /* endwv */


Static long encposition(aparam, element)
parameter *aparam;
long element;
{
  /* this returns the aligned position corresponding to the element in
a parmaeter list (aparam).
the zero element is the first element of the parameter vector. */
  return (aparam->range[(long)start] + element / aparam->wvlength * aparam->wshift);
}  /* encposition */


/* end module encode.functions */

/* begin module encode.evprint */
/* procedures for encoded vector printing */

Static Void evln(afile, element, curparam)
_TEXT *afile;
long element;
parameter *curparam;
{
  /* see evprint for definitions.  this procedure is used to provide
writelns to afile when the user wants to control the display better. */
  /* for end of windows: */
  if (beginwv(curparam, element))
    putc('\n', afile->f);
  else if (curparam == NULL)
    putc('\n', afile->f);
  /* end of vector */

  /* for end of parameters: */
  if (beginpv(curparam, element))
    putc('\n', afile->f);
}  /* evln */


#define codingmax       4   /* the maximum coding level allowed */


/* Local variables for evheader: */
struct LOC_evheader {
  _TEXT *afile;
  Char bases[codingmax];   /* the bases */
} ;

Local Void space(i, LINK)
long i;
struct LOC_evheader *LINK;
{
  /* put out i spaces to afile */
  long index;   /* for i */

  for (index = 1; index <= i; index++)
    putc(' ', LINK->afile->f);
}  /* space */

Local Void advance(level, LINK)
long level;
struct LOC_evheader *LINK;
{
  /* advance the bases at the level.  aa -> ac -> ag -> ... -> tt */
  switch (LINK->bases[level-1]) {

  case 'a':
    LINK->bases[level-1] = 'c';
    break;

  case 'c':
    LINK->bases[level-1] = 'g';
    break;

  case 'g':
    LINK->bases[level-1] = 't';
    break;

  case 't':
    LINK->bases[level-1] = 'a';   /* aha.. */
    advance(level - 1, LINK);
    break;
  }
}  /* advance */


Static Void evheader(afile_, posfield, valfield, firstparam, curparam)
_TEXT *afile_;
long posfield, valfield;
parameter *firstparam, **curparam;
{
  /* see evprint for definitions of the variables and how to use evheader.
evheader must be called if you want to provide your own writelns
to the evprint display, ie, when ln is false */
  struct LOC_evheader V;
  long i, j, n;   /* indices for loops */
  long size;
  /* the number of characters in the encoded object.
                    a=1, aa=2, axa=3, etc */
  spacelist *s;   /* for finding size */
  parameter *WITH;
  long FORLIM, FORLIM1, FORLIM2;

  V.afile = afile_;
  if (*curparam == NULL)   /* start or recycle */
    *curparam = firstparam;

  WITH = *curparam;
  if (WITH->coding > codingmax) {
    printf(" evheader: coding levels must not be higher than %ld\n",
	   (long)codingmax);
    halt();
  }

  /* find the size of the encoded object */
  size = WITH->coding;
  s = (*curparam)->spaces;
  while (s != NULL) {
    size += s->skips;
    s = s->next;
  }

  /* see if there is room for printing */
  if (labs(valfield) < size) {
    printf(
      " evheader: valfield (%ld) cannot be less than coding level and spaces (%ld)\n",
      valfield, size);
    halt();
  }

  FORLIM = WITH->coding;
  /* clear bases */
  for (i = 0; i < FORLIM; i++)
    V.bases[i] = 'a';

  /* print the header */
  /* room for space, ' to ' and two numbers */
  space(posfield * 2 + 5, &V);
  FORLIM = WITH->wvlength;
  for (i = 1; i <= FORLIM; i++) {
    space(labs(valfield) - size + 1, &V);
    s = WITH->spaces;
    FORLIM1 = WITH->coding;
    for (n = 0; n < FORLIM1; n++) {
      putc(V.bases[n], V.afile->f);
      if (s != NULL) {
	FORLIM2 = s->skips;
	for (j = 1; j <= FORLIM2; j++)
	  putc('x', V.afile->f);
	s = s->next;
      }
    }

    /* advance the bases except for the last time */
    if (i != WITH->wvlength)
      advance(WITH->coding, &V);
  }
}  /* evheader */

#undef codingmax


/* Local variables for evprint: */
struct LOC_evprint {
  _TEXT *afile;
  long posfield;
  double val;
  long valfield, valdecimal;
  boolean ln;
  parameter **curparam;
  long *element;
} ;

Local Void sider(LINK)
struct LOC_evprint *LINK;
{
  /* produce the side of the display */
  long position;
  parameter *WITH;

  /* the location that a window begins
                        from the aligned base */
  WITH = *LINK->curparam;
  position = encposition(*LINK->curparam, *LINK->element);
  fprintf(LINK->afile->f, " %*ld to %*ld",
	  (int)LINK->posfield, position, (int)LINK->posfield,
	  position + WITH->window - 1);
}  /* sider */

Local Void middler(LINK)
struct LOC_evprint *LINK;
{
  /* produce the middle of the display */
  putc(' ', LINK->afile->f);
  if (LINK->valfield > 0) {   /* to display or not to display... */
    if (LINK->valdecimal <= 0)   /* chose display */
      fprintf(LINK->afile->f, "%*ld",
	      (int)LINK->valfield, (long)floor(LINK->val + 0.5));
    else
      fprintf(LINK->afile->f, "%*.*f",
	      (int)LINK->valfield, (int)LINK->valdecimal, LINK->val);
    return;
  }
  if (LINK->ln) {
    printf(" evprint: if valfield is negative, then ln must be false\n");
    halt();
  }
}  /* middler */

Local Void step(LINK)
struct LOC_evprint *LINK;
{
  /* step to the next position */
  (*LINK->element)++;   /* move along the vector */

  if (endpv(*LINK->curparam, *LINK->element))
      /* are we off the end of the vector? */
      {  /* yes */
    *LINK->curparam = (*LINK->curparam)->next;
    *LINK->element = 0;
  }
}  /* step */


Static Void evprint(afile_, posfield_, val_, valfield_, valdecimal_, ln_,
		    firstparam, curparam_, element_)
_TEXT *afile_;
long posfield_;
double val_;
long valfield_, valdecimal_;
boolean ln_;
parameter *firstparam, **curparam_;
long *element_;
{
  /* the file to print to */
  /* the number of characters in
                          which to print the position */
  /* the value to print */
  /* the number of characters in
                          in which to print the value */
  /* the number of characters
                          devoted to the digits below the
                          decimal point of val.  if <= 0
                          val is printed as an integer */
  /* control writelns (see below) */
  /* the first parameter */
  /* the current encoding parameter */
  /* the position in the vector */
  /* evprint is an encoded vector printing routine.  it is designed to make
generation of display of encoded vectors or calculated vectors based on
encoded vectors easy.  the basic idea is that evprint will print all the
display associated with one vector value.  so at the beginning it prints
an appropriate header from the vector parameters, and for each line it
prints the beginning and end positions of the window.
1) all fields printed out are separated by one space.
2) the value arguments are to be set by the user before calling evprint.
evprint is then called sequentially with all the vector values placed in val.
this allows the data to come from any data structure; it is not forced
to come from a particular array or file type.
3) to initialize, read in firstparam, set the element to 0 and curparam to nil.
evprint will restart at firstparam after each pass across the vector.
4) once the printing of the display has begun,
do not modify the element or curparam yourself.  it is also not wise to
modify any value parameters except val and possibly valdecimal.  the element
argument is used internally by evprint to keep track of the display.
the value arguments, after being set by the user initially, control the
form of the display.

control of writelns:  if ln is true, then evprint will take care of
writelns.  otherwise the user must do the writeln, as:
        if beginpv(curparam,element) then begin
           evheader(afile,posfield,valfield,firstparam,curparam);
           writeln(afile,[your headers for each column]);
        end;
        evprint(afile,...false...)
        if beginwv(curparam,element)
           then write(afile,[your own material on the same line]);
        evln(afile,curparam,element);
note that evprint had moved the element to the next position, so the
test is for begin of wv, not end of wv.

if valfield is negative (but abs(valfield) is a reasonable size) then
val will not be printed, allowing the user to do so after evprint.
this allows more flexible displays.
note: ln must be false to use this feature.
   */
  struct LOC_evprint V;

  V.afile = afile_;
  V.posfield = posfield_;
  V.val = val_;
  V.valfield = valfield_;
  V.valdecimal = valdecimal_;
  V.ln = ln_;
  V.curparam = curparam_;
  V.element = element_;
  if (V.ln) {
    if (beginpv(*V.curparam, *V.element)) {
      evheader(V.afile, V.posfield, V.valfield, firstparam, V.curparam);
      putc('\n', V.afile->f);
    }
  }

  if (beginwv(*V.curparam, *V.element))
    sider(&V);

  middler(&V);

  step(&V);

  if (V.ln)
    evln(V.afile, *V.element, *V.curparam);
}  /* evprint */


/* end module encode.evprint */


/* begin module encode.test */
Static Void enctest(fin, fout)
_TEXT *fin, *fout;
{
  /* test of encode modules:
        encode.type
        encode.vectorsize
        encode.functions
        encode.evprint
the file fin is expected to contain output from the encode program
the file fout is the results of the test */
  parameter *firstparam;   /* pointer to the parameter list */
  parameter *aparam;   /* one of the parameters */
  long regions;   /* number of separate parameter records */
  long vindex;   /* index to v */
  long vectorlength;   /* length of the vector v */
  vector v;   /* a vector of data */
  long element = 0;   /* a value of v */

  putc('\n', fout->f);
  if (*fin->name != '\0') {
    if (fin->f != NULL)
      fin->f = freopen(fin->name, "r", fin->f);
    else
      fin->f = fopen(fin->name, "r");
  } else
    rewind(fin->f);
  if (fin->f == NULL)
    _EscIO2(FileNotFound, fin->name);
  RESETBUF(fin->f, Char);
  if (BUFEOF(fin->f)) {
    fprintf(fout->f, " if file encseq contains the output from the encode\n");
    fprintf(fout->f, " program, then parameter reading will be tested.\n");
    return;
  }
  firstparam = (parameter *)Malloc(sizeof(parameter));
  readencpar(fin, &firstparam, &regions, &vectorlength);
  fprintf(fout->f,
    " the encode file contains vectors with the following %ld sized components:\n",
    regions);
  aparam = firstparam;
  while (aparam != NULL) {
    fprintf(fout->f, " %5ld", paramsize(aparam));
    aparam = aparam->next;
  }
  fprintf(fout->f, "\n the total size of each vector is %ld\n\n",
	  vectorlength);

  fprintf(fout->f, " the first vector is:\n");
  makevector(&v, vectorlength);
  readvector(fin, &v);
  aparam = NULL;
  for (vindex = 1; vindex <= vectorlength; vindex++)
    evprint(fout, 4L, vget(v, vindex), 5L, 0L, true, firstparam, &aparam,
	    &element);
}  /* enctest */


/* end module encode.test */

/* begin module simpson */
Static Void simpson(upper, answer, tol)
double upper, *answer, *tol;
{
  /* Perform a numerical integration of the normal distribution using
Simpson's rule.  The variable "lower" is the lower bound of the
integration.  The variable "upper" is the the upper bound of the
integration.  "answer" is the value of the integration, and "tol"
is the tolerance of the integration procedure, and thus the error
of the calculation.  Written by Mark Shaner, 1992.  Taken from
Miller, Alan R. PASCAL PROGRAMS FOR SCIENTISTS AND ENGINEERS.
SYBEX Inc., 1981. pgs 260-291

2007 Mar 21:  Tom Schneider:  Using a tolerance of 1/maxint caused an
infinite loop on a 64 bit machine.  Presumably it was attempting to
compute the value to 20 or so decimal places.  To avoid this, the
tolerance is set to a more reasonable value.

*/
  /* a counter */
  long i;

  double x, x1, x2;
  /* the independent variables for calculating the value of
              functions */
  double pi = 4.0 * atan(1.0);   /* the value of pi */
  double val;
  /* the value of 1/sqrt(2*pi), it is used in calculating the value of the
          gaussian for any x value, and is defined as a variable in order to
          speed calculations. */
  double deltax;
      /* the distance between every x value and the subsequent one */
  double evensum = 0.0;
      /* the sum of the area under each of the even numbered parabolas */
  double oddsum;
      /* the sum of the area under each of the odd numbered parabolas */
  double endsum;
      /* the sum of the area under the first and last parabolas */
  double endcor;
      /*  the value of the end correction, it is determined using dgauss */
  double sum1;
  /* a place to store the previous sum so that it can be compared
                with the subsequent sum to determine if the tolerance level has
                been reached */
  long pieces = 2;   /* the number of parabolas under the curve */
  double lower = 0.0;   /* the lower bound of the integration */
  double sum;   /* the value of the area under the curve */

  val = 1 / sqrt(2 * pi);
  /* activate both files to be used */
  upper = fabs(upper);
  /*
    tol := 1/maxint;
  */
  *tol = 1.0e-6;
  deltax = (upper - lower) / pieces;
  x = lower + deltax;
  oddsum = val * exp(-0.5 * x * x);
  x1 = lower;
  x2 = upper;
  endsum = val * exp(-0.5 * x1 * x1) + val * exp(-0.5 * x2 * x2);
  endcor = x2 * val * exp(-0.5 * x2 * x2) - x1 * val * exp(-0.5 * x1 * x1);
  sum = (endsum + 4.0 * oddsum) * deltax / 3.0;

  do {
    pieces *= 2;
    sum1 = sum;
    deltax = (upper - lower) / pieces;
    evensum += oddsum;
    oddsum = 0.0;
    for (i = 1; i <= pieces / 2; i++) {
      x = lower + deltax * (2.0 * i - 1.0);
      oddsum += val * exp(-0.5 * x * x);
    }
    sum = (7.0 * endsum + 14.0 * evensum + 16.0 * oddsum + endcor * deltax) *
	  deltax / 15.0;
  } while (fabs(sum - sum1) > fabs(*tol * sum));
  *answer = 0.5 - sum;
  if (*answer < 0.0)   /* safety for roundoff errors */
    *answer = 0.0;
  /*
    writeln (upper:7:5,'      ',answer:7:5,'+/-',tol)
  */
}


/* end module simpson */

main(argc, argv)
int argc;
Char *argv[];
{
  _TEXT TEMP;

  /* begin module matmod.main */
  PASCAL_MAIN(argc, argv);
  if (setjmp(_JL1))
    goto _L1;
  encseq.f = NULL;
  strcpy(encseq.name, "encseq");
  printf(" mathematics modules %s\n", version);

  TEMP.f = stdout;
  *TEMP.name = '\0';
  vectortest(&TEMP);

  TEMP.f = stdout;
  *TEMP.name = '\0';
  newtontest(&TEMP);

  TEMP.f = stdout;
  *TEMP.name = '\0';
  enctest(&encseq, &TEMP);

  TEMP.f = stdout;
  *TEMP.name = '\0';
  randomtest(&TEMP);
  /* end module matmod.main */
_L1:
  if (encseq.f != NULL)
    fclose(encseq.f);
  exit(EXIT_SUCCESS);
}  /* matmod */



/* End. */