#include "m.h"

void m_print()
{
    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()	/* 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()
{
    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()
{
    Symbolptr a = pop();
    int op = (int) *pc++;

    m_check(a,"unary operator");

    push( unary( a, op) );
}

void m_reduce()
{
    Symbolptr a = pop();
    int op = (int) *pc++;

    m_check(a,"reduce operator");

    push( reduce( a, op) );
}

void m_copy()
{
    Symbolptr a, b;

    a = pop();
    b = pop();
    a = copy( a, b);	/* a <<-- b */
    push(a);
}

void m_pp()	/* 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()	/* ++a(b) or ++a(b,c) where increment is d */
{
    Symbolptr a, b, c, r;
    int d, narg;

    a = pop();
    d = (int) *pc++;
    narg = (int) *pc++;

    m_check(a,"++ or --");
    c_check(a,"++ or --");

    switch( narg) {

	case 0: if( SCALAR(a) )
		  VAL(a) += d;
		else
		  matsdyadic( MAT(a), MAT(a), (double)d, ROWS(a), COLS(a), '+');
		push(a);
		return;

	case 1: b=pop(); push(b); push(b); push(a); m_index1();
		break;

	case 2: c=pop(); b=pop(); push(b); push(c); push(b); push(c);
		m_index2();
		break;

	default: error("m_pp2() called with more than two args","");
    }
}

void m_colon()
{
    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()	/* 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()	/* 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()	/* 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()	/* 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()	/* 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()	/* 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( ROWS(c) == 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));
	}
*/
	if( ROWS(b) != ROWS(c) || COLS(b) != 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(b); i ++)
	for( j = 0; j < COLS(b); j++) {

	    k = (SCALAR(b) ? (int) VAL(b) : (int) B(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()	/* 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()	/* 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()	/* v1 op v2 */
{
    int op = (int) *pc++;
    Symbolptr v1, v2;

    v2 = pop();
    v1 = pop();

    push( dyadic( v1, v2, op) );
}

void m_install()
{
    push( install( NULL, UNDEF, AUTO, 0, 0));
}

void m_uninstall()
{
    npop( (int) *pc++);

    /* clear too ??? */
}