#include "m.h"

#define FORMAT "%15.7g"	/* default format */

#ifndef min
# define min(a,b) ( (a) < (b) ? (a) : (b) )
# define max(a,b) ( (a) > (b) ? (a) : (b) )
#endif

extern double erf( double), erfc( double);	/* not in ANSI C math.h */

Symbol eps, tol, num, format, IO, CLEAN, rowend, LTsym, GTsym, GTGTsym;

/*
 * name for printing temporary variables which have no name
 */
char *NONAME = "--noname--";

void m_init( void)
{
    fin = stdin;

    eps.name = "eps";
    eps.type = NUMBER;
    eps.class = DEFN;
    eps.rows = eps.cols = 1;
    eps.u.val = geteps();

    tol.name = "prtol";
    tol.type = NUMBER;
    tol.class = DEFN;
    tol.rows = tol.cols = 1;
    tol.u.val = 100.0 * eps.u.val;

    num.name = "prnum";
    num.type = NUMBER;
    num.class = DEFN;
    num.rows = num.cols = 1;
    num.u.val = 5;

    format.name = "format";
    format.type = STRING;
    format.class = DEFN;
    format.rows = 1;
    format.cols = strlen(FORMAT);
    if( (format.u.str = strdup(FORMAT)) == NULL)
	error("strdup() failed in m_init()","");

    IO.name = "IO";
    IO.type = NUMBER;
    IO.class = DEFN;
    IO.rows = IO.cols = 1;
    IO.u.val = 1.0;

    CLEAN.name = "cleanup";
    CLEAN.type = NUMBER;
    CLEAN.class = DEFN;
    CLEAN.rows = CLEAN.cols = 1;
    CLEAN.u.val = 10.0;

    rowend.name = "rowend";
    rowend.type = NUMBER;
    rowend.class = DEFN;
    rowend.rows = rowend.cols = 1;
    rowend.u.val = 0.0;

    LTsym.name = "<";
    LTsym.type = NUMBER;
    LTsym.class = DEFN;
    LTsym.rows = LTsym.cols = 1;
    LTsym.u.val = 0.0;

    GTsym.name = ">";
    GTsym.type = NUMBER;
    GTsym.class = DEFN;
    GTsym.rows = GTsym.cols = 1;
    GTsym.u.val = 0.0;

    GTGTsym.name = ">>";
    GTGTsym.type = NUMBER;
    GTGTsym.class = DEFN;
    GTGTsym.rows = GTGTsym.cols = 1;
    GTGTsym.u.val = 0.0;
}

char *str_type( int type)
{
    switch( type) {
	case UNDEF:	return "UNDEF";
	case NUMBER:	return "NUMBER";
	case MATRIX:	return "matrix";
	case STRING:	return "String";
	case FN:	return "FN";
	case CMD:	return "CMD";
	case PROC:	return "PROC";
	case FUNCTION:	return "FUNCTION";
	case PROCEDURE:	return "PROCEDURE";

	case VECTOR:	return "vector";	/* for function args */
	case INT:	return "int";
	case DOUBLE:	return "double";
	default:	return "unknown";
    }
}

char *str_class( int class)
{
    switch( class) {
	case NEW:	return "NEW";
	case EXTERN:	return "extern";
	case AUTO:	return "auto";
	case STATIC:	return "static";
	case DEFN:	return "DEFN";
	case LOOP:	return "LOOP";
	case CONST:	return "CONST";

	case MATRIX:	return "matrix";	/* FUNCTION return type */
	case STRING:	return "String";
	case INT:	return "int";
	case DOUBLE:	return "double";
	case VECTOR:	return "vector";
	case VOID:	return "void";
	default:	return "unknown";
    }
}

void check( Symbolptr a)
{
    if( !a)
	error("null argument","");
    else if( a->type == UNDEF)
	error("undefined variable ", NAME(a) );
}

void c_check( Symbolptr a, char *name)
{
    if( a && (a->class == CONST || a->class == DEFN || a->class == LOOP) )
	error("attempt to alter constant by ", name);
}

void m_check( Symbolptr a, char *name)
{
    check(a);
    if( a->type == STRING)
	error("numeric argument required for ",name);
}

void n_check( Symbolptr a, char *name)
{
    check(a);
    if( !SCALAR(a))
	error("scalar argument required for ",name);
}

void s_check( Symbolptr a, char *name)
{
    check(a);
    if( a->type != STRING)
	error("string argument required for ",name);
}

void print( Symbolptr s)
{
    /* make sure the output starts on a new line */

    if( diary)
	putc( '\n', fdiary);

    if( echo && include_stack_ptr > 0)
	putchar( '\n');

    /* print name and value */

    if( s->name)
	dprint( "%s = ", s->name);

    if( s->type == STRING)

	dprint( "\"%s\"\n", STR(s) ? STR(s) : NONAME );

    else if( s->type == MATRIX && MAT(s)) {

	matprint( stdout, format.u.str, MAT(s), ROWS(s), COLS(s),
			tol.u.val, (int) IO.u.val, (int) num.u.val);
	if( diary)
	  matprint( fdiary, format.u.str, MAT(s), ROWS(s), COLS(s),
			tol.u.val, (int) IO.u.val, (int) num.u.val);

    } else if( s->type == NUMBER) {

	dprint( format.u.str, VAL(s));
	dprint("\n");
    }
}

void arg_check( int narg, char *name)	/* used by zeros(), ones(), eye() () */
{
    npop(narg);

    if( narg != 2)
	error( name, " takes two arguments");

    n_check(sp[0],name);
    n_check(sp[1],name);

    if( (int) VAL(sp[0]) < 1 || (int) VAL(sp[1]) < 1)
	error("arguments out of range for ", name);
}

void unary_check( int narg, char *name)
{
    npop( narg);

    if( narg != 1)
	error( name, " takes one argument");

    m_check( sp[0], name);
}

/*
 * m_create: create a CONST unnamed NUMBER or MATRIX entry in symbol table
 */
Symbolptr m_create( int m, int n)
{
    Symbolptr a;
    int type;

    type = (m == 1 && n == 1) ? NUMBER : MATRIX;

    a = install( NULL, type, CONST, m, n);

    if( type == NUMBER)
	VAL(a) = 0.0;
    else if( (MAT(a) = matcreate(m,n)) == NULL && m > 0 && n > 0 )
	error("matcreate() failed in m_create()","");

    return a;
}

/*
 * s_create: create a CONST unnamed STRING entry in symbol table
 */
Symbolptr s_create( int n)
{
    Symbolptr s;

    s = install( NULL, STRING, CONST, 1, n);

    if( (STR(s) = (char *) malloc(n+1)) == NULL)
	error("malloc() failed in s_create()","");

    return s;
}

Symbolptr m_new( Symbolptr a, int m, int n)
{
    int class;		/* coerce a to be m-by-n initialized to zeros */

    if( !a)
	return m_create( m, n);

    c_check( a, "assignment");

    if( a->type == MATRIX && ROWS(a) == m && COLS(a) == n) {

	zeros( MAT(a), m, n);
	return a;

    } else if( a->type == NUMBER && m == 1 && n == 1) {

	VAL(a) = 0.0;
	return a;
    }

    class = a->class;
    clear(a);

    a->type = (m == 1 && n == 1) ? NUMBER : MATRIX;

    if( SCALAR(a))
	VAL(a) = 0.0;
    else if( (MAT(a) = matcreate(m,n)) == NULL && m > 0 && n > 0)
	error("matcreate() failed in m_new()","");

    ROWS(a) = m;
    COLS(a) = n;

    if( class == NEW)
	a->class = EXTERN;
    else
	a->class = class;

    return a;
}

double sdyadic( double x, double y, int op)	/* return x op y */
{
    double r;

    switch(op) {

	case '+':	r = x +  y; break;
	case '-':	r = x -  y; break;
	case '/':	r = x /  y; break;
	case '*':	r = x *  y; break;
	case OR_OP: 	r = x || y; break;
	case AND_OP:	r = x && y; break;
	case EQ_OP:	r = x == y; break;
	case NE_OP:	r = x != y; break;
	case '<':	r = x <  y; break;
	case LE_OP:	r = x <= y; break;
	case '>':	r = x >  y; break;
	case GE_OP:	r = x >= y; break;
	case '%':	r = fmod( x, y); break;
	case '^':	r = pow( x, y); break;
	case ATAN2:	r = atan2( x, y); break;
	case TRUNC:	r = trunc( x, y); break;
	case ROUND:	r = round( x, y); break;
	case TRUNC2:	r = trunc2( x, y); break;
	case ROUND2:	r = round2( x, y); break;
	case MIN:	r = min( x, y); break;
	case MAX:	r = max( x, y); break;
	default:	error("unknown dyadic operator","");
    }
    return r;
}

void matsdyadic( double *x[], double *y[], double s, int m, int n, int op)
{
    int i, j;				/* x = y op s */

    for( i = 0; i < m; i++)  for( j = 0; j < n; j++)  switch(op) {

	case '+':	x[i][j] = y[i][j] +  s; break;
	case '-':	x[i][j] = y[i][j] -  s; break;
	case '/':	x[i][j] = y[i][j] /  s; break;
	case '*':	x[i][j] = y[i][j] *  s; break;
	case OR_OP: 	x[i][j] = y[i][j] || s; break;
	case AND_OP:	x[i][j] = y[i][j] && s; break;
	case EQ_OP: 	x[i][j] = y[i][j] == s; break;
	case NE_OP: 	x[i][j] = y[i][j] != s; break;
	case '<': 	x[i][j] = y[i][j] <  s; break;
	case LE_OP: 	x[i][j] = y[i][j] <= s; break;
	case '>': 	x[i][j] = y[i][j] >  s; break;
	case GE_OP: 	x[i][j] = y[i][j] >= s; break;
	case '%':	x[i][j] = fmod( y[i][j], s); break;
	case '^':	x[i][j] = pow( y[i][j], s); break;
	case ATAN2:	x[i][j] = atan2( y[i][j], s); break;
	case TRUNC:	x[i][j] = trunc( y[i][j], s); break;
	case ROUND:	x[i][j] = round( y[i][j], s); break;
	case TRUNC2:	x[i][j] = trunc2( y[i][j], s); break;
	case ROUND2:	x[i][j] = round2( y[i][j], s); break;
	case MIN:	x[i][j] = min( y[i][j], s); break;
	case MAX:	x[i][j] = max( y[i][j], s); break;
	default:	error("unknown dyadic operator",""); break;
    }
}

void smatdyadic( double *x[], double s, double *y[], int m, int n, int op)
{
    int i, j;				/* x = s op y */

    for( i = 0; i < m; i++)  for( j = 0; j < n; j++)  switch(op) {

	case '+':	x[i][j] = s +  y[i][j]; break;
	case '-':	x[i][j] = s -  y[i][j]; break;
	case '/':	x[i][j] = s /  y[i][j]; break;
	case '*':	x[i][j] = s *  y[i][j]; break;
	case OR_OP: 	x[i][j] = s || y[i][j]; break;
	case AND_OP:	x[i][j] = s && y[i][j]; break;
	case EQ_OP: 	x[i][j] = s == y[i][j]; break;
	case NE_OP: 	x[i][j] = s != y[i][j]; break;
	case '<': 	x[i][j] = s <  y[i][j]; break;
	case LE_OP: 	x[i][j] = s <= y[i][j]; break;
	case '>': 	x[i][j] = s >  y[i][j]; break;
	case GE_OP: 	x[i][j] = s >= y[i][j]; break;
	case '%':	x[i][j] = fmod( s, y[i][j]); break;
	case '^':	x[i][j] = pow( s, y[i][j]); break;
	case ATAN2:	x[i][j] = atan2( s, y[i][j]); break;
	case TRUNC:	x[i][j] = trunc( s, y[i][j]); break;
	case ROUND:	x[i][j] = round( s, y[i][j]); break;
	case TRUNC2:	x[i][j] = trunc2( s, y[i][j]); break;
	case ROUND2:	x[i][j] = round2( s, y[i][j]); break;
	case MIN:	x[i][j] = min( s, y[i][j]); break;
	case MAX:	x[i][j] = max( s, y[i][j]); break;
	default:	error("unknown dyadic operator",""); break;
    }
}

void matdyadic( double *x[], double *y[], double *z[], int m, int n, int op)
{
    int i, j;				/* x = y op z */

    for( i = 0; i < m; i++)  for( j = 0; j < n; j++)  switch(op) {

	case '+':	x[i][j] = y[i][j] +  z[i][j]; break;
	case '-':	x[i][j] = y[i][j] -  z[i][j]; break;
	case '/':	x[i][j] = y[i][j] /  z[i][j]; break;
	case '*':	x[i][j] = y[i][j] *  z[i][j]; break;
	case OR_OP: 	x[i][j] = y[i][j] || z[i][j]; break;
	case AND_OP:	x[i][j] = y[i][j] && z[i][j]; break;
	case EQ_OP: 	x[i][j] = y[i][j] == z[i][j]; break;
	case NE_OP: 	x[i][j] = y[i][j] != z[i][j]; break;
	case '<': 	x[i][j] = y[i][j] <  z[i][j]; break;
	case LE_OP: 	x[i][j] = y[i][j] <= z[i][j]; break;
	case '>': 	x[i][j] = y[i][j] >  z[i][j]; break;
	case GE_OP: 	x[i][j] = y[i][j] >= z[i][j]; break;
	case '%':	x[i][j] = fmod( y[i][j], z[i][j]); break;
	case '^':	x[i][j] = pow( y[i][j], z[i][j]); break;
	case ATAN2:	x[i][j] = atan2( y[i][j], z[i][j]); break;
	case TRUNC:	x[i][j] = trunc( y[i][j], z[i][j]); break;
	case ROUND:	x[i][j] = round( y[i][j], z[i][j]); break;
	case TRUNC2:	x[i][j] = trunc2( y[i][j], z[i][j]); break;
	case ROUND2:	x[i][j] = round2( y[i][j], z[i][j]); break;
	case MIN: 	x[i][j] = min( y[i][j], z[i][j]); break;
	case MAX:	x[i][j] = max( y[i][j], z[i][j]); break;
	default:	error("unknown dyadic operator",""); break;
    }
}

Symbolptr unary( Symbolptr a, int op)
{
    Symbolptr v;
    int i, j;

    if( a->class == CONST)
	v = a;
    else
	v = m_create( ROWS(a), COLS(a));

    if( SCALAR(a)) switch(op) {

	case '-':	VAL(v) = - VAL(a); break;
	case '!':	VAL(v) = ! VAL(a); break;
	case INT:	VAL(v) = (int) VAL(a); break;
	case SQRT:	VAL(v) = sqrt( VAL(a)); break;
	case ABS:	VAL(v) = fabs( VAL(a)); break;
	case ERF:	VAL(v) = erf( VAL(a)); break;
	case ERFC:	VAL(v) = erfc( VAL(a)); break;
	case SIN:	VAL(v) = sin( VAL(a)); break;
	case COS:	VAL(v) = cos( VAL(a)); break;
	case TAN:	VAL(v) = tan( VAL(a)); break;
	case ASIN:	VAL(v) = asin( VAL(a)); break;
	case ACOS:	VAL(v) = acos( VAL(a)); break;
	case ATAN:	VAL(v) = atan( VAL(a)); break;
	case SINH:	VAL(v) = sinh( VAL(a)); break;
	case COSH:	VAL(v) = cosh( VAL(a)); break;
	case TANH:	VAL(v) = tanh( VAL(a)); break;
	case EXP:	VAL(v) = exp( VAL(a)); break;
	case LOG:	VAL(v) = log( VAL(a)); break;
	case LOG10:	VAL(v) = log10( VAL(a)); break;
	case CEIL:	VAL(v) = ceil( VAL(a)); break;
	case FLOOR:	VAL(v) = floor( VAL(a)); break;
	default: error("unknown unary operator",""); break;

    } else for( i=0; i<ROWS(a); i++)  for( j=0; j<COLS(a); j++)  switch(op) {

	case '-':	V(i,j) = - A(i,j); break;
	case '!':	V(i,j) = ! A(i,j); break;
	case INT:	V(i,j) = (int) A(i,j); break;
	case SQRT:	V(i,j) = sqrt( A(i,j)); break;
	case ABS:	V(i,j) = fabs( A(i,j)); break;
	case ERF:	V(i,j) = erf( A(i,j)); break;
	case ERFC:	V(i,j) = erfc( A(i,j)); break;
	case SIN:	V(i,j) = sin( A(i,j)); break;
	case COS:	V(i,j) = cos( A(i,j)); break;
	case TAN:	V(i,j) = tan( A(i,j)); break;
	case ASIN:	V(i,j) = asin( A(i,j)); break;
	case ACOS:	V(i,j) = acos( A(i,j)); break;
	case ATAN:	V(i,j) = atan( A(i,j)); break;
	case SINH:	V(i,j) = sinh( A(i,j)); break;
	case COSH:	V(i,j) = cosh( A(i,j)); break;
	case TANH:	V(i,j) = tanh( A(i,j)); break;
	case EXP:	V(i,j) = exp( A(i,j)); break;
	case LOG:	V(i,j) = log( A(i,j)); break;
	case LOG10:	V(i,j) = log10( A(i,j)); break;
	case CEIL:	V(i,j) = ceil( A(i,j)); break;
	case FLOOR:	V(i,j) = floor( A(i,j)); break;
	default: error("unknown unary operator",""); break;
    }

    return v;
}

Symbolptr reduce( Symbolptr a, int op)
{
    Symbolptr v;
    int i, j;
    double r;

    v = m_create( 1, ROWS(a));

    if( SCALAR(a))

	VAL(v) = VAL(a);

    else for( i=0; i < ROWS(a); i++) {

	switch(op) {
	    case SUM: case SUB: r = 0; break;
	    case PROD: case DIV: r = 1; break;
	    case MAX: case MIN: r = A(i,COLS(a)-1); break;
	    default: error("unknown reduce operator",""); break;
	}

 	for( j=COLS(a)-1; j >= 0; j--)  switch(op) {

	    case SUM: r = A(i,j) + r; break;
	    case SUB: r = A(i,j) - r; break;
	    case PROD: r = A(i,j) * r; break;
	    case DIV: r = A(i,j) / r; break;
	    case MAX: r = A(i,j) > r ? A(i,j) : r; break;
	    case MIN: r = A(i,j) < r ? A(i,j) : r; break;
	}

	if( SCALAR(v) )
	    VAL(v) = r;
	else
	    V(0,i) = r;
    }
    return v;
}

Symbolptr dyadic( Symbolptr v1, Symbolptr v2, int op)
{
    Symbolptr v;
    int m, n;

    m_check(v1,"dyadic operator");
    m_check(v2,"dyadic operator");

    if( SCALAR(v1) == SCALAR(v2)) {	/* both scalars or matrices */

	m = ROWS(v1);
	n = COLS(v1);
	if( m != ROWS(v2) || n != COLS(v2))
	    error("incompatible sizes for dyadic operator","");

    } else if( SCALAR(v2) ) {		/* && !SCALAR(v1) */

	m = ROWS(v1);
	n = COLS(v1);

    } else {				/* SCALAR(v1) && !SCALAR(v2) */

	m = ROWS(v2);
	n = COLS(v2);
    }

    v = m_create( m, n);

    if( SCALAR(v) )			/* both scalars */

	VAL(v) = sdyadic( VAL(v1), VAL(v2), op);

    else if( SCALAR(v1) == SCALAR(v2))	/* both matrices */

	matdyadic( MAT(v), MAT(v1), MAT(v2), m, n, op);

    else if( SCALAR(v2))		/* && !SCALAR(v1) */

	matsdyadic( MAT(v), MAT(v1), VAL(v2), m, n, op);

    else				/* SCALAR(v1) && !SCALAR(v2) */

	smatdyadic( MAT(v), VAL(v1), MAT(v2), m, n, op);

    return v;
}

Symbolptr copy( Symbolptr a, Symbolptr b)	/* a <- b */
{
    int class;

    check(b);

    c_check( a, "assignment");

    if( b->class == CONST) {

	if(a) {
	    class = a->class;	/* added these three lines, 2/25/92 */
	    clear(a);
	    a->class = class;
	    if( a->class == NEW)
		a->class = EXTERN;
	    a->type = b->type;
	    ROWS(a) = ROWS(b);
	    COLS(a) = COLS(b);
	    a->u = b->u;
	    b->type = UNDEF;
	    return a;
	} else
	    return b;
    }

    if( a == NULL) /* making a temporary copy, as in VAR++ */
	a = install( NULL, UNDEF, CONST, 0, 0);

    if( ROWS(b) < 1 || COLS(b) < 1) {
	clear(a);
	return a;
    }

    if( a->class == NEW)
	a->class = EXTERN;

    if( a->type == b->type) /* if compatible types and sizes, just copy */
	if( SCALAR(a)) {
	    VAL(a) = VAL(b);
	    return a;
	} else if( a->type == STRING && COLS(a) == COLS(b)) {
	    strcpy( STR(a), STR(b));
	    return a;
	} else if( a->type == MATRIX
			&& ROWS(a) == ROWS(b) && COLS(a) == COLS(b)) {
	    matcopy(MAT(a),MAT(b),ROWS(b),COLS(b));
	    return a;
	}

    class = a->class;
    clear(a);
    a->class = class;

    a->type = b->type;
    ROWS(a) = ROWS(b);
    COLS(a) = COLS(b);

    if( b->type == STRING) {
	if( (STR(a) = strdup(STR(b))) == NULL)
	    error("strdup() failed in m_copy()","");
    } else if( SCALAR(b))
	VAL(a) = VAL(b);
    else if( b->type == MATRIX) {
	if( (MAT(a) = matcreate( ROWS(b), COLS(b))) == NULL)
	    error("matcreate() failed in m_copy()","");
	matcopy(MAT(a),MAT(b),ROWS(b),COLS(b));
    } else {
	dprint( "copy(): b->type == %d\n", b->type);
	error("bad type in copy()","");
    }

    return a;
}

Symbolptr diag( Symbolptr a)
{
    Symbolptr r;
    int i, j, m, n, q;

    m = ROWS(a);
    n = COLS(a);

    if( SCALAR(a))				/* scalar */

	r = copy( NULL, a);

    else if( m > 1 && n > 1) {			/* matrix */

	q = m < n ? m : n;
	r = m_create( 1, q);
	for( i = 0; i < q; ++i)
	    R(0,i) = A(i,i);

    } else {					/* vector */

	q = m > n ? m : n;
	r = m_create( q, q);
	for( q = 0, i = 0; i < m; ++i)
	    for( j = 0; j < n; ++j, ++q)
		R(q,q) = A(i,j);
    }

    return r;
}

/*
double round( x, s)
double x, s;
{
    static double ten[] = { 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10,
			    1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18 };
    static int NMAX = 18;
    double y;
    int n, exp;

    n = (int) s;

    if( n < 1)
	return 0.0;
    else if( n > NMAX)
	return x;

    y = (int) (ten[n-1] * frexp( x, &exp) + 0.5);
    y = ldexp(y/ten[n-1],exp);
    return x < 0 ? -y : y;
}
*/

double round( double x, double s)
{
    char line[30];
    int n;

    n = (int) s;

    if( n < 1)
	return 0.0;
    else if( n > 16)
	return x;

    sprintf( line, "%26.*e", n-1, x);
    return atof( line);
}

double trunc( double x, double s)
{
    char *p, line[30];
    int n;

    n = (int) s;

    if( n < 1)
	return 0.0;
    else if( n > 16)
	return x;

    sprintf( line, "%26.18e", x);
    if( (p = strchr( line, 'e')) != NULL) {
	n = 19 - n;
	while( n-- > 0)
	   *--p = '0';
    }
    return atof( line);
}

double round2( double x, double s)
{
    double y;
    int n, exp;

    n = (int) s;

    if( n < 1)
	return 0.0;
    else if( n > 31)
	return x;

    y = (long) ((1L << n) * frexp( x, &exp) + 0.5);
    y = ldexp(y,exp-n);
    return x < 0 ? -y : y;
}

double trunc2( double x, double s)
{
    double y;
    int n, exp;

    n = (int) s;

    if( n < 1)
	return 0.0;
    else if( n > 31)
	return x;

    y = (long) ((1L << n) * frexp( x, &exp));
    y = ldexp(y,exp-n);
    return x < 0 ? -y : y;
}