#include "m.h" void m_print( void) { Symbolptr s, a = pop(); check(a); if( a->name) s = a; else { if( (s = lookup("ans")) == NULL) s = install( "ans", UNDEF, EXTERN, 0, 0); s = copy( s, a); } print(s); } void m_mul( void) /* v(m,q) <- v1(m,n) * v2(p,q) */ { Symbolptr v1, v2; Symbolptr v; int m, n, p, q; v2 = pop(); v1 = pop(); m_check(v1,"*"); m_check(v2,"*"); if( SCALAR(v1) == SCALAR(v2)) { /* both scalars or matrices */ m = ROWS(v1); n = COLS(v1); p = ROWS(v2); q = COLS(v2); if( n != p) error("incompatible sizes for *",""); v = m_create( m, q); if( SCALAR(v1) && SCALAR(v2) ) VAL(v) = VAL(v1) * VAL(v2); else if( SCALAR(v)) VAL(v) = dot( MAT(v1), MAT(v2), n); else matmul( MAT(v), MAT(v1), MAT(v2), m, n, q); push(v); } else push( dyadic( v1, v2, '*') ); } void m_transpose( void) { Symbolptr a = pop(); Symbolptr v; m_check(a,"'"); if( a->class == CONST && COLS(a) == ROWS(a)) v = a; else v = m_create( COLS(a), ROWS(a)); if( SCALAR(a)) VAL(v) = VAL(a); else if( v == a) mattrans1( MAT(v), ROWS(v)); else mattrans( MAT(v), MAT(a), ROWS(a), COLS(a) ); push(v); } void m_unary( void) { Symbolptr a = pop(); int op = (int) *pc++; m_check(a,"unary operator"); push( unary( a, op) ); } void m_pbs( void) { Symbolptr v, a = pop(); double t; int i, j; m_check( a, "+\\"); if( a->class == CONST) v = a; else v = m_create( ROWS(a), COLS(a)); if( SCALAR(v)) VAL(v) = VAL(a); else for( i = 0; i < ROWS(a); ++i) { t = 0; for( j = 0; j < COLS(a); ++j) { t += A(i,j); V(i,j) = t; } } push(v); } void m_tbs( void) { Symbolptr v, a = pop(); double t; int i, j; m_check( a, "*\\"); if( a->class == CONST) v = a; else v = m_create( ROWS(a), COLS(a)); if( SCALAR(v)) VAL(v) = VAL(a); else for( i = 0; i < ROWS(a); ++i) { t = 1; for( j = 0; j < COLS(a); ++j) { t *= A(i,j); V(i,j) = t; } } push(v); } void m_reduce( void) { Symbolptr a = pop(); int op = (int) *pc++; m_check(a,"reduce operator"); push( reduce( a, op) ); } void m_copy( void) { Symbolptr a, b; a = pop(); b = pop(); a = copy( a, b); /* a <<-- b */ push(a); } void m_pp( void) /* a = a + d */ { Symbolptr a = pop(); int d = (int) *pc++; m_check(a,"++ or --"); c_check(a,"++ or --"); if( SCALAR(a) ) VAL(a) += d; else matsdyadic( MAT(a), MAT(a), (double)d, ROWS(a), COLS(a), '+'); push(a); } void m_pp2( void) /* ++a(b) or ++a(b,c) where increment is d */ { Symbolptr a, b = NULL, c = NULL; int narg, i, ia, ja, ib, jb, ic, jc; double d; d = (double) ((int) *pc++); narg = (int) *pc++; a = pop(); m_check(a,"++ or --"); c_check(a,"++ or --"); if( narg > 1) c = pop(); if( narg > 0) b = pop(); if( narg == 1 && ROWS(a) == 1) { c = b; b = NULL; } if( b == NULL) { b = m_create( 1, ROWS(a)); if( SCALAR(b)) VAL(b) = IO.u.val; else for( i = 0; i < ROWS(a); i++) B(0,i) = i + IO.u.val; } else m_check(b,"indexing"); if( c == NULL) { c = m_create( 1, COLS(a)); if( SCALAR(c)) VAL(c) = IO.u.val; else for( i = 0; i < COLS(a); i++) C(0,i) = i + IO.u.val; } else m_check(c,"indexing"); for( ib = 0; ib < ROWS(b); ib++) for( jb = 0; jb < COLS(b); jb++) for( ic = 0; ic < ROWS(c); ic++) for( jc = 0; jc < COLS(c); jc++) { ia = (SCALAR(b) ? (int) VAL(b) : (int) B(ib,jb)) - IO.u.val; ja = (SCALAR(c) ? (int) VAL(c) : (int) C(ic,jc)) - IO.u.val; if( ia >= 0 && ia < ROWS(a) && ja >= 0 && ja < COLS(a)) if( SCALAR(a)) VAL(a) += d; else A(ia,ja) += d; else error("index out of range",""); } push(a); } void m_colon( void) { Symbolptr a, b, r; int i, n; double x, inc; b = pop(); a = pop(); n_check(a,":"); n_check(b,":"); n = 1 + abs((int) VAL(b) - (int) VAL(a)); inc = (int) VAL(a) <= (int) VAL(b) ? 1.0 : -1.0; r = m_create( 1, n); if( SCALAR(r)) VAL(r) = (int) VAL(a); else for( i = 0, x = (int) VAL(a); i < n; i++, x += inc) R(0,i) = x; push(r); } void m_index( void) /* a(b) or a(b,c) */ { int narg = (int) *pc++; switch( narg) { case 0: break; case 1: m_index1(); break; case 2: m_index2(); break; default: error("too many indices for variable",""); break; } } void m_index1( void) /* a(b) */ { Symbolptr a, b, c, r; int i, j, k; m_check( a = pop(),"indexing"); b = pop(); if( b == NULL) { /* a(:) ==>> unravel to single row */ r = copy( NULL, a); COLS(r) = COLS(r) * ROWS(r); ROWS(r) = 1; push(r); return; } if( ROWS(a) != 1 && COLS(a) != 1) { /* matrix */ push(b); push(NULL); push(a); m_index2(); return; } m_check( b, "indexing"); if( SCALAR(b)) { r = m_create( 1, 1); k = ((int) VAL(b)) - IO.u.val; if( ROWS(a) == 1 && k >= 0 && k < COLS(a)) VAL(r) = SCALAR(a) ? VAL(a) : A(0,k); else if( COLS(a) == 1 && k >= 0 && k < ROWS(a)) VAL(r) = SCALAR(a) ? VAL(a) : A(k,0); else error("index out of range",""); } else { c = b; if( (ROWS(b) == 1 && ROWS(a) != 1) || (COLS(b) == 1 && COLS(a) != 1)) { c = m_create( COLS(b), ROWS(b)); mattrans( MAT(c), MAT(b), ROWS(b), COLS(b)); } r = m_create( ROWS(c), COLS(c)); for( i = 0; i < ROWS(c); ++i) for( j = 0; j < COLS(c); ++j) { k = ((int) C(i,j)) - IO.u.val; if( ROWS(a) == 1 && k >= 0 && k < COLS(a)) R(i,j) = SCALAR(a) ? VAL(a) : A(0,k); else if( COLS(a) == 1 && k >= 0 && k < ROWS(a)) R(i,j) = SCALAR(a) ? VAL(a) : A(k,0); else error("index out of range",""); } } push(r); } void m_index2( void) /* a(b,c) */ { Symbolptr a, b, c, r; int ia, ja, ib, jb, ic, jc, i, j, rows, cols; a = pop(); c = pop(); b = pop(); m_check(a,"indexing"); if( b == NULL) { b = m_create( 1, ROWS(a)); if( SCALAR(b)) VAL(b) = IO.u.val; else for( i = 0; i < ROWS(a); i++) B(0,i) = i + IO.u.val; } else m_check(b,"indexing"); if( c == NULL) { c = m_create( 1, COLS(a)); if( SCALAR(c)) VAL(c) = IO.u.val; else for( i = 0; i < COLS(a); i++) C(0,i) = i + IO.u.val; } else m_check(c,"indexing"); rows = ROWS(b) * COLS(b); cols = ROWS(c) * COLS(c); r = m_create( rows, cols); for( i = 0, ib = 0; ib < ROWS(b); ib++) for( jb = 0; jb < COLS(b); jb++, i++) for( j = 0, ic = 0; ic < ROWS(c); ic++) for( jc = 0; jc < COLS(c); jc++, j++) { ia = (SCALAR(b) ? (int) VAL(b) : (int) B(ib,jb)) - IO.u.val; ja = (SCALAR(c) ? (int) VAL(c) : (int) C(ic,jc)) - IO.u.val; if( ia >= 0 && ia < ROWS(a) && ja >= 0 && ja < COLS(a)) if( SCALAR(r)) VAL(r) = SCALAR(a) ? VAL(a) : A(ia,ja); else R(i,j) = SCALAR(a) ? VAL(a) : A(ia,ja); else error("index out of range",""); } push(r); } void m_assign2( void) /* a(b,c) = d */ { Symbolptr a, b, c, d; int ia, ja, ib, jb, ic, jc, i, j, rows, cols; a = pop(); d = pop(); c = pop(); b = pop(); m_check(a,"indexing"); m_check(d,"indexed assignment"); if( b == NULL) { b = m_create( 1, ROWS(a)); if( SCALAR(b)) VAL(b) = IO.u.val; else for( i = 0; i < ROWS(a); i++) B(0,i) = i + IO.u.val; } else m_check(b,"indexing"); if( c == NULL) { c = m_create( 1, COLS(a)); if( SCALAR(c)) VAL(c) = IO.u.val; else for( i = 0; i < COLS(a); i++) C(0,i) = i + IO.u.val; } else m_check(c,"indexing"); rows = ROWS(b) * COLS(b); cols = ROWS(c) * COLS(c); if( !SCALAR(d) && (rows != ROWS(d) || cols != COLS(d)) ) error("dimensions incompatible for indexed assignment",""); if( d == a) /* have to make temporary copy since `a' */ d = copy( NULL, a); /* will be changing during the assignment */ for( i = 0, ib = 0; ib < ROWS(b); ib++) for( jb = 0; jb < COLS(b); jb++, i++) for( j = 0, ic = 0; ic < ROWS(c); ic++) for( jc = 0; jc < COLS(c); jc++, j++) { ia = (SCALAR(b) ? (int) VAL(b) : (int) B(ib,jb)) - IO.u.val; ja = (SCALAR(c) ? (int) VAL(c) : (int) C(ic,jc)) - IO.u.val; if( ia >= 0 && ia < ROWS(a) && ja >= 0 && ja < COLS(a)) if( SCALAR(a)) VAL(a) = SCALAR(d) ? VAL(d) : D(i,j); else A(ia,ja) = SCALAR(d) ? VAL(d) : D(i,j); else error("index out of range",""); } push(a); } void m_assign( void) /* a(b) = c or a(b,c) = d */ { int narg = (int) *pc++; switch( narg) { case 0: m_copy(); break; case 1: m_assign1(); break; case 2: m_assign2(); break; default: error("too many indices for variable",""); break; } } void m_assign1( void) /* a(b) = c */ { Symbolptr a, b, c, d; int i, j, k; a = pop(); m_check(a,"indexing"); if( ROWS(a) != 1 && COLS(a) != 1) { /* matrix */ push(NULL); swap(); push(a); m_assign2(); return; } c = pop(); b = pop(); m_check(b,"indexing"); m_check(c,"indexed assignment"); d = b; if( !SCALAR(c) ) { if( !SCALAR(b) && ( (ROWS(b) == 1 && ROWS(a) != 1) || (COLS(b) == 1 && COLS(a) != 1) ) ) { d = m_create( COLS(b), ROWS(b)); mattrans( MAT(d), MAT(b), ROWS(b), COLS(b)); } if( ROWS(d) != ROWS(c) || COLS(d) != COLS(c) ) error("dimensions incompatible for indexed assignment",""); } if( c == a) /* have to make temporary copy since `a' */ c = copy( NULL, a); /* will be changing during the assignment */ for( i = 0; i < ROWS(d); i ++) for( j = 0; j < COLS(d); j++) { k = (SCALAR(d) ? (int) VAL(d) : (int) D(i,j)) - IO.u.val; if( SCALAR(a) && k == 0) VAL(a) = SCALAR(c) ? VAL(c) : C(i,j); else if( ROWS(a) == 1 && k >= 0 && k < COLS(a)) A(0,k) = SCALAR(c) ? VAL(c) : C(i,j); else if( COLS(a) == 1 && k >= 0 && k < ROWS(a)) A(k,0) = SCALAR(c) ? VAL(c) : C(i,j); else error("index out of range",""); } push(a); } void m_string( void) /* convert matrix to string */ { Symbolptr a = pop(); Symbolptr s; int i, j; char *p; check(a); if( a->type == STRING) { push(a); return; } s = s_create( ROWS(a) * COLS(a)); p = STR(s); if( SCALAR(a)) *p++ = VAL(a); else for( i = 0; i < ROWS(a); i++) for( j = 0; j < COLS(a); j++) *p++ = A(i,j); *p = '\0'; push(s); } void m_matrix( void) /* convert string to matrix */ { Symbolptr s = pop(); Symbolptr a; int i, j; char *p; check(s); if( s->type == NUMBER || s->type == MATRIX) { push(s); return; } a = m_create( ROWS(s), COLS(s)); p = STR(s); if( SCALAR(a)) VAL(a) = *p; else for( i = 0; i < ROWS(a); i++) for( j = 0; j < COLS(a); j++) A(i,j) = *p++; push(a); } void m_dyadic( void) /* v1 op v2 */ { int op = (int) *pc++; Symbolptr v1, v2; v2 = pop(); v1 = pop(); push( dyadic( v1, v2, op) ); } void m_install( void) { push( install( NULL, UNDEF, AUTO, 0, 0)); } void m_uninstall( void) { npop( (int) *pc++); /* clear too ??? */ }