/*
 * m_fn.c: built-in functions
 */

#include "m.h"

#define URAND 1	/* types for m_rand() */
#define IRAND 2
#define ERAND 3
#define GRAND 4

/*
 * functions from <math.h>
 */

void m_erf( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "erf()");
    push( unary( sp[0], ERF) );
}

void m_erfc( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "erfc()");
    push( unary( sp[0], ERFC) );
}

void m_abs( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "abs()");
    push( unary( sp[0], ABS) );
}

void m_acos( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "acos()");
    push( unary( sp[0], ACOS) );
}

void m_asin( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "asin()");
    push( unary( sp[0], ASIN) );
}

void m_atan( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "atan()");
    push( unary( sp[0], ATAN) );
}

void m_ceil( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "ceil()");
    push( unary( sp[0], CEIL) );
}

void m_cos( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "cos()");
    push( unary( sp[0], COS) );
}

void m_cosh( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "cosh()");
    push( unary( sp[0], COSH) );
}

void m_exp( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "exp()");
    push( unary( sp[0], EXP) );
}

void m_floor( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "floor()");
    push( unary( sp[0], FLOOR) );
}

void m_log( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "log()");
    push( unary( sp[0], LOG) );
}

void m_log10( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "log10()");
    push( unary( sp[0], LOG10) );
}

void m_sin( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "sin()");
    push( unary( sp[0], SIN) );
}

void m_sinh( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "sinh()");
    push( unary( sp[0], SINH) );
}


void m_sqrt( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "sqrt()");
    push( unary( sp[0], SQRT) );
}

void m_tan( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "tan()");
    push( unary( sp[0], TAN) );
}

void m_tanh( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "tanh()");
    push( unary( sp[0], TANH) );
}

void m_atan2( void)
{
    int narg = (int) *pc++;
    npop(narg);
    if( narg != 2)
	error("atan2() requires two arguments","");
    push( dyadic( sp[0], sp[1], ATAN2) );
}

/*
 * other functions
 */

void m_round( void)
{
    int narg = (int) *pc++;
    npop(narg);
    if( narg != 2)
	error("round() requires two arguments","");
    push( dyadic( sp[0], sp[1], ROUND) );
}

void m_trunc( void)
{
    int narg = (int) *pc++;
    npop(narg);
    if( narg != 2)
	error("trunc() requires two arguments","");
    push( dyadic( sp[0], sp[1], TRUNC) );
}

void m_round2( void)
{
    int narg = (int) *pc++;
    npop(narg);
    if( narg != 2)
	error("round2() requires two arguments","");
    push( dyadic( sp[0], sp[1], ROUND2) );
}

void m_trunc2( void)
{
    int narg = (int) *pc++;
    npop(narg);
    if( narg != 2)
	error("trunc2() requires two arguments","");
    push( dyadic( sp[0], sp[1], TRUNC2) );
}

void m_null( void)
{
    Symbolptr r;
    int i, narg = (int) *pc++;
    npop( narg);

    r = m_create( 1, narg);

    if( SCALAR(r))
	VAL(r) = sp[0] ? 0.0 : 1.0;
    else for( i = 0; i < narg; ++i)
	R(0,i) = sp[i] ? 0.0 : 1.0;

    push(r);
}

void m_undef( void)
{
    Symbolptr r;
    int i, narg = (int) *pc++;
    npop( narg);

    r = m_create( 1, narg);

    if( SCALAR(r))
	VAL(r) = sp[0] ? ((sp[0]->type == UNDEF) ? 1.0 : 0.0) : 1.0;
    else for( i = 0; i < narg; ++i)
	R(0,i) = sp[i] ? ((sp[i]->type == UNDEF) ? 1.0 : 0.0) : 1.0;

    push(r);
}

void m_format( void)
{
    Symbolptr oldformat;
    int narg = (int) *pc++;
    npop(narg);

    if( narg > 1)
	error("format() takes at most one argument","");

    if( narg == 0) {
	push( &format);
	return;
    }

    s_check(sp[0],"format()");

    oldformat = s_create( format.cols);
    strcpy( STR(oldformat), format.u.str);

    format.cols = 0;
    if( (format.u.str = strdup(STR(sp[0]))) == NULL)
	error("strdup() failed in format()","");
    format.cols = strlen(STR(sp[0]));

    push( oldformat);
}

void m_IO( void)
{
    Symbolptr oldIO;
    int narg = (int) *pc++;
    npop(narg);

    if( narg > 1)
	error("IO() takes at most 1 argument","");

    if( narg == 0) {
	push(&IO);
	return;
    }

    n_check( sp[0], "IO()");

    oldIO = m_create( 1, 1);
    VAL(oldIO) = IO.u.val;
    IO.u.val = (int) VAL(sp[0]);

    push( oldIO);
}

void m_cleanup( void)
{
    Symbolptr oldCLEAN;
    int narg = (int) *pc++;
    npop(narg);

    if( narg > 1)
	error("cleanup() takes at most one argument","");

    if( narg == 0) {
	push(&CLEAN);
	return;
    }

    n_check( sp[0], "cleanup()");

    oldCLEAN = m_create( 1, 1);
    VAL(oldCLEAN) = CLEAN.u.val;
    CLEAN.u.val = (int) VAL(sp[0]);

    push( oldCLEAN);
}

void m_eps( void)
{
    Symbolptr oldeps;
    int narg = (int) *pc++;
    npop(narg);

    if( narg > 1)
	error("eps() takes at most one argument","");

    if( narg == 0) {
	push(&eps);
	return;
    }

    n_check( sp[0], "eps()");

    oldeps = m_create( 1, 1);
    VAL(oldeps) = eps.u.val;
    eps.u.val = VAL(sp[0]);

    push( oldeps);
}

void m_tol( void)
{
    Symbolptr oldtol;
    int narg = (int) *pc++;
    npop(narg);

    if( narg > 1)
	error("prtol() takes at most one argument","");

    if( narg == 0) {
	push( &tol);
	return;
    }

    n_check( sp[0], "prtol()");

    oldtol = m_create( 1, 1);
    VAL(oldtol) = tol.u.val;
    tol.u.val = fabs(VAL(sp[0]));

    push( oldtol);
}

void m_num( void)
{
    Symbolptr oldnum;
    int narg = (int) *pc++;
    npop(narg);

    if( narg > 1)
	error("prnum() takes at most one argument","");

    if( narg == 0) {
	push( &num);
	return;
    }

    n_check( sp[0], "prnum()");

    oldnum = m_create( 1, 1);
    VAL(oldnum) = num.u.val;
    num.u.val = (int) fabs(VAL(sp[0]));

    push( oldnum);
}

void m_sizeof( void)
{
    Symbolptr r;
    int narg = (int) *pc++;
    npop(narg);

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

    check(sp[0]);

    r = m_create( 1, 2);

    R(0,0) = ROWS(sp[0]);
    R(0,1) = COLS(sp[0]);

    push(r);
}

void m_strpaste( void)
{
    Symbolptr r;
    int i, len = 0, narg = (int) *pc++;
    npop(narg);

    if( narg < 2)
	error("strpaste() takes two or more arguments","");

    for( i = 0;  i < narg;  ++i) {
	if( sp[i]) {
	    s_check(sp[i], "strpaste()");
	    len += COLS(sp[i]);
	} else
	    ++len;
    }

    r = s_create( len);

    if( sp[0])
        strcpy( STR(r), STR(sp[0]) );
    else
	strcpy( STR(r), "\n");

    for( i = 1;  i < narg;  ++i)
	if( sp[i])
	    strcat( STR(r), STR( sp[i]) );
	else
	    strcat( STR(r), "\n");

    push(r);
}

void m_rows( void)
{
    Symbolptr r;
    int narg = (int) *pc++;
    npop(narg);

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

    check(sp[0]);

    r = m_create( 1, 1);

    VAL(r) = ROWS(sp[0]);

    push(r);
}

void m_cols( void)
{
    Symbolptr r;
    int narg = (int) *pc++;
    npop(narg);

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

    check(sp[0]);

    r = m_create( 1, 1);

    VAL(r) = COLS(sp[0]);

    push(r);
}

void m_norm( void)
{
    Symbolptr r;
    int narg = (int) *pc++;
    npop(narg);

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

    m_check(sp[0],"norm()");

    r = m_create( 1, 1);

    if( SCALAR(sp[0]))
	VAL(r) = fabs(VAL(sp[0]));
    else
	VAL(r) = norm( MAT(sp[0]), ROWS(sp[0]), COLS(sp[0]) );

    push(r);
}

void m_prod( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "prod()");
    push( reduce( sp[0], PROD) );
}

void m_sum( void)
{
    int narg = (int) *pc++;
    unary_check( narg, "sum()");
    push( reduce( sp[0], SUM) );
}

void m_max( void)
{
/*
    int narg = (int) *pc++;
    unary_check( narg, "max()");
    push( reduce( sp[0], MAX) );
*/
    int narg = (int) *pc++;
    npop( narg);

    if( narg < 1 || narg > 2)
	error("max() takes one or two arguments","");

    if( narg == 1) {
	m_check( sp[0], "max()");
	push( reduce( sp[0], MAX) );
    } else
	push( dyadic( sp[0], sp[1], MAX) );
}

void m_min( void)
{
    int narg = (int) *pc++;
    npop( narg);

    if( narg < 1 || narg > 2)
	error("min() takes one or two arguments","");

    if( narg == 1) {
	m_check( sp[0], "min()");
	push( reduce( sp[0], MIN) );
    } else
	push( dyadic( sp[0], sp[1], MIN) );
}

void m_diag( void)
{
    Symbolptr a;
    int narg = (int) *pc++;
    npop(narg);

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

    m_check( a = sp[0], "diag()");

    push( diag(a) );
}

void m_reshape( void)
{
    Symbolptr a, v;
    int m, n, i, j, k, p;
    int narg = (int) *pc++;
    npop(narg);

    if( narg != 3)
	error("reshape() takes 3 arguments","");

    m_check(a=sp[0],"reshape()");
    n_check(sp[1],"reshape()");
    n_check(sp[2],"reshape()");

    m = (int) VAL(sp[1]);
    n = (int) VAL(sp[2]);

    if( m == ROWS(a) && n == COLS(a)) {
	push(a);
	return;
    }

    if( m < 1 || n < 1)
	error("reshape() arguments out of range","");

    v = m_create( m, n);

    if( SCALAR(v))
	VAL(v) = A(0,0);
    else if( SCALAR(a))
	for( i = 0; i < ROWS(v); i++)
	    for( j = 0; j < COLS(v); j++)
		V(i,j) = VAL(a);
    else for( i = k = p = 0; i < ROWS(v); i++)
	    for( j = 0; j < COLS(v); j++) {
		V(i,j) = A(k,p);
		p = (p+1) % COLS(a);
		if( p == 0)
		    k = (k+1) % ROWS(a);
	    }

    push(v);
}

void m_eye( void)
{
    int narg = (int) *pc++;
    Symbolptr r;

    arg_check( narg, "eye()");

    r = m_create( (int) VAL(sp[0]), (int) VAL(sp[1]) );

    if( SCALAR(r) )
	VAL(r) = 1.0;
    else
	eye( MAT(r), ROWS(r), COLS(r));

    push(r);
}

void m_ones( void)
{
    int narg = (int) *pc++;
    Symbolptr r;

    arg_check( narg, "ones()");

    r = m_create( (int) VAL(sp[0]), (int) VAL(sp[1]) );

    if( SCALAR(r))
	VAL(r) = 1.0;
    else
	ones( MAT(r), ROWS(r), COLS(r));

    push(r);
}

void m_rand( void)
{
    int narg = (int) *pc++;
    Symbolptr r;
    int i, j, type, n;

    npop(narg);

    if( narg < 2)
	error( "rand( m, n, ...) takes at least two arguments","");

    n_check(sp[0],"rand()");
    n_check(sp[1],"rand()");

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

    type = URAND;

    if( narg > 2) {
	s_check( sp[2], "rand()");
	switch( *STR(sp[2]) ) {
	    case 'u': type = URAND; break;
	    case 'i': type = IRAND; break;
	    case 'e': type = ERAND; break;
	    case 'g':
	    case 'n': type = GRAND; break;
	    default:
		error("unrecognized distribution type argument to rand()","");
		break;
	}
    }

    if( type == IRAND) {
	if( narg != 4)
	    error("rand( m, n, \"int\", i) requires fourth argument","");
	n_check( sp[3], "rand()");
	n = (int) VAL(sp[3]);
    }

    r = m_create( (int) VAL(sp[0]), (int) VAL(sp[1]) );

    if( SCALAR(r)) switch( type) {

	case URAND: VAL(r) = urand(); break;
	case IRAND: VAL(r) = irand(n); break;
	case ERAND: VAL(r) = erand(); break;
	case GRAND: VAL(r) = grand(); break;
	default: error("rand() type is not valid",""); break;

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

	case URAND: R(i,j) = urand(); break;
	case IRAND: R(i,j) = irand(n); break;
	case ERAND: R(i,j) = erand(); break;
	case GRAND: R(i,j) = grand(); break;
	default: error("rand() type is not valid",""); break;
    }

    push(r);
}

void m_zeros( void)
{
    int narg = (int) *pc++;

    arg_check( narg, "zeros()");

    push( m_create( (int) VAL(sp[0]), (int) VAL(sp[1])) );
}

void m_any( void)
{
    int narg = (int) *pc++;
    Symbolptr a, r;
    int i, j, more = 1;

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

    m_check( a = pop(), "any()");

    r = m_create( 1, 1);

    if( SCALAR(a)) {
	if( VAL(a))
	    VAL(r) = 1.0;
    } else for( i = 0; i < ROWS(a) && more; ++i)
	     for( j = 0; j < COLS(a) && more; ++j)
		if( A(i,j)) {
		    VAL(r) = 1.0;
		    more = 0;
		}

    push(r);
}

void m_all( void)
{
    int narg = (int) *pc++;
    Symbolptr a, r;
    int i, j, more = 1;

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

    m_check( a = pop(), "all()");

    r = m_create( 1, 1);
    VAL(r) = 1.0;

    if( SCALAR(a)) {
	if( VAL(a) == 0.0)
	    VAL(r) = 0.0;
    } else for( i = 0; i < ROWS(a) && more; ++i)
	    for( j = 0; j < COLS(a) && more; ++j)
		if( A(i,j) == 0.0) {
		    VAL(r) = 0.0;
		    more = 0;
		}

    push(r);
}

void m_find( void)
{
    int narg = (int) *pc++;
    Symbolptr a, r;
    int index, i, j, count = 0;

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

    m_check( a = pop(), "find()");

    if( SCALAR(a)) {
	if( VAL(a) != 0.0)
	    ++count;
    } else for( i = 0; i < ROWS(a); ++i)
	    for( j = 0; j < COLS(a); ++j)
		if( A(i,j) != 0.0)
		    ++count;

    r = m_create( 1, count);
    push(r);
    if( count == 0)
	return;

    index = IO.u.val;

    if( SCALAR(a)) {

	if( VAL(a) != 0.0)
	    if( SCALAR(r)) VAL(r) = index;
	    else R(0,0) = index;

    } else for( count = i = 0; i < ROWS(a); ++i)
	    for( j = 0; j < COLS(a); ++j, ++index)
		if( A(i,j) != 0.0)
		    if( SCALAR(r)) { VAL(r) = index; return; }
		    else { R(0,count) = index; ++count; }

}

void m_svd( void)	/* svd( a, u, s, v) */
{
    Symbolptr a, u, s, S, v, rank;
    int i, m, n;
    int narg = (int) *pc++;
    npop(narg);

    if( narg < 1 || narg > 4)
	error("svd( a, u, s, v) takes one to four arguments","");

    m_check( a = sp[0], "svd()");

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

    if( narg > 1) u = sp[1]; else u = NULL;

    if( narg > 2) S = sp[2]; else S = NULL;

    if( narg > 3) v = sp[3]; else v = NULL;

    if( u) {
	u = m_new( u, m, m);
	if( SCALAR(u))
	    VAL(u) = 1.0;
	else
	    eye( MAT(u), m, m);
    }

    if( S)
	S = m_new( S, m, n);

    v = m_new( v, m > n ? m : n, n);

    if( SCALAR(v))
	VAL(v) = VAL(a);
    else {
	matcopy( MAT(v), MAT(a), m, n);
	s = m_create( 1, n == 1 ? 2 : n);	/* make sure s is MATRIX */
    }

    rank = m_create( 1, 1);
    push(rank);

    if( SCALAR(v)) {
	if( S)
	    VAL(S) = fabs(VAL(v));
	VAL(rank) = VAL(v) == 0.0 ? 0.0 : 1.0;
	VAL(v) = VAL(v) < 0.0 ? -1.0 : 1.0;
	return;
    }

    /*
     * note: u will not be referenced in svd() if m == 1
     */
    VAL(rank) = svd( (u && m > 1) ? MAT(u) : NULL, MAT(s), MAT(v), m, n);

    ROWS(v) = n;	/* resize v */

    if( ROWS(v) == 1 && COLS(v) == 1) {
	double r = MAT(v)[0][0];
	v->type = NUMBER;
	matfree( MAT(v));
	VAL(v) = r;
    }

    if( S) {
	m = m < n ? m : n;		/* set diag(S) = s */
	for( i = 0; i < m; ++i)
	    MAT(S)[i][i] = MAT(s)[0][i];
    }
}