/*
    This file is part of the FElt finite element analysis package.
    Copyright (C) 1993-2000 Jason I. Gobat and Darren C. Atkinson

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/

/************************************************************************
 * File:	mathfunc.c						*
 *									*
 * Description:	This file contains the function definitions for the	*
 *		math library intrinsic functions.			*
 ************************************************************************/

# include <math.h>
# include "debug.h"
# include "coerce.h"
# include "execute.h"
# include "mathfunc.h"
# include "arithmetic.h"


/************************************************************************
 * Function:	ceil_func						*
 *									*
 * Description:	Pops and computes the ceiling of the descriptor on the	*
 *		top of the stack and places the result on the stack.	*
 *		The following types are legal:				*
 *									*
 *		ceil (double) -> double (scalar function call)		*
 *		ceil (matrix) -> matrix (matrix element-wise operation)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int ceil_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	CreateData (result, arg, NULL, T_Double);
	*D_Double (result) = ceil (*D_Double (arg));
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		sdata (b, i, j) = ceil (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("ceil", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("ceil ans =\n");
    d_PrintData (result);

    return type_error == F_True;
}


/************************************************************************
 * Function:	cos_func						*
 *									*
 * Description:	Pops and computes the cosine of the descriptor on the	*
 *		top of the stack and places the result on the stack.	*
 *		The following types are legal:				*
 *									*
 *		cos (double) -> double (scalar function call)		*
 *		cos (matrix) -> matrix (matrix element-wise operation)	*
 *									*
 * 		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int cos_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	CreateData (result, arg, NULL, T_Double);
	*D_Double (result) = cos (*D_Double (arg));
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		sdata (b, i, j) = cos (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("cos", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("cos ans =\n");
    d_PrintData (result);

    return type_error == F_True;
}


/************************************************************************
 * Function:	exp_func						*
 *									*
 * Description:	Pops and computes the exponential of the descriptor on	*
 *		the top of the stack and places the result on the	*
 *		stack.  The following types are legal:			*
 *									*
 *		exp (double) -> double (scalar function call)		*
 *		exp (matrix) -> matrix (matrix element-wise operation)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		a double value.						*
 ************************************************************************/

int exp_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	CreateData (result, arg, NULL, T_Double);
	*D_Double (result) = exp (*D_Double (arg));
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		sdata (b, i, j) = exp (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("exp", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("exp ans =\n");
    d_PrintData (result);

    return type_error == F_True;
}


/************************************************************************
 * Function:	fabs_func						*
 *									*
 * Description:	Pops and computes the absolute value of the descriptor	*
 *		on the top of the stack and places the result on the	*
 *		stack.  The following types are legal:			*
 *									*
 *		fabs (double) -> double	(scalar function call)		*
 *		fabs (matrix) -> matrix (matrix element-wise operation)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int fabs_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	CreateData (result, arg, NULL, T_Double);
	*D_Double (result) = fabs (*D_Double (arg));
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		sdata (b, i, j) = fabs (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("fabs", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("fabs ans =\n");
    d_PrintData (result);

    return type_error == F_True;
}


/************************************************************************
 * Function:	floor_func						*
 *									*
 * Description:	Pops and computes the floor of the descriptor on the	*
 *		top of the stack and places the result on the stack.	*
 *		The following types are legal:				*
 *									*
 *		floor (double) -> double (scalar function call)		*
 *		floor (matrix) -> matrix (matrix element-wise op)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int floor_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	CreateData (result, arg, NULL, T_Double);
	*D_Double (result) = floor (*D_Double (arg));
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		sdata (b, i, j) = floor (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("floor", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("floor ans =\n");
    d_PrintData (result);

    return type_error == F_True;
}


/************************************************************************
 * Function:	fmod_func						*
 *									*
 * Description:	Pops and modulos the two descriptors on the top of the	*
 *		stack and places the result on the stack.  The		*
 *		following types are legal:				*
 *									*
 *		fmod (double, double) -> double (scalar function call)	*
 *		fmod (double, matrix) -> matrix (matrix element-wise)	*
 *		fmod (matrix, double) -> matrix (matrix element-wise)	*
 *		fmod (matrix, matrix) -> matrix (matrix element-wise)	*
 *									*
 *		An attempt is first made to coerce the arguments to	*
 *		double values.  The mod operator is used to compute the	*
 *		result, which will result in a misleading message in	*
 *		the event of an error.					*
 ************************************************************************/

int fmod_func (n)
    int n;
{
    return mod_op ( );
}


/************************************************************************
 * Function:	hypot_func						*
 *									*
 * Description:	Pops and computes the euclidean distance of the two	*
 *		descriptors on the top of the stack and places the	*
 *		result on the stack.  The following types are legal:	*
 *									*
 *		hypot (double, double) -> double (scalar function call)	*
 *		hypot (matrix, matrix) -> matrix (matrix element-wise)	*
 *									*
 *		An attempt is first made to coerce the arguments to	*
 *		double values.						*
 ************************************************************************/

int hypot_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    Matrix	c;
    descriptor *arg1;
    descriptor *arg2;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;
    unsigned	i;
    unsigned	j;


    arg2 = pop ( );
    result = top ( );
    temp = *result;
    arg1 = &temp;


    arg1 = CoerceData (deref (arg1), T_Double);
    arg2 = CoerceData (deref (arg2), T_Double);


    type_error = F_False;
    status = 0;


    switch (D_Type (arg1)) {
    case T_Double:
	switch (D_Type (arg2)) {
	case T_Double:
	    CreateData (result, arg1, arg2, T_Double);
	    *D_Double (result) = hypot (*D_Double (arg1), *D_Double (arg2));
	    break;


	case T_Matrix:
	    a = D_Matrix (arg2);
	    CreateData (result, arg1, arg2, T_Matrix, Mrows (a), Mcols (a));
	    b = D_Matrix (arg1);
	    for (i = 1; i <= Mrows (a); i ++)
		for (j = 1; j <= Mcols (a); j ++)
		    sdata (b, i, j) = hypot (*D_Double (arg1), mdata (a, i, j));
	    break;


	default:
	    type_error = F_True;
	    break;
	}
	break;


    case T_Matrix:
	switch (D_Type (arg2)) {
	case T_Double:
	    a = D_Matrix (arg1);
	    CreateData (result, arg1, arg2, T_Matrix, Mrows (a), Mcols (a));
	    b = D_Matrix (result);
	    for (i = 1; i <= Mrows (a); i ++)
		for (j = 1; j <= Mcols (a); j ++)
		    sdata (b, i, j) = hypot (mdata (a, i, j), *D_Double (arg2));
	    break;


	case T_Matrix:
	    a = D_Matrix (arg1);
	    b = D_Matrix (arg2);
	    if (Mrows (a) == Mrows (b) && Mcols (a) == Mrows (a)) {
		CreateData (result, arg1, arg2, T_Matrix, Mrows (a), Mcols (a));
		c = D_Matrix (result);
		for (i = 1; i <= Mrows (a); i ++)
		    for (j = 1; j <= Mcols (a); j ++)
			sdata (c, i, j) = hypot (mdata (a,i,j), mdata (b,i,j));
	    } else {
		MatrixError ("hypot", a, b, M_SIZEMISMATCH, F_True);
		status = 1;
	    }
	    break;


	default:
	    type_error = F_True;
	    break;
	}
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("hypot", arg1, arg2, NULL, F_True);


    RecycleData (arg1);
    RecycleData (arg2);
    d_printf ("hypot ans=\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	log_func						*
 *									*
 * Description:	Pops and computes the natural logarithm of the		*
 *		descriptor on the top of the stack and places the	*
 *		result on the stack.  The following types are legal:	*
 *									*
 *		log (double) -> double (scalar function call)		*
 *		log (matrix) -> matrix (matrix element-wise operation)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int log_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    double	value;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;
    status = 0;


    switch (D_Type (arg)) {
    case T_Double:
	if ((value = *D_Double (arg)) <= 0.0) {
	    MathException ("logatithm of non-positive number");
	    status = 1;
	} else {
	    CreateData (result, arg, NULL, T_Double);
	    *D_Double (result) = log (*D_Double (arg));
	}
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		if ((value = mdata (a, i, j)) < 0.0) {
		    MathException ("logarithm of non-positive number");
		    status = 1;
		} else
		    sdata (b, i, j) = log (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("log", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("log ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	log10_func						*
 *									*
 * Description:	Pops and computes the base-10 logarithm of the		*
 *		descriptor on the top of the stack and places the	*
 *		result on the stack.  The following types are legal:	*
 *									*
 *		log10 (double) -> double (scalar function call)		*
 *		log10 (matrix) -> matrix (matrix element-wise op)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int log10_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    double	value;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;
    status = 0;


    switch (D_Type (arg)) {
    case T_Double:
	if ((value = *D_Double (arg)) <= 0.0) {
	    MathException ("logatithm of non-positive number");
	    status = 1;
	} else {
	    CreateData (result, arg, NULL, T_Double);
	    *D_Double (result) = log10 (*D_Double (arg));
	}
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		if ((value = mdata (a, i, j)) < 0.0) {
		    MathException ("logarithm of non-positive number");
		    status = 1;
		} else
		    sdata (b, i, j) = log10 (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("log10", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("log10 ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	pow_func						*
 *									*
 * Description:	Pops and exponentiates the two descriptors on the top	*
 *		of the stack and places the result on the stack.  The	*
 *		following types are legal:				*
 *									*
 *		pow (double, double) -> double (scalar function call)	*
 *									*
 *		An attempt is first made to coerce the arguments to	*
 *		double values.  The pow operator is used to compute the	*
 *		result, which will result in a misleading message in	*
 *		the event of an error.					*
 ************************************************************************/

int pow_func (n)
    int n;
{
    return pow_op ( );
}


/************************************************************************
 * Function:	sin_func						*
 *									*
 * Description:	Pops and computes the sine of the descriptor on the top	*
 *		of the stack and places the result on the stack.  The	*
 *		following types are legal:				*
 *									*
 *		sin (double) -> double (scalar function call)		*
 *		sin (matrix) -> matrix (matrix element-wise operation)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int sin_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	CreateData (result, arg, NULL, T_Double);
	*D_Double (result) = sin (*D_Double (arg));
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		sdata (b, i, j) = sin (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("sin", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("sin ans =\n");
    d_PrintData (result);

    return type_error == F_True;
}


/************************************************************************
 * Function:	sqrt_func						*
 *									*
 * Description:	Pops and computes the square root of the descriptor on	*
 *		the top of the stack and places the result on the	*
 *		stack.  The following types are legal:			*
 *									*
 *		sqrt (double) -> double (scalar function call)		*
 *		sqrt (matrix) -> matrix (matrix element-wise operation)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.  A cholesky factorization can be used to	*
 *		compute a true square root for matrices.		*
 ************************************************************************/

int sqrt_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    double	value;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    int		status;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;
    status = 0;


    switch (D_Type (arg)) {
    case T_Double:
	if ((value = *D_Double (arg)) < 0.0) {
	    MathException ("square root of negative number");
	    status = 1;
	} else {
	    CreateData (result, arg, NULL, T_Double);
	    *D_Double (result) = sqrt (*D_Double (arg));
	}
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		if ((value = mdata (a, i, j)) < 0.0) {
		    MathException ("square root of negative number");
		    status = 1;
		} else
		    sdata (b, i, j) = sqrt (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("sqrt", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("sqrt ans =\n");
    d_PrintData (result);

    return type_error == F_True || status != 0;
}


/************************************************************************
 * Function:	tan_func						*
 *									*
 * Description:	Pops and computes the tangent of the descriptor on the	*
 *		top of the stack and places the result on the stack.	*
 *		The following types are legal:				*
 *									*
 *		tan (double) -> double (scalar function call)		*
 *		tan (matrix) -> matrix (matrix element-wise operation)	*
 *									*
 *		An attempt is first made to coerce the argument to a	*
 *		double value.						*
 ************************************************************************/

int tan_func (n)
    int n;
{
    Matrix	a;
    Matrix	b;
    descriptor *arg;
    descriptor *result;
    descriptor	temp;
    int		type_error;
    unsigned	i;
    unsigned	j;


    result = top ( );
    temp = *result;
    arg = &temp;


    arg = CoerceData (deref (arg), T_Double);
    type_error = F_False;


    switch (D_Type (arg)) {
    case T_Double:
	CreateData (result, arg, NULL, T_Double);
	*D_Double (result) = tan (*D_Double (arg));
	break;


    case T_Matrix:
	a = D_Matrix (arg);
	CreateData (result, arg, NULL, T_Matrix, Mrows (a), Mcols (a));
	b = D_Matrix (result);
	for (i = 1; i <= Mrows (a); i ++)
	    for (j = 1; j <= Mcols (a); j ++)
		sdata (b, i, j) = tan (mdata (a, i, j));
	break;


    default:
	type_error = F_True;
	break;
    }


    if (type_error == F_True)
	TypeError ("tan", arg, NULL, NULL, F_True);

    RecycleData (arg);
    d_printf ("tan ans =\n");
    d_PrintData (result);

    return type_error == F_True;
}