diff options
Diffstat (limited to 'generic/tkbltVecMath.C')
| -rw-r--r-- | generic/tkbltVecMath.C | 1609 |
1 files changed, 1609 insertions, 0 deletions
diff --git a/generic/tkbltVecMath.C b/generic/tkbltVecMath.C new file mode 100644 index 0000000..099f5f4 --- /dev/null +++ b/generic/tkbltVecMath.C @@ -0,0 +1,1609 @@ +/* + * 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 "tkbltVecInt.h" +#include "tkbltNsUtil.h" +#include "tkbltParse.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); +} + +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) || std::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 (!std::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 (!std::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 |