/*
 * Smithsonian Astrophysical Observatory, Cambridge, MA, USA
 * This code has been modified under the terms listed below and is made
 * available under the same terms.
 */

/*
 *	Copyright 1995-2004 George A Howlett.
 *
 *	Permission is hereby granted, free of charge, to any person
 *	obtaining a copy of this software and associated documentation
 *	files (the "Software"), to deal in the Software without
 *	restriction, including without limitation the rights to use,
 *	copy, modify, merge, publish, distribute, sublicense, and/or
 *	sell copies of the Software, and to permit persons to whom the
 *	Software is furnished to do so, subject to the following
 *	conditions:
 *
 *	The above copyright notice and this permission notice shall be
 *	included in all copies or substantial portions of the
 *	Software.
 *
 *	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY
 *	KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
 *	WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
 *	PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS
 *	OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
 *	OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
 *	OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 *	SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */

#include <cmath>

#include <float.h>
#include <stdlib.h>
#include <errno.h>
#include <ctype.h>
#include <cmath>

#include "tkbltInt.h"
#include "tkbltVecInt.h"
#include "tkbltNsUtil.h"
#include "tkbltParse.h"

using namespace std;
using namespace Blt;

/*
 * Three types of math functions:
 *
 *	ComponentProc		Function is applied in multiple calls to
 *				each component of the vector.
 *	VectorProc		Entire vector is passed, each component is
 *				modified.
 *	ScalarProc		Entire vector is passed, single scalar value
 *				is returned.
 */

typedef double (ComponentProc)(double value);
typedef int (VectorProc)(Vector *vPtr);
typedef double (ScalarProc)(Vector *vPtr);

/*
 * Built-in math functions:
 */
typedef int (GenericMathProc) (void*, Tcl_Interp*, Vector*);

/*
 * MathFunction --
 *
 *	Contains information about math functions that can be called
 *	for vectors.  The table of math functions is global within the
 *	application.  So you can't define two different "sqrt"
 *	functions.
 */
typedef struct {
  const char *name;		/* Name of built-in math function.  If
				 * NULL, indicates that the function
				 * was user-defined and dynamically
				 * allocated.  Function names are
				 * global across all interpreters. */

  void *proc;			/* Procedure that implements this math
				 * function. */

  ClientData clientData;	/* Argument to pass when invoking the
				 * function. */

} MathFunction;

/* The data structure below is used to describe an expression value,
 * which can be either a double-precision floating-point value, or a
 * string.  A given number has only one value at a time.  */

#define STATIC_STRING_SPACE 150

/*
 * Tokens --
 *
 *	The token types are defined below.  In addition, there is a
 *	table associating a precedence with each operator.  The order
 *	of types is important.  Consult the code before changing it.
 */
enum Tokens {
  VALUE, OPEN_PAREN, CLOSE_PAREN, COMMA, END, UNKNOWN,
  MULT = 8, DIVIDE, MOD, PLUS, MINUS,
  LEFT_SHIFT, RIGHT_SHIFT,
  LESS, GREATER, LEQ, GEQ, EQUAL, NEQ,
  OLD_BIT_AND, EXPONENT, OLD_BIT_OR, OLD_QUESTY, OLD_COLON,
  AND, OR, UNARY_MINUS, OLD_UNARY_PLUS, NOT, OLD_BIT_NOT
};

typedef struct {
  Vector *vPtr;
  char staticSpace[STATIC_STRING_SPACE];
  ParseValue pv;		/* Used to hold a string value, if any. */
} Value;

/*
 * ParseInfo --
 *
 *	The data structure below describes the state of parsing an
 *	expression.  It's passed among the routines in this module.
 */
typedef struct {
  const char *expr;		/* The entire right-hand side of the
				 * expression, as originally passed to
				 * Blt_ExprVector. */

  const char *nextPtr;	/* Position of the next character to
			 * be scanned from the expression
			 * string. */

  enum Tokens token;		/* Type of the last token to be parsed
				 * from nextPtr.  See below for
				 * definitions.  Corresponds to the
				 * characters just before nextPtr. */

} ParseInfo;

/*
 * Precedence table.  The values for non-operator token types are ignored.
 */
static int precTable[] =
  {
    0, 0, 0, 0, 0, 0, 0, 0,
    12, 12, 12,			/* MULT, DIVIDE, MOD */
    11, 11,			/* PLUS, MINUS */
    10, 10,			/* LEFT_SHIFT, RIGHT_SHIFT */
    9, 9, 9, 9,			/* LESS, GREATER, LEQ, GEQ */
    8, 8,			/* EQUAL, NEQ */
    7,				/* OLD_BIT_AND */
    13,				/* EXPONENTIATION */
    5,				/* OLD_BIT_OR */
    4,				/* AND */
    3,				/* OR */
    2,				/* OLD_QUESTY */
    1,				/* OLD_COLON */
    14, 14, 14, 14		/* UNARY_MINUS, OLD_UNARY_PLUS, NOT,
				 * OLD_BIT_NOT */
  };


/*
 * Forward declarations.
 */

static int NextValue(Tcl_Interp* interp, ParseInfo *piPtr, int prec, 
		     Value *valuePtr);

static int Sort(Vector *vPtr)
{
  size_t* map = Vec_SortMap(&vPtr, 1);
  double* values = (double*)malloc(sizeof(double) * vPtr->length);
  for(int ii = vPtr->first; ii <= vPtr->last; ii++)
    values[ii] = vPtr->valueArr[map[ii]];

  free(map);
  for (int ii = vPtr->first; ii <= vPtr->last; ii++)
    vPtr->valueArr[ii] = values[ii];

  free(values);
  return TCL_OK;
}

static double Length(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  return (double)(vPtr->last - vPtr->first + 1);
}

double Blt_VecMax(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  return Vec_Max(vPtr);
}

double Blt_VecMin(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  return Vec_Min(vPtr);
}

int Blt_ExprVector(Tcl_Interp* interp, char *string, Blt_Vector *vector)
{
  return ExprVector(interp,string,vector);
}

static double Product(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  double prod;
  double *vp, *vend;

  prod = 1.0;
  for(vp = vPtr->valueArr + vPtr->first,
	vend = vPtr->valueArr + vPtr->last; vp <= vend; vp++) {
    prod *= *vp;
  }
  return prod;
}

static double Sum(Blt_Vector *vectorPtr)
{
  // Kahan summation algorithm

  Vector *vPtr = (Vector *)vectorPtr;
  double* vp = vPtr->valueArr + vPtr->first;
  double sum = *vp++;
  double c = 0.0;			/* A running compensation for lost
				 * low-order bits.*/
  for (double* vend = vPtr->valueArr + vPtr->last; vp <= vend; vp++) {
    double y = *vp - c;		/* So far, so good: c is zero.*/
    double t = sum + y;		/* Alas, sum is big, y small, so
				 * low-order digits of y are lost.*/
    c = (t - sum) - y;	/* (t - sum) recovers the high-order
			 * part of y; subtracting y recovers
			 * -(low part of y) */
    sum = t;
  }

  return sum;
}

static double Mean(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  double sum = Sum(vectorPtr);
  int n = vPtr->last - vPtr->first + 1;

  return sum / (double)n;
}

// var = 1/N Sum( (x[i] - mean)^2 )
static double Variance(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  double mean = Mean(vectorPtr);
  double var = 0.0;
  int count = 0;
  for(double *vp=vPtr->valueArr+vPtr->first, *vend=vPtr->valueArr+vPtr->last;
      vp <= vend; vp++) {
    double dx = *vp - mean;
    var += dx * dx;
    count++;
  }

  if (count < 2)
    return 0.0;

  var /= (double)(count - 1);
  return var;
}

// skew = Sum( (x[i] - mean)^3 ) / (var^3/2)
static double Skew(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  double mean = Mean(vectorPtr);
  double var = 0;
  double skew = 0;
  int count = 0;
  for(double *vp=vPtr->valueArr+vPtr->first, *vend=vPtr->valueArr+vPtr->last;
      vp <= vend; vp++) {
    double diff = *vp - mean;
    diff = fabs(diff);
    double diffsq = diff * diff;
    var += diffsq;
    skew += diffsq * diff;
    count++;
  }

  if (count < 2)
    return 0.0;

  var /= (double)(count - 1);
  skew /= count * var * sqrt(var);
  return skew;
}

static double StdDeviation(Blt_Vector *vectorPtr)
{
  double var;

  var = Variance(vectorPtr);
  if (var > 0.0) {
    return sqrt(var);
  }
  return 0.0;
}

static double AvgDeviation(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  double mean = Mean(vectorPtr);
  double avg = 0.0;
  int count = 0;
  for(double *vp=vPtr->valueArr+vPtr->first, *vend=vPtr->valueArr+vPtr->last;
      vp <= vend; vp++) {
    double diff = *vp - mean;
    avg += fabs(diff);
    count++;
  }

  if (count < 2)
    return 0.0;

  avg /= (double)count;
  return avg;
}

static double Kurtosis(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  double mean = Mean(vectorPtr);
  double var = 0;
  double kurt = 0;
  int count = 0;
  for(double *vp=vPtr->valueArr+vPtr->first, *vend=vPtr->valueArr+vPtr->last;
      vp <= vend; vp++) {
    double diff = *vp - mean;
    double diffsq = diff * diff;
    var += diffsq;
    kurt += diffsq * diffsq;
    count++;
  }

  if (count < 2)
    return 0.0;

  var /= (double)(count - 1);

  if (var == 0.0)
    return 0.0;

  kurt /= (count * var * var);
  return kurt - 3.0;		/* Fisher Kurtosis */
}

static double Median(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  size_t *map;
  double q2;
  int mid;

  if (vPtr->length == 0) {
    return -DBL_MAX;
  }
  map = Vec_SortMap(&vPtr, 1);
  mid = (vPtr->length - 1) / 2;

  /*  
   * Determine Q2 by checking if the number of elements [0..n-1] is
   * odd or even.  If even, we must take the average of the two
   * middle values.  
   */
  if (vPtr->length & 1) { /* Odd */
    q2 = vPtr->valueArr[map[mid]];
  } else {			/* Even */
    q2 = (vPtr->valueArr[map[mid]] + 
	  vPtr->valueArr[map[mid + 1]]) * 0.5;
  }
  free(map);
  return q2;
}

static double Q1(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  double q1;
  size_t *map;

  if (vPtr->length == 0) {
    return -DBL_MAX;
  } 
  map = Vec_SortMap(&vPtr, 1);

  if (vPtr->length < 4) {
    q1 = vPtr->valueArr[map[0]];
  } else {
    int mid, q;

    mid = (vPtr->length - 1) / 2;
    q = mid / 2;

    /* 
     * Determine Q1 by checking if the number of elements in the
     * bottom half [0..mid) is odd or even.   If even, we must
     * take the average of the two middle values.
     */
    if (mid & 1) {		/* Odd */
      q1 = vPtr->valueArr[map[q]]; 
    } else {		/* Even */
      q1 = (vPtr->valueArr[map[q]] + 
	    vPtr->valueArr[map[q + 1]]) * 0.5; 
    }
  }
  free(map);
  return q1;
}

static double Q3(Blt_Vector *vectorPtr)
{
  Vector *vPtr = (Vector *)vectorPtr;
  double q3;
  size_t *map;

  if (vPtr->length == 0) {
    return -DBL_MAX;
  } 

  map = Vec_SortMap(&vPtr, 1);

  if (vPtr->length < 4) {
    q3 = vPtr->valueArr[map[vPtr->length - 1]];
  } else {
    int mid, q;

    mid = (vPtr->length - 1) / 2;
    q = (vPtr->length + mid) / 2;

    /* 
     * Determine Q3 by checking if the number of elements in the
     * upper half (mid..n-1] is odd or even.   If even, we must
     * take the average of the two middle values.
     */
    if (mid & 1) {		/* Odd */
      q3 = vPtr->valueArr[map[q]];
    } else {		/* Even */
      q3 = (vPtr->valueArr[map[q]] + 
	    vPtr->valueArr[map[q + 1]]) * 0.5; 
    }
  }
  free(map);
  return q3;
}

static int Norm(Blt_Vector *vector)
{
  Vector *vPtr = (Vector *)vector;
  double norm, range, min, max;
  int i;

  min = Vec_Min(vPtr);
  max = Vec_Max(vPtr);
  range = max - min;
  for (i = 0; i < vPtr->length; i++) {
    norm = (vPtr->valueArr[i] - min) / range;
    vPtr->valueArr[i] = norm;
  }
  return TCL_OK;
}

static double Nonzeros(Blt_Vector *vector)
{
  Vector *vPtr = (Vector *)vector;
  int count;
  double *vp, *vend;

  count = 0;
  for(vp = vPtr->valueArr + vPtr->first, vend = vPtr->valueArr + vPtr->last; vp <= vend; vp++) {
    if (*vp == 0.0)
      count++;
  }
  return (double) count;
}

static double Fabs(double value)
{
  if (value < 0.0)
    return -value;
  return value;
}

static double Round(double value)
{
  if (value < 0.0)
    return ceil(value - 0.5);
  else
    return floor(value + 0.5);
}

static double Fmod(double x, double y)
{
  if (y == 0.0)
    return 0.0;
  return x - (floor(x / y) * y);
}

/*
 *---------------------------------------------------------------------------
 *
 * MathError --
 *
 *	This procedure is called when an error occurs during a
 *	floating-point operation.  It reads errno and sets
 *	interp->result accordingly.
 *
 * Results:
 *	Interp->result is set to hold an error message.
 *
 * Side effects:
 *	None.
 *
 *---------------------------------------------------------------------------
 */
static void MathError(Tcl_Interp* interp, double value)
{
  if ((errno == EDOM) || (value != value)) {
    Tcl_AppendResult(interp, "domain error: argument not in valid range",
		     (char *)NULL);
    Tcl_SetErrorCode(interp, "ARITH", "DOMAIN", 
		     Tcl_GetStringResult(interp), (char *)NULL);
  }
  else if ((errno == ERANGE) || isinf(value)) {
    if (value == 0.0) {
      Tcl_AppendResult(interp, 
		       "floating-point value too small to represent",
		       (char *)NULL);
      Tcl_SetErrorCode(interp, "ARITH", "UNDERFLOW", 
		       Tcl_GetStringResult(interp), (char *)NULL);
    }
    else {
      Tcl_AppendResult(interp, 
		       "floating-point value too large to represent",
		       (char *)NULL);
      Tcl_SetErrorCode(interp, "ARITH", "OVERFLOW", 
		       Tcl_GetStringResult(interp), (char *)NULL);
    }
  }
  else {
    Tcl_AppendResult(interp, "unknown floating-point error, ",
		     "errno = ", Itoa(errno), (char *)NULL);
    Tcl_SetErrorCode(interp, "ARITH", "UNKNOWN", 
		     Tcl_GetStringResult(interp), (char *)NULL);
  }
}

static int ParseString(Tcl_Interp* interp, const char *string, Value *valuePtr)
{
  const char *endPtr;
  double value;

  errno = 0;

  /*   
   * The string can be either a number or a vector.  First try to
   * convert the string to a number.  If that fails then see if
   * we can find a vector by that name.
   */

  value = strtod(string, (char **)&endPtr);
  if ((endPtr != string) && (*endPtr == '\0')) {
    if (errno != 0) {
      Tcl_ResetResult(interp);
      MathError(interp, value);
      return TCL_ERROR;
    }
    /* Numbers are stored as single element vectors. */
    if (Vec_ChangeLength(interp, valuePtr->vPtr, 1) != TCL_OK) {
      return TCL_ERROR;
    }
    valuePtr->vPtr->valueArr[0] = value;
    return TCL_OK;
  } else {
    Vector *vPtr;

    while (isspace((unsigned char)(*string))) {
      string++;		/* Skip spaces leading the vector name. */    
    }
    vPtr = Vec_ParseElement(interp, valuePtr->vPtr->dataPtr, 
				string, &endPtr, NS_SEARCH_BOTH);
    if (vPtr == NULL) {
      return TCL_ERROR;
    }
    if (*endPtr != '\0') {
      Tcl_AppendResult(interp, "extra characters after vector", 
		       (char *)NULL);
      return TCL_ERROR;
    }
    /* Copy the designated vector to our temporary. */
    Vec_Duplicate(valuePtr->vPtr, vPtr);
  }
  return TCL_OK;
}

static int ParseMathFunction(Tcl_Interp* interp, const char *start,
			     ParseInfo *piPtr, Value *valuePtr)
{
  Tcl_HashEntry *hPtr;
  MathFunction *mathPtr;	/* Info about math function. */
  char *p;
  VectorInterpData *dataPtr;	/* Interpreter-specific data. */
  GenericMathProc *proc;

  /*
   * Find the end of the math function's name and lookup the
   * record for the function.
   */
  p = (char *)start;
  while (isspace((unsigned char)(*p))) {
    p++;
  }
  piPtr->nextPtr = p;
  while (isalnum((unsigned char)(*p)) || (*p == '_')) {
    p++;
  }
  if (*p != '(') {
    return TCL_RETURN;	/* Must start with open parenthesis */
  }
  dataPtr = valuePtr->vPtr->dataPtr;
  *p = '\0';
  hPtr = Tcl_FindHashEntry(&dataPtr->mathProcTable, piPtr->nextPtr);
  *p = '(';
  if (hPtr == NULL) {
    return TCL_RETURN;	/* Name doesn't match any known function */
  }
  /* Pick up the single value as the argument to the function */
  piPtr->token = OPEN_PAREN;
  piPtr->nextPtr = p + 1;
  valuePtr->pv.next = valuePtr->pv.buffer;
  if (NextValue(interp, piPtr, -1, valuePtr) != TCL_OK) {
    return TCL_ERROR;	/* Parse error */
  }
  if (piPtr->token != CLOSE_PAREN) {
    Tcl_AppendResult(interp, "unmatched parentheses in expression \"",
		     piPtr->expr, "\"", (char *)NULL);
    return TCL_ERROR;	/* Missing right parenthesis */
  }
  mathPtr = (MathFunction*)Tcl_GetHashValue(hPtr);
  proc = (GenericMathProc*)mathPtr->proc;
  if ((*proc) (mathPtr->clientData, interp, valuePtr->vPtr) != TCL_OK) {
    return TCL_ERROR;	/* Function invocation error */
  }
  piPtr->token = VALUE;
  return TCL_OK;
}

static int NextToken(Tcl_Interp* interp, ParseInfo *piPtr, Value *valuePtr)
{
  const char *p;
  const char *endPtr;
  const char *var;
  int result;

  p = piPtr->nextPtr;
  while (isspace((unsigned char)(*p))) {
    p++;
  }
  if (*p == '\0') {
    piPtr->token = END;
    piPtr->nextPtr = p;
    return TCL_OK;
  }
  /*
   * Try to parse the token as a floating-point number. But check
   * that the first character isn't a "-" or "+", which "strtod"
   * will happily accept as an unary operator.  Otherwise, we might
   * accidently treat a binary operator as unary by mistake, which
   * will eventually cause a syntax error.
   */
  if ((*p != '-') && (*p != '+')) {
    double value;

    errno = 0;
    value = strtod(p, (char **)&endPtr);
    if (endPtr != p) {
      if (errno != 0) {
	MathError(interp, value);
	return TCL_ERROR;
      }
      piPtr->token = VALUE;
      piPtr->nextPtr = endPtr;

      /*
       * Save the single floating-point value as an 1-component vector.
       */
      if (Vec_ChangeLength(interp, valuePtr->vPtr, 1) != TCL_OK) {
	return TCL_ERROR;
      }
      valuePtr->vPtr->valueArr[0] = value;
      return TCL_OK;
    }
  }
  piPtr->nextPtr = p + 1;
  switch (*p) {
  case '$':
    piPtr->token = VALUE;
    var = Tcl_ParseVar(interp, p, &endPtr);
    if (var == NULL) {
      return TCL_ERROR;
    }
    piPtr->nextPtr = endPtr;
    Tcl_ResetResult(interp);
    result = ParseString(interp, var, valuePtr);
    return result;

  case '[':
    piPtr->token = VALUE;
    result = ParseNestedCmd(interp, p + 1, 0, &endPtr, &valuePtr->pv);
    if (result != TCL_OK) {
      return result;
    }
    piPtr->nextPtr = endPtr;
    Tcl_ResetResult(interp);
    result = ParseString(interp, valuePtr->pv.buffer, valuePtr);
    return result;

  case '"':
    piPtr->token = VALUE;
    result = ParseQuotes(interp, p + 1, '"', 0, &endPtr, &valuePtr->pv);
    if (result != TCL_OK) {
      return result;
    }
    piPtr->nextPtr = endPtr;
    Tcl_ResetResult(interp);
    result = ParseString(interp, valuePtr->pv.buffer, valuePtr);
    return result;

  case '{':
    piPtr->token = VALUE;
    result = ParseBraces(interp, p + 1, &endPtr, &valuePtr->pv);
    if (result != TCL_OK) {
      return result;
    }
    piPtr->nextPtr = endPtr;
    Tcl_ResetResult(interp);
    result = ParseString(interp, valuePtr->pv.buffer, valuePtr);
    return result;

  case '(':
    piPtr->token = OPEN_PAREN;
    break;

  case ')':
    piPtr->token = CLOSE_PAREN;
    break;

  case ',':
    piPtr->token = COMMA;
    break;

  case '*':
    piPtr->token = MULT;
    break;

  case '/':
    piPtr->token = DIVIDE;
    break;

  case '%':
    piPtr->token = MOD;
    break;

  case '+':
    piPtr->token = PLUS;
    break;

  case '-':
    piPtr->token = MINUS;
    break;

  case '^':
    piPtr->token = EXPONENT;
    break;

  case '<':
    switch (*(p + 1)) {
    case '<':
      piPtr->nextPtr = p + 2;
      piPtr->token = LEFT_SHIFT;
      break;
    case '=':
      piPtr->nextPtr = p + 2;
      piPtr->token = LEQ;
      break;
    default:
      piPtr->token = LESS;
      break;
    }
    break;

  case '>':
    switch (*(p + 1)) {
    case '>':
      piPtr->nextPtr = p + 2;
      piPtr->token = RIGHT_SHIFT;
      break;
    case '=':
      piPtr->nextPtr = p + 2;
      piPtr->token = GEQ;
      break;
    default:
      piPtr->token = GREATER;
      break;
    }
    break;

  case '=':
    if (*(p + 1) == '=') {
      piPtr->nextPtr = p + 2;
      piPtr->token = EQUAL;
    } else {
      piPtr->token = UNKNOWN;
    }
    break;

  case '&':
    if (*(p + 1) == '&') {
      piPtr->nextPtr = p + 2;
      piPtr->token = AND;
    } else {
      piPtr->token = UNKNOWN;
    }
    break;

  case '|':
    if (*(p + 1) == '|') {
      piPtr->nextPtr = p + 2;
      piPtr->token = OR;
    } else {
      piPtr->token = UNKNOWN;
    }
    break;

  case '!':
    if (*(p + 1) == '=') {
      piPtr->nextPtr = p + 2;
      piPtr->token = NEQ;
    } else {
      piPtr->token = NOT;
    }
    break;

  default:
    piPtr->token = VALUE;
    result = ParseMathFunction(interp, p, piPtr, valuePtr);
    if ((result == TCL_OK) || (result == TCL_ERROR)) {
      return result;
    } else {
      Vector *vPtr;

      while (isspace((unsigned char)(*p))) {
	p++;		/* Skip spaces leading the vector name. */    
      }
      vPtr = Vec_ParseElement(interp, valuePtr->vPtr->dataPtr, 
				  p, &endPtr, NS_SEARCH_BOTH);
      if (vPtr == NULL) {
	return TCL_ERROR;
      }
      Vec_Duplicate(valuePtr->vPtr, vPtr);
      piPtr->nextPtr = endPtr;
    }
  }
  return TCL_OK;
}

static int NextValue(Tcl_Interp* interp, ParseInfo *piPtr,
		     int prec, Value *valuePtr)
{
  Value value2;		/* Second operand for current operator.  */
  int oper;		/* Current operator (either unary or binary). */
  int gotOp;			/* Non-zero means already lexed the operator
				 * (while picking up value for unary operator).
				 * Don't lex again. */
  int result;
  Vector *vPtr, *v2Ptr;
  int i;

  /*
   * There are two phases to this procedure.  First, pick off an initial
   * value.  Then, parse (binary operator, value) pairs until done.
   */

  vPtr = valuePtr->vPtr;
  v2Ptr = Vec_New(vPtr->dataPtr);
  gotOp = 0;
  value2.vPtr = v2Ptr;
  value2.pv.buffer = value2.pv.next = value2.staticSpace;
  value2.pv.end = value2.pv.buffer + STATIC_STRING_SPACE - 1;
  value2.pv.expandProc = ExpandParseValue;
  value2.pv.clientData = NULL;

  result = NextToken(interp, piPtr, valuePtr);
  if (result != TCL_OK) {
    goto done;
  }
  if (piPtr->token == OPEN_PAREN) {

    /* Parenthesized sub-expression. */

    result = NextValue(interp, piPtr, -1, valuePtr);
    if (result != TCL_OK) {
      goto done;
    }
    if (piPtr->token != CLOSE_PAREN) {
      Tcl_AppendResult(interp, "unmatched parentheses in expression \"",
		       piPtr->expr, "\"", (char *)NULL);
      result = TCL_ERROR;
      goto done;
    }
  } else {
    if (piPtr->token == MINUS) {
      piPtr->token = UNARY_MINUS;
    }
    if (piPtr->token >= UNARY_MINUS) {
      oper = piPtr->token;
      result = NextValue(interp, piPtr, precTable[oper], valuePtr);
      if (result != TCL_OK) {
	goto done;
      }
      gotOp = 1;
      /* Process unary operators. */
      switch (oper) {
      case UNARY_MINUS:
	for(i = 0; i < vPtr->length; i++) {
	  vPtr->valueArr[i] = -(vPtr->valueArr[i]);
	}
	break;

      case NOT:
	for(i = 0; i < vPtr->length; i++) {
	  vPtr->valueArr[i] = (double)(!vPtr->valueArr[i]);
	}
	break;
      default:
	Tcl_AppendResult(interp, "unknown operator", (char *)NULL);
	goto error;
      }
    } else if (piPtr->token != VALUE) {
      Tcl_AppendResult(interp, "missing operand", (char *)NULL);
      goto error;
    }
  }
  if (!gotOp) {
    result = NextToken(interp, piPtr, &value2);
    if (result != TCL_OK) {
      goto done;
    }
  }
  /*
   * Got the first operand.  Now fetch (operator, operand) pairs.
   */
  for (;;) {
    oper = piPtr->token;

    value2.pv.next = value2.pv.buffer;
    if ((oper < MULT) || (oper >= UNARY_MINUS)) {
      if ((oper == END) || (oper == CLOSE_PAREN) || 
	  (oper == COMMA)) {
	result = TCL_OK;
	goto done;
      } else {
	Tcl_AppendResult(interp, "bad operator", (char *)NULL);
	goto error;
      }
    }
    if (precTable[oper] <= prec) {
      result = TCL_OK;
      goto done;
    }
    result = NextValue(interp, piPtr, precTable[oper], &value2);
    if (result != TCL_OK) {
      goto done;
    }
    if ((piPtr->token < MULT) && (piPtr->token != VALUE) &&
	(piPtr->token != END) && (piPtr->token != CLOSE_PAREN) &&
	(piPtr->token != COMMA)) {
      Tcl_AppendResult(interp, "unexpected token in expression",
		       (char *)NULL);
      goto error;
    }
    /*
     * At this point we have two vectors and an operator.
     */

    if (v2Ptr->length == 1) {
      double *opnd;
      double scalar;

      /*
       * 2nd operand is a scalar.
       */
      scalar = v2Ptr->valueArr[0];
      opnd = vPtr->valueArr;
      switch (oper) {
      case MULT:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] *= scalar;
	}
	break;

      case DIVIDE:
	if (scalar == 0.0) {
	  Tcl_AppendResult(interp, "divide by zero", (char *)NULL);
	  goto error;
	}
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] /= scalar;
	}
	break;

      case PLUS:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] += scalar;
	}
	break;

      case MINUS:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] -= scalar;
	}
	break;

      case EXPONENT:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = pow(opnd[i], scalar);
	}
	break;

      case MOD:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = Fmod(opnd[i], scalar);
	}
	break;

      case LESS:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] < scalar);
	}
	break;

      case GREATER:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] > scalar);
	}
	break;

      case LEQ:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] <= scalar);
	}
	break;

      case GEQ:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] >= scalar);
	}
	break;

      case EQUAL:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] == scalar);
	}
	break;

      case NEQ:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] != scalar);
	}
	break;

      case AND:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] && scalar);
	}
	break;

      case OR:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] || scalar);
	}
	break;

      case LEFT_SHIFT:
	{
	  int offset;

	  offset = (int)scalar % vPtr->length;
	  if (offset > 0) {
	    double *hold;
	    int j;

	    hold = (double*)malloc(sizeof(double) * offset);
	    for (i = 0; i < offset; i++) {
	      hold[i] = opnd[i];
	    }
	    for (i = offset, j = 0; i < vPtr->length; i++, j++) {
	      opnd[j] = opnd[i];
	    }
	    for (i = 0, j = vPtr->length - offset;
		 j < vPtr->length; i++, j++) {
	      opnd[j] = hold[i];
	    }
	    free(hold);
	  }
	}
	break;

      case RIGHT_SHIFT:
	{
	  int offset;

	  offset = (int)scalar % vPtr->length;
	  if (offset > 0) {
	    double *hold;
	    int j;
			
	    hold = (double*)malloc(sizeof(double) * offset);
	    for (i = vPtr->length - offset, j = 0; 
		 i < vPtr->length; i++, j++) {
	      hold[j] = opnd[i];
	    }
	    for (i = vPtr->length - offset - 1, 
		   j = vPtr->length - 1; i >= 0; i--, j--) {
	      opnd[j] = opnd[i];
	    }
	    for (i = 0; i < offset; i++) {
	      opnd[i] = hold[i];
	    }
	    free(hold);
	  }
	}
	break;

      default:
	Tcl_AppendResult(interp, "unknown operator in expression",
			 (char *)NULL);
	goto error;
      }

    } else if (vPtr->length == 1) {
      double *opnd;
      double scalar;

      /*
       * 1st operand is a scalar.
       */
      scalar = vPtr->valueArr[0];
      Vec_Duplicate(vPtr, v2Ptr);
      opnd = vPtr->valueArr;
      switch (oper) {
      case MULT:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] *= scalar;
	}
	break;

      case PLUS:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] += scalar;
	}
	break;

      case DIVIDE:
	for(i = 0; i < vPtr->length; i++) {
	  if (opnd[i] == 0.0) {
	    Tcl_AppendResult(interp, "divide by zero",
			     (char *)NULL);
	    goto error;
	  }
	  opnd[i] = (scalar / opnd[i]);
	}
	break;

      case MINUS:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = scalar - opnd[i];
	}
	break;

      case EXPONENT:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = pow(scalar, opnd[i]);
	}
	break;

      case MOD:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = Fmod(scalar, opnd[i]);
	}
	break;

      case LESS:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(scalar < opnd[i]);
	}
	break;

      case GREATER:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(scalar > opnd[i]);
	}
	break;

      case LEQ:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(scalar >= opnd[i]);
	}
	break;

      case GEQ:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(scalar <= opnd[i]);
	}
	break;

      case EQUAL:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] == scalar);
	}
	break;

      case NEQ:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] != scalar);
	}
	break;

      case AND:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] && scalar);
	}
	break;

      case OR:
	for(i = 0; i < vPtr->length; i++) {
	  opnd[i] = (double)(opnd[i] || scalar);
	}
	break;

      case LEFT_SHIFT:
      case RIGHT_SHIFT:
	Tcl_AppendResult(interp, "second shift operand must be scalar",
			 (char *)NULL);
	goto error;

      default:
	Tcl_AppendResult(interp, "unknown operator in expression",
			 (char *)NULL);
	goto error;
      }
    } else {
      double *opnd1, *opnd2;
      /*
       * Carry out the function of the specified operator.
       */
      if (vPtr->length != v2Ptr->length) {
	Tcl_AppendResult(interp, "vectors are different lengths",
			 (char *)NULL);
	goto error;
      }
      opnd1 = vPtr->valueArr, opnd2 = v2Ptr->valueArr;
      switch (oper) {
      case MULT:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] *= opnd2[i];
	}
	break;

      case DIVIDE:
	for (i = 0; i < vPtr->length; i++) {
	  if (opnd2[i] == 0.0) {
	    Tcl_AppendResult(interp,
			     "can't divide by 0.0 vector component",
			     (char *)NULL);
	    goto error;
	  }
	  opnd1[i] /= opnd2[i];
	}
	break;

      case PLUS:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] += opnd2[i];
	}
	break;

      case MINUS:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] -= opnd2[i];
	}
	break;

      case MOD:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = Fmod(opnd1[i], opnd2[i]);
	}
	break;

      case EXPONENT:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = pow(opnd1[i], opnd2[i]);
	}
	break;

      case LESS:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = (double)(opnd1[i] < opnd2[i]);
	}
	break;

      case GREATER:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = (double)(opnd1[i] > opnd2[i]);
	}
	break;

      case LEQ:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = (double)(opnd1[i] <= opnd2[i]);
	}
	break;

      case GEQ:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = (double)(opnd1[i] >= opnd2[i]);
	}
	break;

      case EQUAL:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = (double)(opnd1[i] == opnd2[i]);
	}
	break;

      case NEQ:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = (double)(opnd1[i] != opnd2[i]);
	}
	break;

      case AND:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = (double)(opnd1[i] && opnd2[i]);
	}
	break;

      case OR:
	for (i = 0; i < vPtr->length; i++) {
	  opnd1[i] = (double)(opnd1[i] || opnd2[i]);
	}
	break;

      case LEFT_SHIFT:
      case RIGHT_SHIFT:
	Tcl_AppendResult(interp, "second shift operand must be scalar",
			 (char *)NULL);
	goto error;

      default:
	Tcl_AppendResult(interp, "unknown operator in expression",
			 (char *)NULL);
	goto error;
      }
    }
  }
 done:
  if (value2.pv.buffer != value2.staticSpace) {
    free(value2.pv.buffer);
  }
  Vec_Free(v2Ptr);
  return result;

 error:
  if (value2.pv.buffer != value2.staticSpace) {
    free(value2.pv.buffer);
  }
  Vec_Free(v2Ptr);
  return TCL_ERROR;
}

static int EvaluateExpression(Tcl_Interp* interp, char *string,
			      Value *valuePtr)
{
  ParseInfo info;
  int result;
  Vector *vPtr;
  double *vp, *vend;

  info.expr = info.nextPtr = string;
  valuePtr->pv.buffer = valuePtr->pv.next = valuePtr->staticSpace;
  valuePtr->pv.end = valuePtr->pv.buffer + STATIC_STRING_SPACE - 1;
  valuePtr->pv.expandProc = ExpandParseValue;
  valuePtr->pv.clientData = NULL;

  result = NextValue(interp, &info, -1, valuePtr);
  if (result != TCL_OK) {
    return result;
  }
  if (info.token != END) {
    Tcl_AppendResult(interp, ": syntax error in expression \"",
		     string, "\"", (char *)NULL);
    return TCL_ERROR;
  }
  vPtr = valuePtr->vPtr;

  /* Check for NaN's and overflows. */
  for (vp = vPtr->valueArr, vend = vp + vPtr->length; vp < vend; vp++) {
    if (!isfinite(*vp)) {
      /*
       * IEEE floating-point error.
       */
      MathError(interp, *vp);
      return TCL_ERROR;
    }
  }
  return TCL_OK;
}

static int ComponentFunc(ClientData clientData, Tcl_Interp* interp,
			 Vector *vPtr)
{
  ComponentProc *procPtr = (ComponentProc *) clientData;
  double *vp, *vend;

  errno = 0;
  for(vp = vPtr->valueArr + vPtr->first, 
	vend = vPtr->valueArr + vPtr->last; vp <= vend; vp++) {
    *vp = (*procPtr) (*vp);
    if (errno != 0) {
      MathError(interp, *vp);
      return TCL_ERROR;
    }
    if (!isfinite(*vp)) {
      /*
       * IEEE floating-point error.
       */
      MathError(interp, *vp);
      return TCL_ERROR;
    }
  }
  return TCL_OK;
}

static int ScalarFunc(ClientData clientData, Tcl_Interp* interp, Vector *vPtr)
{
  double value;
  ScalarProc *procPtr = (ScalarProc *) clientData;

  errno = 0;
  value = (*procPtr) (vPtr);
  if (errno != 0) {
    MathError(interp, value);
    return TCL_ERROR;
  }
  if (Vec_ChangeLength(interp, vPtr, 1) != TCL_OK) {
    return TCL_ERROR;
  }
  vPtr->valueArr[0] = value;
  return TCL_OK;
}

static int VectorFunc(ClientData clientData, Tcl_Interp* interp, Vector *vPtr)
{
  VectorProc *procPtr = (VectorProc *) clientData;

  return (*procPtr) (vPtr);
}


static MathFunction mathFunctions[] =
  {
    {"abs",     (void*)ComponentFunc, (ClientData)Fabs},
    {"acos",	(void*)ComponentFunc, (ClientData)(double (*)(double))acos},
    {"asin",	(void*)ComponentFunc, (ClientData)(double (*)(double))asin},
    {"atan",	(void*)ComponentFunc, (ClientData)(double (*)(double))atan},
    {"adev",	(void*)ScalarFunc,    (ClientData)AvgDeviation},
    {"ceil",	(void*)ComponentFunc, (ClientData)(double (*)(double))ceil},
    {"cos",	(void*)ComponentFunc, (ClientData)(double (*)(double))cos},
    {"cosh",	(void*)ComponentFunc, (ClientData)(double (*)(double))cosh},
    {"exp",	(void*)ComponentFunc, (ClientData)(double (*)(double))exp},
    {"floor",	(void*)ComponentFunc, (ClientData)(double (*)(double))floor},
    {"kurtosis",(void*)ScalarFunc,    (ClientData)Kurtosis},
    {"length",	(void*)ScalarFunc,    (ClientData)Length},
    {"log",	(void*)ComponentFunc, (ClientData)(double (*)(double))log},
    {"log10",	(void*)ComponentFunc, (ClientData)(double (*)(double))log10},
    {"max",	(void*)ScalarFunc,    (ClientData)Blt_VecMax},
    {"mean",	(void*)ScalarFunc,    (ClientData)Mean},
    {"median",	(void*)ScalarFunc,    (ClientData)Median},
    {"min",	(void*)ScalarFunc,    (ClientData)Blt_VecMin},
    {"norm",	(void*)VectorFunc,    (ClientData)Norm},
    {"nz",	(void*)ScalarFunc,    (ClientData)Nonzeros},
    {"q1",	(void*)ScalarFunc,    (ClientData)Q1},
    {"q3",	(void*)ScalarFunc,    (ClientData)Q3},
    {"prod",	(void*)ScalarFunc,    (ClientData)Product},
    {"random",	(void*)ComponentFunc, (ClientData)drand48},
    {"round",	(void*)ComponentFunc, (ClientData)Round},
    {"sdev",	(void*)ScalarFunc,    (ClientData)StdDeviation},
    {"sin",	(void*)ComponentFunc, (ClientData)(double (*)(double))sin},
    {"sinh",	(void*)ComponentFunc, (ClientData)(double (*)(double))sinh},
    {"skew",	(void*)ScalarFunc,    (ClientData)Skew},
    {"sort",	(void*)VectorFunc,    (ClientData)Sort},
    {"sqrt",	(void*)ComponentFunc, (ClientData)(double (*)(double))sqrt},
    {"sum",	(void*)ScalarFunc,    (ClientData)Sum},
    {"tan",	(void*)ComponentFunc, (ClientData)(double (*)(double))tan},
    {"tanh",	(void*)ComponentFunc, (ClientData)(double (*)(double))tanh},
    {"var",	(void*)ScalarFunc,    (ClientData)Variance},
    {(char *)NULL,},
  };

void Blt::Vec_InstallMathFunctions(Tcl_HashTable *tablePtr)
{
  MathFunction *mathPtr;

  for (mathPtr = mathFunctions; mathPtr->name != NULL; mathPtr++) {
    Tcl_HashEntry *hPtr;
    int isNew;

    hPtr = Tcl_CreateHashEntry(tablePtr, mathPtr->name, &isNew);
    Tcl_SetHashValue(hPtr, (ClientData)mathPtr);
  }
}

void Blt::Vec_UninstallMathFunctions(Tcl_HashTable *tablePtr)
{
  Tcl_HashEntry *hPtr;
  Tcl_HashSearch cursor;

  for (hPtr = Tcl_FirstHashEntry(tablePtr, &cursor); hPtr != NULL; 
       hPtr = Tcl_NextHashEntry(&cursor)) {
    MathFunction *mathPtr = (MathFunction*)Tcl_GetHashValue(hPtr);
    if (mathPtr->name == NULL)
      free(mathPtr);
  }
}

static void InstallIndexProc(Tcl_HashTable *tablePtr, const char *string,
			     Blt_VectorIndexProc *procPtr)
{
  Tcl_HashEntry *hPtr;
  int dummy;

  hPtr = Tcl_CreateHashEntry(tablePtr, string, &dummy);
  if (procPtr == NULL)
    Tcl_DeleteHashEntry(hPtr);
  else
    Tcl_SetHashValue(hPtr, (ClientData)procPtr);
}

void Blt::Vec_InstallSpecialIndices(Tcl_HashTable *tablePtr)
{
  InstallIndexProc(tablePtr, "min",  Blt_VecMin);
  InstallIndexProc(tablePtr, "max",  Blt_VecMax);
  InstallIndexProc(tablePtr, "mean", Mean);
  InstallIndexProc(tablePtr, "sum",  Sum);
  InstallIndexProc(tablePtr, "prod", Product);
}

int Blt::ExprVector(Tcl_Interp* interp, char *string, Blt_Vector *vector)
{
  VectorInterpData *dataPtr;	/* Interpreter-specific data. */
  Vector *vPtr = (Vector *)vector;
  Value value;

  dataPtr = (vector != NULL) ? vPtr->dataPtr : Vec_GetInterpData(interp);
  value.vPtr = Vec_New(dataPtr);
  if (EvaluateExpression(interp, string, &value) != TCL_OK) {
    Vec_Free(value.vPtr);
    return TCL_ERROR;
  }
  if (vPtr != NULL) {
    Vec_Duplicate(vPtr, value.vPtr);
  } else {
    Tcl_Obj *listObjPtr;
    double *vp, *vend;

    /* No result vector.  Put values in interp->result.  */
    listObjPtr = Tcl_NewListObj(0, (Tcl_Obj **) NULL);
    for (vp = value.vPtr->valueArr, vend = vp + value.vPtr->length; 
	 vp < vend; vp++) {
      Tcl_ListObjAppendElement(interp, listObjPtr, Tcl_NewDoubleObj(*vp));
    }
    Tcl_SetObjResult(interp, listObjPtr);
  }
  Vec_Free(value.vPtr);
  return TCL_OK;
}

#ifdef _WIN32
double drand48(void)
{
  return (double)rand() / (double)RAND_MAX;
}

void srand48(long int seed)
{
  srand(seed);
}
#endif