diff --git a/CHANGES b/CHANGES index d311d63..94e17f3 100644 --- a/CHANGES +++ b/CHANGES @@ -26,6 +26,20 @@ The following are the changes from calc version 2.11.0t10 to date: domul() function in zmil.c thanks to a patch by Ernest Bowen . + Added zero dimensional matrices. A zero dimensional matrix is defined as: + + mat A[] or A = mat[] + + Updated the help/mat file to reflect the current status of matrices + including zero dimensional matrices. + + Added indices() builtin function as written by Ernest Bowen + developed from an idea of Klaus Seistrup + . See help/indices for details. + + Fixed a number of insure warnings as reported by Michel van der List + . + The following are the changes from calc version 2.11.0t8.9.1 to 2.11.0t9.4.5: diff --git a/assocfunc.c b/assocfunc.c index b58b10d..caa2da3 100644 --- a/assocfunc.c +++ b/assocfunc.c @@ -47,8 +47,8 @@ associndex(ASSOC *ap, BOOL create, long dim, VALUE *indices) QCKHASH hash; int i; - if (dim <= 0) { - math_error("No dimensions for indexing association"); + if (dim < 0) { + math_error("Negative dimension for indexing association"); /*NOTREACHED*/ } @@ -217,6 +217,27 @@ assocfindex(ASSOC *ap, long index) } +/* + * Returns the list of indices for an association element with specified + * double-bracket index. + */ +LIST * +associndices(ASSOC *ap, long index) +{ + ASSOCELEM *ep; + LIST *lp; + int i; + + ep = elemindex(ap, index); + if (ep == NULL) + return NULL; + lp = listalloc(); + for (i = 0; i < ep->e_dim; i++) + insertlistlast(lp, &ep->e_indices[i]); + return lp; +} + + /* * Compare two associations to see if they are identical. * Returns TRUE if they are different. diff --git a/calcerr.tbl b/calcerr.tbl index 63b2a47..b1fcf77 100644 --- a/calcerr.tbl +++ b/calcerr.tbl @@ -336,3 +336,5 @@ E_STRCPY Bad argument type for strcpy E_STRNCPY Bad argument type for strncpy E_BACKSLASH Bad argument type for unary backslash E_SETMINUS Bad argument type for setminus +E_INDICES1 Bad first argument type for indices +E_INDICES2 Bad second argument for indices diff --git a/codegen.c b/codegen.c index e880081..26cc121 100644 --- a/codegen.c +++ b/codegen.c @@ -1135,6 +1135,21 @@ getonematrix(int symtype) } rescantoken(); + if (gettoken() == T_LEFTPAREN) { + if (isrvalue(getexprlist())) { + scanerror(T_SEMICOLON, "Lvalue expected"); + return; + } + if (gettoken() != T_RIGHTPAREN) { + scanerror(T_SEMICOLON, "Missing right parenthesis"); + return; + } + getonematrix(symtype); + addop(OP_ASSIGN); + return; + } + rescantoken(); + if (gettoken() != T_LEFTBRACKET) { rescantoken(); scanerror(T_SEMICOLON, "Left-bracket expected"); @@ -1150,23 +1165,32 @@ getonematrix(int symtype) * will patch the correct value back into the opcode. */ if (gettoken() == T_RIGHTBRACKET) { - clearopt(); - patchpc = curfunc->f_opcodecount + 1; - addopone(OP_NUMBER, (long) -1); - clearopt(); - addop(OP_ZERO); - addopone(OP_MATCREATE, dim); - addop(OP_ZERO); - addop(OP_INITFILL); - count = 0; - if (gettoken() == T_ASSIGN) + if (gettoken() == T_ASSIGN) { + clearopt(); + patchpc = curfunc->f_opcodecount + 1; + addopone(OP_NUMBER, (long) -1); + clearopt(); + addop(OP_ZERO); + addopone(OP_MATCREATE, dim); + addop(OP_ZERO); + addop(OP_INITFILL); + count = 0; count = getinitlist(); - else + index = addqconstant(itoq(count)); + if (index < 0) + math_error("Cannot allocate constant"); + curfunc->f_opcodes[patchpc] = index; + return; + } + rescantoken(); + addopone(OP_MATCREATE, 0); + if (gettoken() == T_LEFTBRACKET) { + creatematrix(); + } else { rescantoken(); - index = addqconstant(itoq(count)); - if (index < 0) - math_error("Cannot allocate constant"); - curfunc->f_opcodes[patchpc] = index; + addop(OP_ZERO); + } + addop(OP_INITFILL); return; } @@ -1186,41 +1210,45 @@ creatematrix(void) { long dim; - dim = 1; + dim = 0; - while (TRUE) { + for (;;) { + if (gettoken() == T_RIGHTBRACKET) { + addopone(OP_MATCREATE, dim); + if (gettoken() == T_LEFTBRACKET) { + creatematrix(); + } else { + rescantoken(); + addop(OP_ZERO); + } + addop(OP_INITFILL); + return; + } + rescantoken(); + if (++dim > MAXDIM) { + scanerror(T_SEMICOLON, "Only %ld dimensions allowed", MAXDIM); + return; + } (void) getopassignment(); switch (gettoken()) { case T_RIGHTBRACKET: - case T_COMMA: rescantoken(); + case T_COMMA: addop(OP_ONE); addop(OP_SUB); addop(OP_ZERO); break; case T_COLON: (void) getopassignment(); - break; + switch(gettoken()) { + case T_RIGHTBRACKET: + rescantoken(); + case T_COMMA: + continue; + } + /*FALLTHRU*/ default: rescantoken(); - } - switch (gettoken()) { - case T_RIGHTBRACKET: - addopone(OP_MATCREATE, dim); - if (gettoken() == T_LEFTBRACKET) { - creatematrix(); - } else { - rescantoken(); - addop(OP_ZERO); - } - addop(OP_INITFILL); - return; - case T_COMMA: - if (++dim <= MAXDIM) - break; - scanerror(T_SEMICOLON, "Only %ld dimensions allowed", MAXDIM); - return; - default: scanerror(T_SEMICOLON, "Illegal matrix definition"); return; } @@ -2191,8 +2219,14 @@ getmatargs(void) * finds that the element will be referenced for writing, then * it will call writeindexop to change the flag in the opcode. */ - dim = 1; + dim = 0; + if (gettoken() == T_RIGHTBRACKET) { + addoptwo(OP_INDEXADDR, (long) dim, (long) FALSE); + return; + } + rescantoken(); for (;;) { + ++dim; (void) getopassignment(); switch (gettoken()) { case T_RIGHTBRACKET: @@ -2200,7 +2234,6 @@ getmatargs(void) (long) FALSE); return; case T_COMMA: - dim++; break; default: rescantoken(); diff --git a/func.c b/func.c index 2de27eb..1f22ef1 100644 --- a/func.c +++ b/func.c @@ -91,6 +91,8 @@ extern CONST char *error_table[E__COUNT+2]; /* calc coded error messages */ extern void matrandperm(MATRIX *M); extern void listrandperm(LIST *lp); extern int idungetc(FILEID id, int ch); +extern LIST* associndices(ASSOC *ap, long index); +extern LIST* matindices(MATRIX *mp, long index); extern int stoponerror; @@ -3558,7 +3560,7 @@ f_mattrans(VALUE *vp) if (vp->v_type != V_MAT) return error_value(E_MATTRANS1); - if (vp->v_mat->m_dim != 2) + if (vp->v_mat->m_dim > 2) return error_value(E_MATTRANS2); result.v_type = V_MAT; result.v_mat = mattrans(vp->v_mat); @@ -3569,15 +3571,8 @@ f_mattrans(VALUE *vp) static VALUE f_det(VALUE *vp) { - MATRIX *m; - if (vp->v_type != V_MAT) return error_value(E_DET1); - m = vp->v_mat; - if (m->m_dim != 2) - return error_value(E_DET2); - if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) - return error_value(E_DET3); return matdet(vp->v_mat); } @@ -4464,6 +4459,36 @@ f_assoc(int count, VALUE **vals) } +static VALUE +f_indices(VALUE *v1, VALUE *v2) +{ + VALUE result; + LIST *lp; + + if (v2->v_type != V_NUM || zge31b(v2->v_num->num)) + return error_value(E_INDICES2); + + switch (v1->v_type) { + case V_ASSOC: + lp = associndices(v1->v_assoc, qtoi(v2->v_num)); + break; + case V_MAT: + lp = matindices(v1->v_mat, qtoi(v2->v_num)); + break; + default: + return error_value(E_INDICES1); + } + + result.v_type = V_NULL; + result.v_subtype = V_NOSUBTYPE; + if (lp) { + result.v_type = V_LIST; + result.v_list = lp; + } + return result; +} + + static VALUE f_listinsert(int count, VALUE **vals) { @@ -6763,10 +6788,9 @@ f_blk(int count, VALUE **vals) result.v_type = V_BLOCK; result.v_subtype = V_NOSUBTYPE; - vp = *vals; type = 0; - result.v_subtype = V_NOSUBTYPE; if (count > 0) { + vp = *vals; type = vp->v_type; if (type == V_STR || type == V_NBLOCK || type == V_BLOCK) { vals++; @@ -7552,6 +7576,8 @@ static CONST struct builtin builtins[] = { "integral log of a number base 2"}, {"im", 1, 1, 0, OP_IM, 0, 0, "imaginary part of complex number"}, + {"indices", 2, 2, 0, OP_NOP, 0, f_indices, + "indices of a specified assoc or mat value"}, {"inputlevel", 0, 0, 0, OP_NOP, 0, f_inputlevel, "current input depth"}, {"insert", 2, IN, FA, OP_NOP, 0, f_listinsert, diff --git a/help/Makefile b/help/Makefile index 8eb2662..0d25c65 100644 --- a/help/Makefile +++ b/help/Makefile @@ -104,16 +104,16 @@ BLT_HELP_FILES= ${BLT_HELP_FILES_3} ${BLT_HELP_FILES_5} \ # This list is prodiced by the detaillist rule when no WARNINGS are detected. # DETAIL_HELP= abs access acos acosh acot acoth acsc acsch address agd append \ - appr arg arrow asec asech asin asinh assign atan atan2 atanh avg \ - base bit blk blkcpy blkfree blocks bround btrunc calclevel ceil \ - cfappr cfsim char cmdbuf cmp comb conj cos cosh cot coth count cp \ - csc csch ctime delete den dereference det digit digits dp epsilon \ - errcount errmax errno error eval exp fact factor fclose fcnt feof \ - ferror fflush fgetc fgetfield fgetline fgets fgetstr fib files \ - floor fopen forall fprintf fputc fputs fputstr frac free freeglobals \ - freeredc freestatics frem freopen fscan fscanf fseek fsize ftell gcd \ - gcdrem gd getenv hash head highbit hmean hnrmod hypot ilog ilog10 \ - ilog2 im inputlevel insert int inverse iroot isassoc isatty isblk \ + appr arg arrow asec asech asin asinh assign atan atan2 atanh avg base \ + bit blk blkcpy blkfree blocks bround btrunc calclevel ceil cfappr \ + cfsim char cmdbuf cmp comb conj cos cosh cot coth count cp csc csch \ + ctime delete den dereference det digit digits dp epsilon errcount \ + errmax errno error eval exp fact factor fclose fcnt feof ferror \ + fflush fgetc fgetfield fgetline fgets fgetstr fib files floor fopen \ + forall fprintf fputc fputs fputstr frac free freeglobals freeredc \ + freestatics frem freopen fscan fscanf fseek fsize ftell gcd gcdrem \ + gd getenv hash head highbit hmean hnrmod hypot ilog ilog10 ilog2 \ + im indices inputlevel insert int inverse iroot isassoc isatty isblk \ isconfig isdefined iserror iseven isfile ishash isident isint islist \ ismat ismult isnull isnum isobj isobjtype isodd isprime isptr isqrt \ isrand israndom isreal isrel issimple issq isstr istype jacobi join \ diff --git a/help/indices b/help/indices new file mode 100644 index 0000000..ed7acb6 --- /dev/null +++ b/help/indices @@ -0,0 +1,58 @@ +NAME + indices - indices for specified matrix or association element + +SYNOPSIS + indices(V, index) + +TYPES + V matrix or association + index integer + + return list with up to 4 elements + +DESCRIPTION + For 0 <= index < size(V), indices(V, index) returns list(i_0, i_1, ...) + for which V[i_0, i_1, ...] is the same lvalue as V[[index]]. + + For other values of index, a null value is returned. + + This function can be useful for determining those elements for which + the indices satisfy some condition. This is particularly so for + associations since these have no simple relation between the + double-bracket index and the single-bracket indices, which may be + non-integer numbers or strings or other types of value. The + information provided by indices() is often required after the use + of search() or rsearch() which, when successful, return the + double-bracket index of the item found. + +EXAMPLE + > mat M[2,3,1:5] + + > indices(M, 11) + list (3 elements, 2 nonzero): + [[0]] = 0 + [[1]] = 2 + [[2]] = 2 + + > A = assoc(); + + > A["cat", "dog"] = "fight"; + > A[2,3,5,7] = "primes"; + > A["square", 3] = 9 + + > indices(A, search(A, "primes")) + list (4 elements, 4 nonzero): + [[0]] = 2 + [[1]] = 3 + [[2]] = 5 + [[3]] = 7 + +LIMITS + abs(index) < 2^31 + +LIBRARY + LIST* associndices(ASSOC *ap, long index) + LIST* matindices(MATRIX *mp, long index) + +SEE ALSO + assoc, mat diff --git a/help/mat b/help/mat index 414c538..df252ef 100644 --- a/help/mat +++ b/help/mat @@ -1,102 +1,397 @@ -Using matrices +NAME + mat - keyword to create a matrix value - Matrices can have from 1 to 4 dimensions, and are indexed by a - normal-sized integer. The lower and upper bounds of a matrix can - be specified at runtime. The elements of a matrix are defaulted - to zeroes, but can be assigned to be of any type. Thus matrices - can hold complex numbers, strings, objects, etc. Matrices are - stored in memory as an array so that random access to the elements - is easy. +SYNOPSIS + mat [index-range-list] [ = {value_0. ...} ] + mat [] [= {value_0, ...}] + mat variable_1 ... [index-range-list] [ = {value_0, ...} ] + mat variable_1 ... [] [ = {value_0, ...} ] - Matrices are normally indexed using square brackets. If the matrix - is multi-dimensional, then an element can be indexed either by - using multiple pairs of square brackets (as in C), or else by - separating the indexes by commas. Thus the following two statements - reference the same matrix element: + mat [index-range-list_1[index-ranges-list_2] ... [ = { { ...} ...} ] - x = name[3][5]; - x = name[3,5]; + decl id_1 id_2 ... [index-range-list] ... - The double-square bracket operator can be used on any matrix to - make references to the elements easy and efficient. This operator - bypasses the normal indexing mechanism, and treats the array as if - it was one-dimensional and with a lower bound of zero. In this - indexing mode, elements correspond to the normal indexing mode where - the rightmost index increases most frequently. For example, when - using double-square bracket indexing on a two-dimensional matrix, - increasing indexes will reference the matrix elements left to right, - row by row. Thus in the following example, 'x' and 'y' are copied - from the same matrix element: +TYPES + index-range-list range_1 [, range_2, ...] up to 4 ranges + range_1, ... integer, or integer_1 : integer_2 + value, value_1, ... any + variable_1 ... lvalue + decl declarator = global, static or local + id_1, ... identifier - mat m[1:2, 1:3]; - x = m[2,1]; - y = m[[3]]; +DESCRIPTION + The expression mat [index-range-list] returns a matrix value. + This may be assigned to one or more lvalues A, B, ... by either - There are functions which return information about a matrix. - The 'size' functions returns the total number of elements. - The 'matdim', 'matmin', and 'matmax' functions return the number - of dimensions of a matrix, and the lower and upper index bounds - for a dimension of a matrix. For square matrices, the 'det' - function calculates the determinant of the matrix. + mat A B ... [index-range-list] - Some functions return matrices as their results. These functions - do not affect the original matrix argument, but instead return - new matrices. For example, the 'mattrans' function returns the - transpose of a matrix, and 'inverse' returns the inverse of a - matrix. So to invert a matrix called 'x', you could use: + or - x = inverse(x); + A = B = ... = mat[index-range-list] - The 'matfill' function fills all elements of a matrix with the - specified value, and optionally fills the diagonal elements of a - square matrix with a different value. For example: + If a variable is specified by an expression that is not a symbol with + possibly object element specifiers, the expression should be enclosed + in parentheses. For example, parentheses are required in + mat (A[2]) [3] and mat (*p) [3] but mat P.x [3] is acceptable. - matfill(x,1); + When an index-range is specified as integer_1 : integer_2, where + integer_1 and integer_2 are expressions which evaluate to integers, + the index-range consists of all integers from the minimum of the + two integers to the maximum of the two integers. For example, + mat[2:5, 0:4] and mat[5:2, 4:0] return the same matrix value. - will fill any matrix with ones, and: + If an index-range is an expression which evaluates to an integer, + the range is as if specified by 0 : integer - 1. For example, + mat[4] and mat[0:3] return the same 4-element matrix; mat[-2] and + mat[-3:0] return the same 4-element matrix. - matfill(x, 0, 1); + If the variable A has a matrix value, then for integer indices + i_1, i_2, ..., equal in number to the number of ranges specified at + its creation, and such that each index is in the corresponding range, + the matrix element associated with those index list is given as an + lvalue by the expressions A[i_1, i_2, ...]. - will create an identity matrix out of any square matrix. Note that - unlike most matrix functions, this function does not return a matrix - value, but manipulates the matrix argument itself. + The elements of the matrix are stored internally as a linear array + in which locations are arranged in order of increasing indices. + For example, in order of location, the six element of A = mat [2,3] + are - Matrices can be multiplied by numbers, which multiplies each element - by the number. Matrices can also be negated, conjugated, shifted, - rounded, truncated, fractioned, and modulo'ed. Each of these - operations is applied to each element. + A[0,0], A[0,1], A[0,2], A[1,0], A[1,,1], A[1,2]. - Matrices can be added or multiplied together if the operation is - legal. Note that even if the dimensions of matrices are compatible, - operations can still fail because of mismatched lower bounds. The - lower bounds of two matrices must either match, or else one of them - must have a lower bound of zero. Thus the following code: + These elements may also be specified using the double-bracket operator + with a single integer index as in A[[0]], A[[1]], ..., A[[5]]. + If p is assigned the value &A[0.0], the address of A[[i]] for 0 <= i < 6 + is p + i as long as A exists and a new value is not assigned to A. - mat x[3:3]; - mat y[4:4]; - z = x + y; + When a matrix is created, each element is initially assigned the + value zero. Other values may be assigned then or later using the + "= {...}" assignment operation. Thus - fails because the calculator does not have a way of knowing what - the bounds should be on the resulting matrix. If the bounds match, - then the resulting matrix has the same bounds. If exactly one of - the lower bounds is zero, then the resulting matrix will have the - nonzero lower bounds. Thus means that the bounds of a matrix are - preserved when operated on by matrices with lower bounds of zero. - For example: + A = {value_0, value_1, ...} - mat x[3:7]; - mat y[5]; - z = x + y; + assigns the values value_0, value_1, ... to the elements A[[0]], + A[[1]], ... Any blank "value" is passed over. For example, - will succeed and assign the variable 'z' a matrix whose - bounds are 3-7. + A = {1, , 2} - Vectors are matrices of only a single dimension. The 'dp' and 'cp' - functions calculate the dot product and cross product of a vector - (cross product is only defined for vectors of size 3). + will assign the value 1 to A[[0]], 2 to A[[2]] and leave all other + elements unchanged. Values may also be assigned to elements by + simple assignments, as in A[0,0] = 1, A[0,2] = 2; - Matrices can be searched for particular values by using the 'search' - and 'rsearch' functions. They return the element number of the - found value (zero based), or null if the value does not exist in the - matrix. Using the element number in double-bracket indexing will - then refer to the found element. + If the index-range is left blank but an initializer list is specified + as in + + mat A[] = {1, 2 } + B = mat[] = {1, , 3, } + + the matrix created is one-dimensional. If the list contains a + positive number n of values or blanks, the result is as if the + range were specified by [n], i.e. the range of indices is from + 0 to n - 1. In the above examples, A is of size 2 with A[0] = 1 + and A[1] = 2; B is of size 4 with B[0] = 1, B[1] = B[3] = 0, + B[2] = 3. The specification mat[] = { } creates the same as mat[1]. + + If the index-range is left blank and no initializer list is specified, + as in mat C[] or C = mat[], the matrix assigned to C has zero + dimension; this has one element C[]. To assign a value using "= { ...}" + at the same time as creating C, parentheses are required as in + (mat[]) = {value} or (mat C[]) = {value}. Later a value may be + assigned to C[] by C[] = value or C = {value}. + + The value assigned at any time to any element of a matrix can be of + any type - number, string, list, matrix, object of previously specified + type, etc. For some matrix operations there are of course conditions + that elements may have to satisfy: for example, addition of matrices + requires that addition of corresponding elements be possible. + If an element of a matrix is a structure for which indices or an + object element specifier is required, an element of that structure is + referred to by appropriate uses of [ ] or ., and so on if an element + of that element is required. For example, one may have an expressions + like + + A[1,2][3].alpha[2]; + + if A[1,2][3].alpha is a list with at least three elements, A[1,2][3] is + an object of a type like obj {alpha, beta}, A[1,2] is a matrix of + type mat[4] and A is a mat[2,3] matrix. When an element of a matrix + is a matrix and the total number of indices does not exceed 4, the + indices can be combined into one list, e.g. the A[1,2][3] in the + above example can be shortened to A[1,2,3]. (Unlike C, A[1,2] cannot + be expressed as A[1][2].) + + The function ismat(V) returns 1 if V is a matrix, 0 otherwise. + + isident(V) returns 1 if V is a square matrix with diagonal elements 1, + off-diagonal elements zero, or a zero- or one-dimensional matrix with + every element 1; otherwise zero is returned. Thus isident(V) = 1 + indicates that for V * A and A * V where A is any matrix of + for which either product is defined and the elements of A are real + or complex numbers, that product will equal A. + + If V is matrix-valued, test(V) returns 0 if every element of V tests + as zero; otherwise 1 is returned. + + The dimension of a matrix A, i.e. the number of index-ranges in the + initial creation of the matrix, is returned by the function matdim(A). + For 1 <= i <= matdim(A), the minimum and maximum values for the i-th + index range are returned by matmin(A, i) and matmax(A,i), respectively. + The total number of elements in the matrix is returned by size(A). + The sum of the elements in the matrix is returned by matsum(A). + + The default method of printing matrices is to give a line of information + about the matrix, and to list on separate lines up to 15 elements, + the indices and either the value (for numbers, strings, objects) or + some descriptive information for lists or matrices, etc. + Numbers are displayed in the current number-printing mode. + The maximum number of elements to be printed can be assigned + any nonnegative integer value m by config("maxprint", m). + + Users may define another method of printing matrices by defining a + function mat_print(M); for example, for a not too big 2-dimensional + matrix A it is a common practice to use a loop like: + + for (i = matmin(A,1); i <= matmax(A,1); i++) { + for (j = matmin(A,2); j <= matmax(A,2); j++) + printf("%8d", A[i,j]; + print; + } + + The default printing may be restored by + + undefine mat_print; + + + The keyword "mat" followed by two or more index-range-lists returns a + matrix with indices specified by the first list, whose elements are + matrices as determined by the later index-range-lists. For + example mat[2][3] is a 2-element matrix, each of whose elements has + as its value a 3-element matrix. Values may be assigned to the + elements of the innermost matrices by nested = {...} operations as in + + mat [2][3] = {{1,2,3},{4,5,6}} + + An example of the use of mat with a declarator is + + global mat A B [2,3], C [4] + + This creates, if they do not already exist, three global variables with + names A, B, C, and assigns to A and B the value mat[2,3] and to C mat[4]. + + Some operations are defined for matrices. + + A == B + Returns 1 if A and B are of the same "shape" and "corresponding" + elements are equal; otherwise 0 is returned. Being of the same + shape means they have the same dimension d, and for each i <= d, + + matmax(A,i) - matmin(A,i) == matmax(B,i) - matmin(B,i), + + One consequence of being the same shape is that the matrices will + have the same size. Elements "correspond" if they have the same + double-bracket indices; thus A == B implies that A[[i]] == B[[i]] + for 0 <= i < size(A) == size(B). + + A + B + A - B + These are defined A and B have the same shape, the element + with double-bracket index j being evaluated by A[[j]] + B[[j]] and + A[[j]] - B[[j]], respectively. The index-ranges for the results + are those for the matrix A. + + A[i,j] + If A is two-dimensional, it is customary to speak of the indices + i, j in A[i,j] as referring to rows and columns; the number of + rows is matmax(A,1) - matmin(A,1) + 1; the number of columns if + matmax(A,2) - matmin(A,2) + 1. A matrix is said to be square + if it is two-dimensional and the number of rows is equal to the + number of columns. + + A * B + Multiplication is defined provided certain conditions by the + dimensions and shapes of A and B are satisfied. If both have + dimension 2 and the column-index-list for A is the same as + the row-index-list for B, C = A * B is defined in the usual + way so that for i in the row-index-list of A and j in the + column-index-list for B, + + C[i,j] = Sum A[i,k] * B[k,j] + + the sum being over k in the column-index-list of A. The same + formula is used so long as the number of columns in A is the same + as the number of rows in B and k is taken to refer to the offset + from matmin(A,2) and matmin(B,1), respectively, for A and B. + If the multiplications and additions required cannot be performed, + an execution error may occur or the result for C may contain + one or more error-values as elements. + + If A or B has dimension zero, the result for A * B is simply + that of multiplying the elements of the other matrix on the + left by A[] or on the right by B[]. + + If both A and B have dimension 1, A * B is defined if A and B + have the same size; the result has the same index-list as A + and each element is the product of corresponding elements of + A and B. If A and B have the same index-list, this multiplication + is consistent with multiplication of 2D matrices if A and B are + taken to represent 2D matrices for which the off-diagonal elements + are zero and the diagonal elements are those of A and B. + the real and complex numbers. + + If A is of dimension 1 and B is of dimension 2, A * B is defined + if the number of rows in B is the same as the size of A. The + result has the same index-lists as B; each row of B is multiplied + on the left by the corresponding element of A. + + If A is of dimension 2 and B is of dimension 1, A * B is defined + if number of columns in A is the same as the size of A. The + result has the same index-lists as A; each column of A is + multiplied on the right by the corresponding element of B. + + The algebra of additions and multiplications involving both one- + and two-dimensional matrices is particularly simple when all the + elements are real or complex numbers and all the index-lists are + the same, as occurs, for example, if for some positive integer n, + all the matrices start as mat [n] or mat [n,n]. + + det(A) + If A is a square, det(A) is evaluated by an algorithm that returns + the determinant of A if the elements of A are real or complex + numbers, and if such an A is non-singular, inverse(A) returns + the inverse of A indexed in the same way as A. For matrix A of + dimension 0 or 1, det(A) is defined as the product of the elements + of A in the order in which they occur in A, inverse(A) returns + a matrix indexed in the same way as A with each element inverted. + + + The following functions are defined to return matrices with the same + index-ranges as A and the specified operations performed on all + elements of A. Here num is an arbitrary complex number (nonzero + when it is a divisor), int an integer, rnd a rounding-type + specifier integer, real a real number. + + num * A + A * num + A / num + - A + conj(A) + A << int, A >> int + scale(A, int) + round(A, int, rnd) + bround(A, int, rnd) + appr(A, real, rnd) + int(A) + frac(A) + A // real + A % real + A ^ int + + If A and B are one-dimensional of the same size dp(A, B) returns + their dot-product, i.e. the sum of the products of corresponding + elements. + + If A and B are one-dimension and of size 3, cp(A, B) returns their + cross-product. + + randperm(A) returns a matrix indexed the same as A in which the elements + of A have been randomly permuted. + + sort(A) returns a matrix indexed the same as A in which the elements + of A have been sorted. + + If A is an lvalue whose current value is a matrix, matfill(A, v) + assigns the value v to every element of A, and if also, A is + square, matfill(A, v1, v2) assigns v1 to the off-diagonal elements, + v2 to the diagonal elements. To create and assign to A the unit + n * n matrix, one may use matfill(mat A[n,n], 0, 1). + + For a square matrix A, mattrace(A) returns the trace of A, i.e. the + sum of the diagonal elements. For zero- or one-dimensional A, + mattrace(A) returns the sum of the elements of A. + + For a two-dimensional matrix A, mattrans(A) returns the transpose + of A, i.e. if A is mat[m,n], it returns a mat[n,m] matrix with + [i,j] element equal to A[j,i]. For zero- or one-dimensional A, + mattrace(A) returns a matrix with the same value as A. + + The functions search(A, value, start, end]) and + rsearch(A, value, start, end]) return the first or last index i + for which A[[i]] == value and start <= i < end, or if there is + no such index, the null value. For further information on default + values and the use of an "accept" function, see the help files for + search and rsearch. + + reverse(A) returns a matrix with the same index-lists as A but the + elements in reversed order. + + The copy and blkcpy functions may be used to copy data to a matrix from + a matrix or list, or from a matrix to a list. In copying from a + matrix to a matrix the matrices need not have the same dimension; + in effect they are treated as linear arrays. + +EXAMPLE + > obj point {x,y} + > mat A[5] = {1, 2+3i, "ab", mat[2] = {4,5}. obj point = {6,7}} + > A + mat [5] (5 elements, 5 nonzero): + [0] = 1 + [1] = 2+3i + [2] = "ab" + [3] = mat [2] (2 elements, 2 nonzero) + [4] = obj point {6, 7} + + > print A[0], A[1], A[2], A[3][0], A[4].x + 1 2+3i ab 4 6 + + > define point_add(a,b) = obj point = {a.x + b.x, a.y + b.y} + point_add(a,b) defined + + > mat [B] = {8, , "cd", mat[2] = {9,10}, obj point = {11,12}} + > A + B + + mat [5] (5 elements, 5 nonzero): + [0] = 9 + [1] = 2+3i + [2] = "abcd" + [3] = mat [2] (2 elements, 2 nonzero) + [4] = obj point {17, 19} + + > mat C[2,2] = {1,2,3,4} + > C^10 + + mat [2,2] (4 elements, 4 nonzero): + [0,0] = 4783807 + [0,1] = 6972050 + [1,0] = 10458075 + [1,1] = 15241882 + + > C^-10 + + mat [2,2] (4 elements, 4 nonzero): + [0,0] = 14884.650390625 + [0,1] = -6808.642578125 + [1,0] = -10212.9638671875 + [1,1] = 4671.6865234375 + + > mat A[4] = {1,2,3,4}, A * reverse(A); + + mat [4] (4 elements, 4 nonzero): + [0] = 4 + [1] = 6 + [2] = 6 + [3] = 4 + +LIMITS + The theoretical upper bound for the absolute values of indices is + 2^31 - 1, but the size of matrices that can be handled in practice will + be limited by the availability of memory and what is an acceptable + runtime. For example, although it may take only a fraction of a + second to invert a 10 * 10 matrix, it will probably take about 1000 + times as long to invert a 100 * 100 matrix. + +LIBRARY + n/a + +SEE ALSO + ismat, matdim, matmax, matmin, mattrans, mattrace, matsum, det, inverse, + isident, test, config, search, rsearch, reverse, copy, blkcpy, dp, cp, + randperm, sort diff --git a/help/todo b/help/todo index 87308f5..4efd6dd 100644 --- a/help/todo +++ b/help/todo @@ -48,7 +48,6 @@ Very High priority items: history command history interrupt how interrupts are handled list using lists - mat using matrices obj user defined data types operator math, relational, logic and variable access ... statement flow control and declaration statements diff --git a/help/wishlist b/help/wishlist index 9a671ef..1cc8ab0 100644 --- a/help/wishlist +++ b/help/wishlist @@ -11,7 +11,7 @@ Calc Enhancement Wish List: The following items are in the calc wish list. Programs like this can be extended and improved forever. - Calc bug repoers, however, should be sent to: + Calc bug reports, however, should be sent to: calc-bugs at postofc dot corp dot sgi dot com diff --git a/matfunc.c b/matfunc.c index 5243e7b..4ab8e0c 100644 --- a/matfunc.c +++ b/matfunc.c @@ -44,7 +44,8 @@ matadd(MATRIX *m1, MATRIX *m2) max1 = m1->m_max[dim]; min2 = m2->m_min[dim]; max2 = m2->m_max[dim]; - if ((min1 && min2 && (min1 != min2)) || ((max1-min1) != (max2-min2))) { + if ((min1 && min2 && (min1 != min2)) || + ((max1-min1) != (max2-min2))) { math_error("Incompatible matrix bounds for add"); /*NOTREACHED*/ } @@ -85,7 +86,8 @@ matsub(MATRIX *m1, MATRIX *m2) max1 = m1->m_max[dim]; min2 = m2->m_min[dim]; max2 = m2->m_max[dim]; - if ((min1 && min2 && (min1 != min2)) || ((max1-min1) != (max2-min2))) { + if ((min1 && min2 && (min1 != min2)) || + ((max1-min1) != (max2-min2))) { math_error("Incompatible matrix bounds for sub"); /*NOTREACHED*/ } @@ -131,15 +133,89 @@ matmul(MATRIX *m1, MATRIX *m2) { register MATRIX *res; long i1, i2, max1, max2, index, maxindex; - VALUE *v1, *v2; + VALUE *v1, *v2, *vres; VALUE sum, tmp1, tmp2; + if (m1->m_dim == 0) { + i2 = m2->m_size; + v2 = m2->m_table; + res = matalloc(i2); + *res = *m2; + vres = res->m_table; + while (i2-- > 0) + mulvalue(m1->m_table, v2++, vres++); + return res; + } + if (m2->m_dim == 0) { + i1 = m1->m_size; + v1 = m1->m_table; + res = matalloc(i1); + *res = *m1; + vres = res->m_table; + while (i1-- > 0) + mulvalue(v1++, m2->m_table, vres++); + return res; + } + if (m1->m_dim == 1 && m2->m_dim == 1) { + if (m1->m_max[0]-m1->m_min[0] != m2->m_max[0]-m2->m_min[0]) { + math_error("Incompatible bounds for 1D * 1D matmul"); + /*NOTREACHED*/ + } + res = matalloc(m1->m_size); + *res = *m1; + v1 = m1->m_table; + v2 = m2->m_table; + vres = res->m_table; + for (index = m1->m_size; index > 0; index--) + mulvalue(v1++, v2++, vres++); + return res; + } + if (m1->m_dim == 1 && m2->m_dim == 2) { + if (m1->m_max[0]-m1->m_min[0] != m2->m_max[0]-m2->m_min[0]) { + math_error("Incompatible bounds for 1D * 2D matmul"); + /*NOTREACHED*/ + } + res = matalloc(m2->m_size); + *res = *m2; + i1 = m1->m_max[0] - m1->m_min[0] + 1; + max2 = m2->m_max[1] - m2->m_min[1] + 1; + v1 = m1->m_table; + v2 = m2->m_table; + vres = res->m_table; + while (i1-- > 0) { + i2 = max2; + while (i2-- > 0) + mulvalue(v1, v2++, vres++); + v1++; + } + return res; + } + if (m1->m_dim == 2 && m2->m_dim == 1) { + if (m1->m_max[1]-m1->m_min[1] != m2->m_max[0]-m2->m_min[0]) { + math_error("Incompatible bounds for 2D * 1D matmul"); + /*NOTREACHED*/ + } + res = matalloc(m1->m_size); + *res = *m1; + i1 = m1->m_max[0] - m1->m_min[0] + 1; + max1 = m1->m_max[1] - m1->m_min[1] + 1; + v1 = m1->m_table; + vres = res->m_table; + while (i1-- > 0) { + v2 = m2->m_table; + i2 = max1; + while (i2-- > 0) + mulvalue(v1++, v2++, vres++); + } + return res; + } + if ((m1->m_dim != 2) || (m2->m_dim != 2)) { - math_error("Matrix dimension must be two for mul"); + math_error("Matrix dimensions not compatible for mul"); /*NOTREACHED*/ } - if ((m1->m_max[1] - m1->m_min[1]) != (m2->m_max[0] - m2->m_min[0])) { - math_error("Incompatible bounds for matrix mul"); + if ((m1->m_max[1]-m1->m_min[1]) != (m2->m_max[0]-m2->m_min[0])) { + math_error("Incompatible bounds for 2D * 2D matrix mul"); /*NOTREACHED*/ } max1 = (m1->m_max[0] - m1->m_min[0] + 1); @@ -164,7 +240,8 @@ matmul(MATRIX *m1, MATRIX *m2) freevalue(&sum); sum = tmp2; v1++; - v2 += max2; + if (index+1 < maxindex) + v2 += max2; } index = (i1 * max2) + i2; res->m_table[index] = sum; @@ -185,8 +262,17 @@ matsquare(MATRIX *m) VALUE *v1, *v2; VALUE sum, tmp1, tmp2; + if (m->m_dim < 2) { + res = matalloc(m->m_size); + *res = *m; + v1 = m->m_table; + v2 = res->m_table; + for (index = m->m_size; index > 0; index--) + squarevalue(v1++, v2++); + return res; + } if (m->m_dim != 2) { - math_error("Matrix dimension must be two for square"); + math_error("Matrix dimension exceeds two for square"); /*NOTREACHED*/ } if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) { @@ -224,7 +310,8 @@ matsquare(MATRIX *m) /* - * Compute the result of raising a square matrix to an integer power. + * Compute the result of raising a matrix to an integer power if + * dimension <= 2 and for dimension == 2, the matrix is square. * Negative powers mean the positive power of the inverse. * Note: This calculation could someday be improved for large powers * by using the characteristic polynomial of the matrix. @@ -240,12 +327,13 @@ matpowi(MATRIX *m, NUMBER *q) long power; /* power to raise to */ FULL bit; /* current bit value */ - if (m->m_dim != 2) { - math_error("Matrix dimension must be two for power"); + if (m->m_dim > 2) { + math_error("Matrix dimension greater than 2 for power"); /*NOTREACHED*/ } - if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) { - math_error("Raising non-square matrix to a power"); + if (m->m_dim == 2 && (m->m_max[0] - m->m_min[0] != + m->m_max[1] - m->m_min[1])) { + math_error("Raising non-square 2D matrix to a power"); /*NOTREACHED*/ } if (qisfrac(q)) { @@ -544,6 +632,10 @@ mattrace(MATRIX *m) VALUE tmp; long i, j; + if (m->m_dim < 2) { + matsum(m, &sum); + return sum; + } if (m->m_dim != 2) return error_value(E_MATTRACE2); i = (m->m_max[0] - m->m_min[0] + 1); @@ -574,6 +666,8 @@ mattrans(MATRIX *m) long row, col; /* current row and column */ MATRIX *res; + if (m->m_dim < 2) + return matcopy(m); res = matalloc(m->m_size); res->m_dim = 2; res->m_min[0] = m->m_min[1]; @@ -742,13 +836,15 @@ matindex(MATRIX *mp, BOOL create, long dim, VALUE *indices) long offset; /* current offset into array */ int i; /* loop counter */ - if (dim <= 0) { - math_error("Bad dimension %ld for matrix", dim); + if (dim < 0) { + math_error("Negative dimension %ld for matrix", dim); /*NOTREACHED*/ } for (;;) { if (dim < mp->m_dim) { - math_error("Indexing a %ldd matrix as a %ldd matrix", mp->m_dim, dim); + math_error( + "Indexing a %ldd matrix as a %ldd matrix", + mp->m_dim, dim); /*NOTREACHED*/ } offset = 0; @@ -763,7 +859,8 @@ matindex(MATRIX *mp, BOOL create, long dim, VALUE *indices) /*NOTREACHED*/ } index = qtoi(q); - if (zge31b(q->num) || (index < mp->m_min[i]) || (index > mp->m_max[i])) { + if (zge31b(q->num) || (index < mp->m_min[i]) || + (index > mp->m_max[i])) { math_error("Index out of bounds for matrix"); /*NOTREACHED*/ } @@ -785,6 +882,36 @@ matindex(MATRIX *mp, BOOL create, long dim, VALUE *indices) } +/* + * Returns the list of indices for a matrix element with specified + * double-bracket index. + */ +LIST * +matindices(MATRIX *mp, long index) +{ + LIST *lp; + int j; + long d; + VALUE val; + + if (index < 0 || index >= mp->m_size) + return NULL; + + lp = listalloc(); + val.v_type = V_NUM; + j = mp->m_dim; + + while (--j >= 0) { + d = mp->m_max[j] - mp->m_min[j] + 1; + val.v_num = itoq(index % d + mp->m_min[j]); + insertlistfirst(lp, &val); + qfree(val.v_num); + index /= d; + } + return lp; +} + + /* * Search a matrix for the specified value, starting with the specified index. * Returns 0 and stores index if value found; otherwise returns 1. @@ -898,7 +1025,8 @@ matident(MATRIX *m) MATRIX *res; /* resulting matrix */ if (m->m_dim != 2) { - math_error("Matrix dimension must be two for setting to identity"); + math_error( + "Matrix dimension must be two for setting to identity"); /*NOTREACHED*/ } if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) { @@ -912,7 +1040,8 @@ matident(MATRIX *m) for (row = 0; row < rows; row++) { for (col = 0; col < rows; col++) { val->v_type = V_NUM; - val->v_num = ((row == col) ? qlink(&_qone_) : qlink(&_qzero_)); + val->v_num = ((row == col) ? qlink(&_qone_) : + qlink(&_qzero_)); val++; } } @@ -934,11 +1063,21 @@ matinv(MATRIX *m) long cur; /* current row being worked on */ long row, col; /* temp row and column values */ VALUE *val; /* current value in matrix*/ + VALUE *vres; /* current value in result for dim < 2 */ VALUE mulval; /* value to multiply rows by */ VALUE tmpval; /* temporary value */ + if (m->m_dim < 2) { + res = matalloc(m->m_size); + *res = *m; + val = m->m_table; + vres = res->m_table; + for (cur = m->m_size; cur > 0; cur--) + invertvalue(val++, vres++); + return res; + } if (m->m_dim != 2) { - math_error("Matrix dimension must be two for inverse"); + math_error("Matrix dimension exceeds two for inverse"); /*NOTREACHED*/ } if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) { @@ -994,26 +1133,30 @@ matinv(MATRIX *m) matswaprow(res, row, cur); } /* - * Now for every other nonzero entry in the current column, subtract - * the appropriate multiple of the current row to force that entry - * to become zero. + * Now for every other nonzero entry in the current column, + * subtract the appropriate multiple of the current row to + * force that entry to become zero. */ val = &m->m_table[cur]; - /* ignore Saber-C warning #26 - storing bad pointer in val */ - /* ok to ignore on name matinv`val */ - for (row = 0; row < rows; row++, val += rows) { - if ((row == cur) || (testvalue(val) == 0)) + for (row = 0; row < rows; row++) { + if ((row == cur) || (testvalue(val) == 0)) { + if (row+1 < rows) + val += rows; continue; + } mulvalue(val, &mulval, &tmpval); matsubrow(m, row, cur, &tmpval); matsubrow(res, row, cur, &tmpval); freevalue(&tmpval); + if (row+1 < rows) + val += rows; } freevalue(&mulval); } /* - * Now the original matrix has nonzero entries only on its main diagonal. - * Scale the rows of the result matrix by the inverse of those entries. + * Now the original matrix has nonzero entries only on its main + * diagonal. Scale the rows of the result matrix by the inverse + * of those entries. */ val = m->m_table; for (row = 0; row < rows; row++) { @@ -1022,7 +1165,8 @@ matinv(MATRIX *m) matmulrow(res, row, &mulval); freevalue(&mulval); } - val += (rows + 1); + if (row+1 < rows) + val += (rows + 1); } matfree(m); return res; @@ -1044,6 +1188,24 @@ matdet(MATRIX *m) VALUE tmp1, tmp2, tmp3; BOOL neg; /* whether to negate determinant */ + if (m->m_dim < 2) { + vp = m->m_table; + i = m->m_size; + copyvalue(vp, &tmp1); + + while (--i > 0) { + mulvalue(&tmp1, ++vp, &tmp2); + freevalue(&tmp1); + tmp1 = tmp2; + } + return tmp1; + } + + if (m->m_dim != 2) + return error_value(E_DET2); + if ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1])) + return error_value(E_DET3); + /* * Loop over each row, and eliminate all lower entries in the * corresponding column by using row operations. Copy the original @@ -1113,7 +1275,8 @@ matdet(MATRIX *m) } } div = pivot; - pivot += n + 1; + if (k > 0) + pivot += n + 1; } if (neg) negvalue(div, &tmp1); @@ -1298,7 +1461,8 @@ matalloc(long size) m = (MATRIX *) malloc(matsize(size)); if (m == NULL) { - math_error("Cannot get memory to allocate matrix of size %d", size); + math_error("Cannot get memory to allocate matrix of size %d", + size); /*NOTREACHED*/ } m->m_size = size; @@ -1384,7 +1548,8 @@ matcmp(MATRIX *m1, MATRIX *m2) if ((m1->m_dim != m2->m_dim) || (m1->m_size != m2->m_size)) return TRUE; for (i = 0; i < m1->m_dim; i++) { - if ((m1->m_max[i] - m1->m_min[i]) != (m2->m_max[i] - m2->m_min[i])) + if ((m1->m_max[i] - m1->m_min[i]) != + (m2->m_max[i] - m2->m_min[i])) return TRUE; } v1 = m1->m_table; @@ -1509,10 +1674,20 @@ matisident(MATRIX *m) register VALUE *val; /* current value */ long row, col; /* row and column numbers */ + val = m->m_table; + if (m->m_dim == 0) { + return (val->v_type == V_NUM && qisone(val->v_num)); + } + if (m->m_dim == 1) { + for (row = m->m_min[0]; row <= m->m_max[0]; row++, val++) { + if (val->v_type != V_NUM || !qisone(val->v_num)) + return FALSE; + } + return TRUE; + } if ((m->m_dim != 2) || ((m->m_max[0] - m->m_min[0]) != (m->m_max[1] - m->m_min[1]))) return FALSE; - val = m->m_table; for (row = m->m_min[0]; row <= m->m_max[0]; row++) { /* * We could use col = m->m_min[1]; col < m->m_max[1] @@ -1558,15 +1733,22 @@ matprint(MATRIX *m, long max_print) fullsize *= (m->m_max[i] - m->m_min[i] + 1); } msg = ((max_print > 0) ? "\nmat [" : "mat ["); - for (i = 0; i < dim; i++) { - if (m->m_min[i]) - math_fmt("%s%ld:%ld", msg, m->m_min[i], m->m_max[i]); - else - math_fmt("%s%ld", msg, m->m_max[i] + 1); - msg = ","; + if (dim) { + for (i = 0; i < dim; i++) { + if (m->m_min[i]) { + math_fmt("%s%ld:%ld", msg, + m->m_min[i], m->m_max[i]); + } else { + math_fmt("%s%ld", msg, m->m_max[i] + 1); + } + msg = ","; + } + } else { + math_str("mat ["); } - if (max_print > fullsize) + if (max_print > fullsize) { max_print = fullsize; + } vp = m->m_table; count = 0; for (index = 0; index < fullsize; index++) { @@ -1588,10 +1770,15 @@ matprint(MATRIX *m, long max_print) for (index = 0; index < max_print; index++) { msg = " ["; j = index; - for (i = 0; i < dim; i++) { - math_fmt("%s%ld", msg, m->m_min[i] + (j / sizes[i])); - j %= sizes[i]; - msg = ","; + if (dim) { + for (i = 0; i < dim; i++) { + math_fmt("%s%ld", msg, + m->m_min[i] + (j / sizes[i])); + j %= sizes[i]; + msg = ","; + } + } else { + math_str(msg); } math_str("] = "); printvalue(vp++, PRINT_SHORT | PRINT_UNAMBIG); diff --git a/opcodes.c b/opcodes.c index 7d3a97e..c5c5ece 100644 --- a/opcodes.c +++ b/opcodes.c @@ -320,7 +320,7 @@ o_matcreate(FUNC *fp, long dim) long tmp; /* temporary */ long size; /* size of matrix */ - if ((dim <= 0) || (dim > MAXDIM)) { + if ((dim < 0) || (dim > MAXDIM)) { math_error("Bad dimension %ld for matrix", dim); /*NOTREACHED*/ } @@ -489,8 +489,8 @@ o_indexaddr(FUNC *fp, long dim, long writeflag) BLOCK *blk; flag = (writeflag != 0); - if (dim <= 0) { - math_error("Zero or negative dimensions for indexing"); + if (dim < 0) { + math_error("Negative dimension for indexing"); /*NOTREACHED*/ } val = &stack[-dim]; diff --git a/qmod.c b/qmod.c index c46c1b0..70d7d9f 100644 --- a/qmod.c +++ b/qmod.c @@ -327,7 +327,7 @@ qfindredc(NUMBER *q) /* * First try for an exact pointer match in the table. */ - for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) { + for (rcp = redc_cache; rcp <= &redc_cache[MAXREDC-1]; rcp++) { if (q == rcp->rnum) { rcp->age = ++redc_age; return rcp->redc; @@ -337,7 +337,7 @@ qfindredc(NUMBER *q) /* * Search the table again looking for a value which matches. */ - for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) { + for (rcp = redc_cache; rcp <= &redc_cache[MAXREDC-1]; rcp++) { if (rcp->age && (qcmp(q, rcp->rnum) == 0)) { rcp->age = ++redc_age; return rcp->redc; @@ -355,7 +355,7 @@ qfindredc(NUMBER *q) } bestrcp = NULL; - for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) { + for (rcp = redc_cache; rcp <= &redc_cache[MAXREDC-1]; rcp++) { if ((bestrcp == NULL) || (rcp->age < bestrcp->age)) bestrcp = rcp; } diff --git a/quickhash.c b/quickhash.c index 075b402..9526c6d 100644 --- a/quickhash.c +++ b/quickhash.c @@ -271,12 +271,14 @@ mathash(MATRIX *m, QCKHASH val) * hash 10 more elements if they exist */ i = 16; - vp = &m->m_table[16]; - skip = (m->m_size / 11) + 1; - while (i < m->m_size) { - val = hashvalue(vp, val); - i += skip; - vp += skip; + if (i < m->m_size) { + vp = (VALUE *)&m->m_table[i]; + skip = (m->m_size / 11) + 1; + while (i < m->m_size) { + val = hashvalue(vp, val); + i += skip; + vp += skip; + } } return val; } diff --git a/string.c b/string.c index 89d9afe..97433cc 100644 --- a/string.c +++ b/string.c @@ -584,8 +584,20 @@ stringsegment(STRING *s1, long n1, long n2) s->s_str = c; c1 = s1->s_str + n1; if (n1 >= n2) { - while (len-- > 0) + /* + * We prevent the c1 pointer from walking behind s1_s_str + * by stopping one short of the end and running the loop one + * more time. + * + * We could stop the loop with just len-- > 0, but stopping + * short and running the loop one last time manually helps make + * code checkers such as insure happy. + */ + while (len-- > 1) { *c++ = *c1--; + } + /* run the loop manually one last time */ + *c++ = *c1; } else { while (len-- > 0) *c++ = *c1++; diff --git a/version.c b/version.c index e52869b..3fd4ad7 100644 --- a/version.c +++ b/version.c @@ -18,7 +18,7 @@ static char *program; #define MAJOR_VER 2 /* major version */ #define MINOR_VER 11 /* minor version */ #define MAJOR_PATCH 0 /* patch level or 0 if no patch */ -#define MINOR_PATCH "10.2" /* test number or empty string if no patch */ +#define MINOR_PATCH "10.3" /* test number or empty string if no patch */ /* * calc version constants diff --git a/zio.c b/zio.c index 9e23ff1..fad1747 100644 --- a/zio.c +++ b/zio.c @@ -394,7 +394,7 @@ zprintb(ZVALUE z, long width) didprint = 0; PUTSTR("0b"); while (len-- >= 0) { - val = *hp--; + val = ((len >= 0) ? *hp-- : *hp); mask = ((HALF)1 << (BASEB - 1)); while (mask) { ch = '0' + ((mask & val) != 0); @@ -481,15 +481,17 @@ zprinto(ZVALUE z, long width) break; } len -= rem; - hp -= rem; - while (len > 0) { /* finish in groups of 3 words */ - PRINTF4("%08lo%08lo%08lo%08lo", - (PRINT) ((hp[0]) >> 8), - (PRINT) (((hp[0] & 0xff) << 16) + (hp[-1] >> 16)), - (PRINT) (((hp[-1] & 0xffff) << 8) + (hp[-2] >> 24)), - (PRINT) (hp[-2] & 0xffffff)); - hp -= 3; - len -= 3; + if (len > 0) { + hp -= rem; + while (len > 0) { /* finish in groups of 3 words */ + PRINTF4("%08lo%08lo%08lo%08lo", + (PRINT) ((hp[0]) >> 8), + (PRINT) (((hp[0] & 0xff) << 16) + (hp[-1] >> 16)), + (PRINT) (((hp[-1] & 0xffff) << 8) + (hp[-2] >> 24)), + (PRINT) (hp[-2] & 0xffffff)); + hp -= 3; + len -= 3; + } } #else switch (rem) { /* handle odd amounts first */ @@ -513,13 +515,15 @@ zprinto(ZVALUE z, long width) PRINTF1("0%lo", num2); } len -= rem; - hp -= rem; - while (len > 0) { /* finish in groups of 3 halfwords */ - PRINTF2("%08lo%08lo", - ((((FULL) hp[0]) << 8) + (((FULL) hp[-1]) >> 8)), - ((((FULL) (hp[-1] & 0xff)) << 16) + ((FULL) hp[-2]))); - hp -= 3; - len -= 3; + if (len > 0) { + hp -= rem; + while (len > 0) { /* finish in groups of 3 halfwords */ + PRINTF2("%08lo%08lo", + ((((FULL) hp[0]) << 8) + (((FULL) hp[-1]) >> 8)), + ((((FULL) (hp[-1] & 0xff))<<16) + ((FULL) hp[-2]))); + hp -= 3; + len -= 3; + } } #endif } diff --git a/zmath.c b/zmath.c index d028999..3e2a35b 100644 --- a/zmath.c +++ b/zmath.c @@ -1618,15 +1618,14 @@ ztest(ZVALUE z) /* - * Compare two numbers to see which is larger. - * Returns -1 if first number is smaller, 0 if they are equal, and 1 if - * first number is larger. This is the same result as ztest(z2-z1). + * Return the sign of z1 - z2, i.e. 1 if the first integer is greater, + * 0 if they are equal, -1 otherwise. */ FLAG zrel(ZVALUE z1, ZVALUE z2) { - register HALF *h1, *h2; - register FULL len1, len2; + HALF *h1, *h2; + LEN len; int sign; sign = 1; @@ -1636,66 +1635,47 @@ zrel(ZVALUE z1, ZVALUE z2) return -1; if (z2.sign) sign = -1; - len1 = z1.len; - len2 = z2.len; - h1 = z1.v + z1.len - 1; - h2 = z2.v + z2.len - 1; - while (len1 > len2) { - if (*h1--) - return sign; - len1--; - } - while (len2 > len1) { - if (*h2--) - return -sign; - len2--; - } - while (len1--) { - if (*h1-- != *h2--) + if (z1.len != z2.len) + return (z1.len > z2.len) ? sign : -sign; + len = z1.len; + h1 = z1.v + len; + h2 = z2.v + len; + + while (len > 0) { + if (*--h1 != *--h2) break; + len--; } - if ((len1 = *++h1) > (len2 = *++h2)) - return sign; - if (len1 < len2) - return -sign; + if (len > 0) + return (*h1 > *h2) ? sign : -sign; return 0; } /* - * Compare the absolute value two numbers to see which is larger. - * Returns -1 if first number is smaller, 0 if they are equal, and 1 if - * first number is larger. This is the same result as ztest(abs(z2)-abs(z1)) - * or zrel(abs(z1), abs(z2)). + * Return the sign of abs(z1) - abs(z2), i.e. 1 if the first integer + * has greater absolute value, 0 is they have equal absolute value, + * -1 otherwise. */ FLAG zabsrel(ZVALUE z1, ZVALUE z2) { - register HALF *h1, *h2; - register FULL len1, len2; + HALF *h1, *h2; + LEN len; - len1 = z1.len; - len2 = z2.len; - h1 = z1.v + z1.len - 1; - h2 = z2.v + z2.len - 1; - while (len1 > len2) { - if (*h1--) - return 1; - len1--; - } - while (len2 > len1) { - if (*h2--) - return -1; - len2--; - } - while (len1--) { - if (*h1-- != *h2--) + if (z1.len != z2.len) + return (z1.len > z2.len) ? 1 : -1; + + len = z1.len; + h1 = z1.v + len; + h2 = z2.v + len; + while (len > 0) { + if (*--h1 != *--h2) break; + len--; } - if ((len1 = *++h1) > (len2 = *++h2)) - return 1; - if (len1 < len2) - return -1; + if (len > 0) + return (*h1 > *h2) ? 1 : -1; return 0; } @@ -1715,8 +1695,8 @@ zcmp(ZVALUE z1, ZVALUE z2) len = z1.len; h1 = z1.v; h2 = z2.v; - while (len-- > 0) { - if (*h1++ != *h2++) + while (--len > 0) { + if (*++h1 != *++h2) return TRUE; } return FALSE; diff --git a/zmod.c b/zmod.c index a34fd74..72da090 100644 --- a/zmod.c +++ b/zmod.c @@ -543,7 +543,7 @@ zpowermod(ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res) } /* zzz */ - for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++) { + for (pp = &lowpowers[2]; pp <= &lowpowers[POWNUMS-1]; pp++) { pp->len = 0; pp->v = NULL; } @@ -558,16 +558,17 @@ zpowermod(ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res) curshift -= POWBITS; /* - * Calculate the result by examining the power POWBITS bits at a time, - * and use the table of low powers at each iteration. + * Calculate the result by examining the power POWBITS bits at + * a time, and use the table of low powers at each iteration. */ for (;;) { curpow = (curhalf >> curshift) & (POWNUMS - 1); pp = &lowpowers[curpow]; /* - * If the small power is not yet saved in the table, then - * calculate it and remember it in the table for future use. + * If the small power is not yet saved in the table, + * then calculate it and remember it in the table for + * future use. */ if (pp->v == NULL) { if (curpow & 0x1) @@ -575,10 +576,13 @@ zpowermod(ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res) else modpow = _one_; - for (curbit = 0x2; curbit <= curpow; curbit *= 2) { + for (curbit = 0x2; + curbit <= curpow; + curbit *= 2) { pp = &lowpowers[curbit]; if (pp->v == NULL) { - zsquare(lowpowers[curbit/2], &temp); + zsquare(lowpowers[curbit/2], + &temp); zmod5(&temp); zcopy(temp, pp); zfree(temp); @@ -599,8 +603,8 @@ zpowermod(ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res) } /* - * If the power is nonzero, then accumulate the small power - * into the result. + * If the power is nonzero, then accumulate the small + * power into the result. */ if (curpow) { zmul(ans, *pp, &temp); @@ -611,20 +615,20 @@ zpowermod(ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res) } /* - * Select the next POWBITS bits of the power, if there is - * any more to generate. + * Select the next POWBITS bits of the power, if + * there is any more to generate. */ curshift -= POWBITS; if (curshift < 0) { - if (hp-- == z2.v) + if (hp == z2.v) break; - curhalf = *hp; + curhalf = *--hp; curshift = BASEB - POWBITS; } /* - * Square the result POWBITS times to make room for the next - * chunk of bits. + * Square the result POWBITS times to make room for + * the next chunk of bits. */ for (i = 0; i < POWBITS; i++) { zsquare(ans, &temp); @@ -635,7 +639,7 @@ zpowermod(ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res) } } - for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++) { + for (pp = &lowpowers[2]; pp <= &lowpowers[POWNUMS-1]; pp++) { if (pp->v != NULL) freeh(pp->v); } @@ -669,7 +673,7 @@ zpowermod(ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res) * Modulus or power is small enough to perform the power raising * directly. Initialize the table of powers. */ - for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++) { + for (pp = &lowpowers[2]; pp <= &lowpowers[POWNUMS-1]; pp++) { pp->len = 0; pp->v = NULL; } @@ -757,7 +761,7 @@ zpowermod(ZVALUE z1, ZVALUE z2, ZVALUE z3, ZVALUE *res) } } - for (pp = &lowpowers[2]; pp < &lowpowers[POWNUMS]; pp++) { + for (pp = &lowpowers[2]; pp <= &lowpowers[POWNUMS-1]; pp++) { if (pp->v != NULL) freeh(pp->v); } diff --git a/zrand.c b/zrand.c index 9c78e1e..020ef2e 100644 --- a/zrand.c +++ b/zrand.c @@ -1724,7 +1724,7 @@ randcopy(CONST RAND *state) math_error("can't allocate RAND state"); /*NOTREACHED*/ } - *ret = *state; + memcpy(ret, state, sizeof(RAND)); /* * return copy