From ffedd811557a92c57955cd5f8bbb03a82ddf9b66 Mon Sep 17 00:00:00 2001 From: William Joye Date: Tue, 31 May 2016 13:40:27 -0400 Subject: Rename bltVecMath.C to tkbltVecMath.C --- src/bltVecMath.C | 1603 ---------------------------------------------------- src/tkbltVecMath.C | 1603 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 1603 insertions(+), 1603 deletions(-) delete mode 100644 src/bltVecMath.C create mode 100644 src/tkbltVecMath.C diff --git a/src/bltVecMath.C b/src/bltVecMath.C deleted file mode 100644 index 1e79a7b..0000000 --- a/src/bltVecMath.C +++ /dev/null @@ -1,1603 +0,0 @@ -/* - * 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 -#include -#include -#include -#include - -#include "bltVecInt.h" -#include "bltNsUtil.h" -#include "bltParse.h" - -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); -} - -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)acos}, - {"asin", (void*)ComponentFunc, (ClientData)asin}, - {"atan", (void*)ComponentFunc, (ClientData)atan}, - {"adev", (void*)ScalarFunc, (ClientData)AvgDeviation}, - {"ceil", (void*)ComponentFunc, (ClientData)ceil}, - {"cos", (void*)ComponentFunc, (ClientData)cos}, - {"cosh", (void*)ComponentFunc, (ClientData)cosh}, - {"exp", (void*)ComponentFunc, (ClientData)exp}, - {"floor", (void*)ComponentFunc, (ClientData)floor}, - {"kurtosis",(void*)ScalarFunc, (ClientData)Kurtosis}, - {"length", (void*)ScalarFunc, (ClientData)Length}, - {"log", (void*)ComponentFunc, (ClientData)log}, - {"log10", (void*)ComponentFunc, (ClientData)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)sin}, - {"sinh", (void*)ComponentFunc, (ClientData)sinh}, - {"skew", (void*)ScalarFunc, (ClientData)Skew}, - {"sort", (void*)VectorFunc, (ClientData)Sort}, - {"sqrt", (void*)ComponentFunc, (ClientData)sqrt}, - {"sum", (void*)ScalarFunc, (ClientData)Sum}, - {"tan", (void*)ComponentFunc, (ClientData)tan}, - {"tanh", (void*)ComponentFunc, (ClientData)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 diff --git a/src/tkbltVecMath.C b/src/tkbltVecMath.C new file mode 100644 index 0000000..1e79a7b --- /dev/null +++ b/src/tkbltVecMath.C @@ -0,0 +1,1603 @@ +/* + * 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 +#include +#include +#include +#include + +#include "bltVecInt.h" +#include "bltNsUtil.h" +#include "bltParse.h" + +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); +} + +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)acos}, + {"asin", (void*)ComponentFunc, (ClientData)asin}, + {"atan", (void*)ComponentFunc, (ClientData)atan}, + {"adev", (void*)ScalarFunc, (ClientData)AvgDeviation}, + {"ceil", (void*)ComponentFunc, (ClientData)ceil}, + {"cos", (void*)ComponentFunc, (ClientData)cos}, + {"cosh", (void*)ComponentFunc, (ClientData)cosh}, + {"exp", (void*)ComponentFunc, (ClientData)exp}, + {"floor", (void*)ComponentFunc, (ClientData)floor}, + {"kurtosis",(void*)ScalarFunc, (ClientData)Kurtosis}, + {"length", (void*)ScalarFunc, (ClientData)Length}, + {"log", (void*)ComponentFunc, (ClientData)log}, + {"log10", (void*)ComponentFunc, (ClientData)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)sin}, + {"sinh", (void*)ComponentFunc, (ClientData)sinh}, + {"skew", (void*)ScalarFunc, (ClientData)Skew}, + {"sort", (void*)VectorFunc, (ClientData)Sort}, + {"sqrt", (void*)ComponentFunc, (ClientData)sqrt}, + {"sum", (void*)ScalarFunc, (ClientData)Sum}, + {"tan", (void*)ComponentFunc, (ClientData)tan}, + {"tanh", (void*)ComponentFunc, (ClientData)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 -- cgit v0.12